Assessment of dynamic instability of laminated composite-sandwich plates

Accepted Manuscript Assessment of Dynamic Instability of Laminated Composite-Sandwich Plates

Rosalin Sahoo, B.N. Singh

PII: DOI: Reference:

S1270-9638(17)30355-3 https://doi.org/10.1016/j.ast.2018.07.041 AESCTE 4692

To appear in:

Aerospace Science and Technology

Received date: Revised date: Accepted date:

1 March 2017 7 May 2018 22 July 2018

Please cite this article in press as: R. Sahoo, B.N. Singh, Assessment of Dynamic Instability of Laminated Composite-Sandwich Plates, Aerosp. Sci. Technol. (2018), https://doi.org/10.1016/j.ast.2018.07.041

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.

Assessment of Dynamic Instability of Laminated Composite-Sandwich Plates Rosalin Sahoo* a, B. N. Singh b a

Department of Civil Engineering, Indian Institute of Technology (Banaras Hindu University) Varanasi, 221005, UP, India b

Department of Aerospace Engineering, Indian Institute of Technology Kharagpur, 721302, West Bengal, India

Abstract The current work deals with the assessment of dynamic stability behavior of laminated composite and sandwich plates subjected to in-plane static and periodic compressive loads based on a recently developed zigzag theory by the authors. This theory satisfies the traction-free boundary conditions at top and bottom surfaces of the laminate as well as the inter-laminar stress continuity at layer interfaces. Also, it obviates the need of artificial shear correction factor. The theory is based upon shear strain shape function assuming non-linear distribution of transverse shear stresses. An efficient C0 continuous, eight-noded isoparametric element with seven field variable is employed for the dynamic stability analysis of laminated composite and sandwich plates. The boundaries of principal instability domains are obtained following Bolotin’s approach and are represented either in the non-dimensional load amplitude-excitation frequency plane or load amplitude-load frequency plane. A series of numerical examples on the dynamic stability analysis of laminated composite and sandwich plates are studied to demonstrate the effects of modular ratio, span to thickness ratio, boundary conditions, thickness ratio, static load factor and various load parameters on the principal instability regions. The predicted results are compared with the available existing results in order to ensure the performance of the proposed model. Keywords: Dynamic instability; Zigzag theory; Finite element method; Laminated composite; Sandwich plate

*

Corresponding Author, email: [email protected] 1

1. Introduction Excellent specific modulus/strength properties as well as unique flexibility in design tailoring have enhanced the growth of composite materials in various weight sensitive laminated plate structures in aerospace, nuclear, marine and automotive industries. These structures are often experience to periodic in-plane loads during the service period and due to these periodic inplane loads, the plate structures may undergo unstable transverse vibrations, leading to parametric resonance. The resonance occurred due to parametric excitation, is different from the usual resonance of forced vibration as in forced vibration, resonance occurs when the frequency of the transverse forcing system coincides with the natural frequency of the structure and hence, the structure is unstable at constant frequencies of the transverse loads. However, the dynamic instability in case of parametric resonance occurs over a range of frequencies of the periodic inplane load rather than a single value which increases the complexity zone. The static instability of the structure occurs at the static buckling load values whereas the dynamic instability may occur at earlier load level than the static buckling load. Hence, the structure may not survive in the periodic in-plane loading atmosphere if it is designed to withstand static buckling load and thus, the analysis of structure subjected to periodic loading should consider the parametric resonance phenomenon. On the other hand, composite structures are weak in shear due to their low shear modulus compared to extensional rigidity and thus, the effect of shear deformation is quite significant in case of sandwich structure due to the wide variation in material properties among core and face sheet. Hence, the development of an appropriate mathematical model to analyze the structural behavior of laminated composite and sandwich plates accurately, has drawn a considerable amount of attention. The extensive reviews on development of modern plate theories are well documented [1-5]. The first order shear deformation theories (FSDT) over estimates the fundamental frequency, buckling load parameters and so on excitation frequency due to not considering the transverse shear deformation effect and assuming the transverse shear stress to be constant. Hence, FSDT requires a shear correction factor to compensate with the actual parabolic variation of the shear stresses. However, the shear correction factors are sensitive to geometrical parameter, loading and boundary conditions according to Pai [6] and thus, FSDT fails to analyze the structural behavior of plate structures accurately. 2

The further improvement comes in the form of higher order shear deformation theory (HSDT) where the need of shear correction factor has been avoided. The HSDT predicts the transverse shear strain to be continuous across the thickness, which shows discontinuity in the shear stress distribution at layer interfaces due to different values of shear rigidity at the adjacent layers. However, the actual case is just the reverse, i.e., the shear stress is continuous at the layer interfaces and the shear strain is discontinuous. The above mentioned theories belong to equivalent single layer (ESL) theory, where the deformation of the plate is expressed in terms of unknown parameters of the mid plane of the plate. However, the ESL theories are not efficient to predict accurate result for thick plates with layer in-homogeneity across the thickness. To overcome these disparities, Toledano and Murakami [7], Lu and Liu [8] and Robbins and Reddy [9] developed the layer-wise (LW) theories where every single layer is treated as a plate itself and the kinematic description is given for every single layer. The trigonometric sinus function [10] was used in LWT by Arya et al. [11], Shimpi and Ghugal [12], Shimpi and Aynapure [13] for the analysis of composite beams. Mantari et al. [14] and Roque et al. [15] modelled composite plates by trigonometric layerwise deformation theory. The finite element approach of Vidal and Polit [16], Gaudenzi et al. [17] and Pai and Palazotto [18] using mesh less method is quite appreciable. The static, buckling and vibration responses of the plate are analyzed by Ferreira et al. [19] using radial basis collocation function. Plagianakos and Saravanos [20] studied the composite and sandwich structure using higher order LWT. Though the LW theory predicts well in terms of solution accuracy, it requires a huge computational effort as the number of unknown variables increase with number of layers. The above drawbacks are overcome by the zigzag (ZZ) theory developed by Di Sciuva [21] and Murakami [22] where a zigzag concept is introduced in the in-plane displacement field and the condition of inter-laminar continuity (IC) has been satisfied at the layer interfaces. The ZZ theories have been refined further by Bhasker and Varadan [23], Di Sciuva [24] and Cho and Parmerter [25] by satisfying the zero transverse shear stress free boundary conditions at top and bottom of the plates. Furthermore, Demasi [26] considered higher order terms in the displacement field, using Murakami’s ZZ function [22]. Rodrigues et al. [27] and Neves et al. [28] also developed an approximate solution based on Murkami's ZZ theory [22] for modelling of the laminated composite plates. Moreover, Cho and Parmerter [29] included the Heaviside 3

step functions to handle the zigzag variation. Icardi [29] developed higher-order zig-zag model for analysis of thick composite beams with inclusion of transverse normal stress and sub-laminates approximations. Also, he has analysed [30] the local effects in laminated and sandwich composites using layer-wise mixed element with sub-laminates approximation and 3D zig-zag field. An extensive review on ZZ theories has been published by Carrera [31]. However, the analytical methods are restricted to some specific loading and boundary conditions and hence, Di Sciuva [32], Cho and Parmerter [33], Chakrabarti and Sheikh [34], Pandit et al. [35], Singh et al. [36], Chalak et al. [37] and Sahoo and Singh [38, 39] used finite element method (FEM) to approach the ZZ theories. On the contrary of that, many researchers have shown their major concern on various structural responses specifically the dynamic instability of structures based on different deformation theories. Bolotin [40] approached the general theory on the dynamic stability of various elastic systems. A brief review on dynamic stability analysis has well documented [4144]. The finite element (FE) dynamic instability analysis of isotropic plates based on FSDT has demonstrated by Hutt and Salam [45] and Deolasi and Datta [46]. Srinivasan and Chellapandi [47] used finite strip method (FSM) to study the dynamic instability of laminated rectangular plates subjected to uniaxial loading based on classical laminated plate theory. Moorthy et al. [48] investigated the instability of uniaxially loaded laminated composite plate based on FSDT, without considering the static component of load. Cederbaum [49] used method of multiple scales for the analysis of dynamic stability of laminated composite plate subjected to in-plane loads. Kwon [50] used FE approach and Chattopadhay and Radu [51] chosen analytical method to study the dynamic instability of laminated composite plate based on HSDT. Sahu and Datta [52] used FEM for the parametric instability of FSDT based laminated composite plates subjected to non-uniform in-plane periodic loads. Wang and Dawe [53] approached B-spline FSM for the dynamic instability analysis of composite rectangular and prismatic plate structures based on FSDT. The instability, considering the geometric non-linearity of laminated composite plate has also presented [54, 55]. The dynamic instability of anisotropic composite plates on elastic foundations has studied Patel et al. [56]. Chakrabarti and Sheikh [57] reported the FE representation of dynamic stability analysis of sandwich plate with interfacial slips based on ZZ theory. The dynamic stability analysis of laminated composite plates under harmonic axial inplane loads and with non-uniform in-plane loads has been analyzed by Fazilati and Ovesy [58] 4

and Ramachandra and Panda [59] respectively. Noh and Lee [60] studied the dynamic stability of HSDT based delaminated composite skew structures under various periodic in-plane loads. Based upon the literature review, it is observed that the developed zigzag theories where the higher order variation of in-plane displacement varies with power series expansion, predicts quite accurate results at the cost of high computational effort. Thus, there is a genuine requirement of computationally efficient ZZ plate model that needs less number of unknown variables. Also, to the author’s knowledge, no work has been reported on the FE representation of dynamic stability analysis of laminated composite and sandwich plates based on the theory that combines the zigzag concept with trigonometric shear strain functions. Keeping these viewpoints in mind, an attempt has been made in the present work to study the dynamic stability analysis of laminated composite and sandwich plates based on newly developed zigzag plate theory where the trigonometric function specifically the secant function is assumed as the shear strain shape function. The theory gives non-linear distribution of transverse shear stresses and zero transverse normal strain. Moreover, this secant function based zigzag theory (ZZTSF) satisfies the inter-laminar stress continuity at the layer interfaces as well as the traction free boundary conditions on the top and bottom surfaces of the laminated structures obviating the need of shear correction factor. The number of unknown variables considered in this model is less (seven) which makes the analysis computationally efficient. The number of unknown variables remains constant throughout the thickness. Also, the transverse displacement is assumed to be constant in the present displacement field. An efficient C0 continuous finite element is employed for the complete analysis. To avoid the usual difficulties associated with C1 continuity, the continuity requirement is reduced to C0, by assuming the derivatives of the transverse displacement as separate independent field variable. The whole domain is discretized by implementing an eight noded isoparametric biquadratic serendipity element with 56 unknown variables. The Hill’s method of infinite determinants is used to solve a system of Mathieu-type equation to obtain the dynamic instability region (DIR). Numbers of numerical examples are solved in the MATLAB environment on the dynamic stability studies of laminated composite and sandwich plates subjected to uniform in-plane periodic loads. The effects of various parameters such as span-thickness ratio, modular ratio, thickness ratio, static load factor, boundary conditions etc. on the instability region are also 5

investigated. Isotropic plate is also analyzed for the validation purposes. It is observed from the present study that, ZZTSF is more efficient in terms of mesh convergence and solution accuracy for the prediction of dynamic stability behavior of laminated composite and sandwich plates as compared to other deformation theories. 2. Mathematical formulation The basic configuration of the problem considered here is a rectangular laminated plate of dimensions (a × b × h) is subjected to uni-axial in-plane periodic loads as shown in Fig. 1. The plate is composed of finite number of linearly elastic, orthotropic layers bonded together as given in Fig. 1. A Cartesian co-ordinate system (x-y-z) is associated with the reference plane of the plate. 2.1. Displacement model Eq. (1) signifies the same displacement field as in recently developed plate model ZZTSF [39], which is assumed to be the combination of mid-plane displacements, rotations of the normal to the mid-plane, a linear zigzag and shear strain shape functions varying with layers. Heaviside step function is introduced to handle the zigzag effect. The transverse displacement is considered to be constant throughout the thickness. Fig.2 indicates the displacement configuration. The through thickness variation of the in-plane displacements are considered to be the combination of a linear zigzag model with different slopes and the parameter Ω ( k ) in each layer and the displacement model is represented as follows:

U ( x, y , z , t ) = u0 − z V ( x, y, z , t ) = v0 − z

W ( x, y, z, t ) = w0

nl −1 ∂w0 nu−1 ª º § rz · i + ¦ z − ziu H z − ziu α xu + ¦ z − z lj H − z + z lj α xlj + « z sec ¨ ¸ + Ω( k ) z » β x ∂x i =1 ©h¹ j =1 ¬ ¼

(

) (

)

(

) (

)

nl −1 ∂w0 nu−1 ª § rz · k º i + ¦ z − ziu H z − ziu α yu + ¦ z − z lj H − z + z lj α ylj + « z sec ¨ ¸ + Ω( ) z » β y ∂y i =1 ©h¹ j =1 ¬ ¼ (1)

(

) (

)

(

) (

)

where u0 , v0 , w0 are the mid-plane displacements and β x , β y are the shear deformations at the mid plane about the y- and x-axis respectively. nu and nl are number of upper and lower layers 6

i i , α yu , α xlj , α ylj are the slopes of ith and j th layer corresponding to upper and lower respectively. α xu

layers respectively and H ( z − ziu ) and H (− z + z lj ) are the Heaviside step functions. The prescribed displacement field inherently satisfies zero tangential stress-free boundary conditions at top and bottom surfaces of the plates (τ xz ,τ yz = 0 at z = ± h / 2) . The value of parameter Ω is

§ r · ª § r · § r ·º assumed to be − sec ¨ ¸ «1 + ¨ ¸ tan ¨ ¸ » for the top and bottom surfaces of the laminates. The © 2 ¹ ¬ © 2 ¹ © 2 ¹¼ i i variation of parameter Ω ( k ) , α xu , α yu , α xlj , α ylj occurs with the change in layers and is evaluated by

implementing the inter-laminar continuity conditions (τ xzi = τ xzi +1 and τ yzi = τ yzi +1 ) at the layer interfaces. The optimized value of transverse shear stress parameter r is ascertained by the inverse method [61] in the post processing step and is chosen to be 0.1 [62]. The actual displacement field shown in Eq. (1) contains the first order derivatives of the transverse displacement which requires C1 continuity at the element interfaces during its finite element implementation and dealing with the C1 continuity is computationally expensive. Thus, the derivatives of the transverse displacement are considered as separate field variables

§ § ∂w0 · § ∂w · · , θ y ¨ = 0 ¸ ¸ in the modified displacement field to avoid the usual complexities ¨θ x ¨ = ¸ © ∂y ¹ ¹ © © ∂x ¹ associated with C1 continuity. Thus, C1 continuity with five DOF is transformed into C0 continuity with seven DOF and the modified displacement field can be presented as follows: nu −1

(

) (

)

nl −1

(

) (

)

)

nl −1

(

) (

)

i + ¦ z − z lj H − z + z lj α xlj + ª¬ g ( z ) + Ω( k ) z º¼ β x U ( x, y, z, t ) = u0 − θ x + ¦ z − ziu H z − ziu α xu i =1 j =1 nu −1

(

) (

i + ¦ z − z lj H − z + z lj α ylj + ª¬ g ( z ) + Ω( k ) z º¼ β y V ( x, y, z, t ) = v0 − θ y + ¦ z − ziu H z − ziu α yu i =1 j =1

W ( x, y, z, t ) = w0

(2)

The considered value of the function g ( z ) is z sec ( rz / h ) . The substituted artificial constraints in Eq. (2) are imposed variationaly incorporating the following constraint equation using penalty approach [63]. 7

∂w0 − θ x = 0; ∂x

∂w0 −θy = 0 ∂y

(3)

2.2. Constitutive equations The stress-strain relation of k th layer of an orthotropic lamina with respect to structural axis system (x-y-z) may be expressed as

{σ }

ij 5×1

= ª¬Qij º¼

(k )

{ε } ij

(4) 5×1

{ }

where σ ij

= ª¬σ xx σ yy τ xy τ yz τ xz º¼ , {ε ij } = ª¬ε xx ε yy γ xy γ yz γ xz º¼ and ª¬Qij º¼ is (k ) 5×1 5×1 T

T

the reduced transformed rigidity matrix of k th layer of multi-layered plate. The rigidity matrix

ª¬Qij º¼ can be evaluated in terms of material properties and fibre orientation α . (k ) 3.

Finite element method

3.1. Plate element An eight noded isoparametric biquadratic quadrilateral serendipity element (Fig. 3) with

(

(

)

)

seven DOF u0 , v0 , w0 , θ x ( = w0, x ) , θ y = w0, y , β x , β y per node is employed for the discretization of the plate domain. However, the unknown field variables and element geometry can be expressed in terms of the shape function N i [55] associated with the node i , as follows: n

n

n

i =1

i =1

i =1

δ = ¦ Niδ i ; x = ¦ Ni xi ; y = ¦ Ni yi

(5)

where δ is the generalized field variable, x, y are the generalized Cartesian coordinate and δi ,

xi , yi are the corresponding field variable and coordinate value at ith node respectively and n is the number of nodes per element.

8

3.2. Linear strain-displacement relations

The linear strain-displacement relationships are presented as

ª ∂U {ε } = « ¬ ∂x

∂V ∂y

∂U ∂V + ∂y ∂x

∂V ∂W + ∂z ∂y

∂U ∂W º + ∂z ∂x »¼

T

(6)

Further, the linear strains given in Eq. (6) can be represented in terms of generalized strains which are the functions of generalized unknown field variables.

{ε } = [ H ]5 x14 {ε }14 x1

(7)

Where, the elements of [ H ] are the functions of z , unit step functions and the parameter

Ω which

are shown in Appendix A.

{

Further, {ε } = ε10

ε 20 ε 60 k11 k21 k61 k12 k23 k64 k65 ε 40 ε 50 k46 k57 }

T

(8)

Moreover, the strain vector {ε } provided in Eq. (6) may be presented in terms of unknown field variables as

{ε } = [ B]{δ }

(9)

where, [ B] is the strain-displacement matrix in the Cartesian coordinate system and {δ } is the displacement vector consists of the fundamental DOF and thus it may be reported as

{δ } = {u0

v0

w0 θ x θ y

βx

βy}

T

(10)

3.3. Geometric strain-displacement relations

The generalized geometric strain vector may be expressed as

9

ª 1 2 1 2 1 2 º « 2 W, x + 2 U , x + 2 V, x » « » 1 2 1 2 1 2 » 1 « {ε G } = « W, y + U, y + V, y » = [ AG ]{θG } 2 2 2 2 «W W + U U + V V » ,x , y ,x , y » « ,x , y «¬ »¼

ªW, x 0 U , x 0 V, x 0 º « » where [ AG ] = « 0 W, y 0 U , y 0 V, y » «W, y W, x U , y U , x V, y V, x » ¬ ¼

{θG } = ª¬W, x

W, y U , x U , y V, x V, y º¼

(11)

(12)

T

(13)

where, comma (,) shows the partial differentiation with respect to the coordinate subscripts that follow. Again, the above equation may be expressed with the use of Eq. (9), as

{θG } = [ H G ]{ε } = [ H G ][ BG ]{δ }

(14)

where, [ H G ] is also a function of z , unit step functions and the parameter Ω like that of [ H ] which are indicated in Appendix A.

3.4. Energy Expressions

The total potential energy Π of the system may be presented as

∏ = T − U S − U in − Pλ

(15)

where T is the total kinetic energy, U S is the strain energy, U in is the potential energy due to applied in-plane compressive loads and Pλ shows the penalty term.

10

3.4.1. Kinetic energy of laminate The kinetic energy (T ) of the laminated plate may be presented as T=

{ } ρ {L } dv = 12 ³³ {L} [ I ]{L}dxdy

1 n L ¦ ³ 2 k =1 v

T

T

k

(16)

where ρ k is the mass density of the k th layer and the dot (.) shows the derivative with respect to time ' t '. The displacement vector may be expressed in terms of displacement vector using Eq. (5), as

{L} = [Q ]{δ }

(17)

3.4.2. Strain energy of laminate The strain energy of composite plate is expressed by

US =

1 n T {ε } [σ ] dv ¦ ³ 2 k =1 v

(18)

The strain energy can be presented using Eqs. (2) and (6), as

US =

1 n 1 n T T ª º ε Q ε dxdydz = { } ¬ ¼(k ) { } {ε } [ D ]{ε } dxdy ¦ ¦ ³³³ ³³ 2 k =1 2 k =1

where [ D ] =

(19)

1 n T [ H ] ª¬Q º¼ ( k ) [ H ] dz ¦ ³ 2 k =1

3.4.3. Potential energy due to applied in-plane loads

The potential energy due to applied in-plane compressive loads with the use of Eqs. (11)- (14), is presented as 11

1 n 1 n T T k ª º U in = ¦ ³³³ [θG ] ¬ S ¼ [θG ] dxdydz = ¦ ³³ {ε } [G ]{ε } dxdy 2 k =1 2 k =1 n

(20)

where [G ] = ¦ ³ [ H G ] ª¬ S k º¼ [ H G ] dz T

k =1

The stress matrix ª¬ S k º¼ may be expressed in terms of in-plane stress components

(σ11 , σ 22 , τ12 ) of the k th

layer as

0 0 0 0 º ªσ 11 τ 12 «τ 0 0 0 »» « 12 σ 22 0 «0 0 σ 11 τ 12 0 0 » ª¬ S k º¼ = « » 0 τ 12 σ 22 0 0 » «0 «0 0 0 0 σ 11 τ 12 » « » 0 0 0 τ 12 σ 22 ¼» ¬« 0 3.4.4. Penalty approach

The penalty term is expressed in terms of penalty parameter, ‘ γ ’ (taken as 106 in the present work) and field variables associated with artificial constraints using Eq. (3), is represented as T T ª§ ∂w · § ∂w0 ·º · § ∂w0 · § ∂w0 0 −θx ¸ ¨ −θx ¸ + ¨ −θ y ¸ ¨ − θ y ¸ »dxdy Pλ = ³³ «¨ 2 «© ∂x ¹ © ∂x ¹ © ∂y ¹ © ∂y ¹ »¼ ¬

γ

(21)

Eq. (21) can be represented using Eq. (6), as follows:

­ ∂w ½ ­ ∂w0 ½ − θ x ¾ = Pxδ ; ® 0 − θ y ¾ = Pyδ ® ¯ ∂x ¿ ¯ ∂y ¿

(22)

The total elemental potential energy as shown in Eq. (15) may be rewritten considering the above equations, as

∏e =

1  δ 2

{ } [ M ]{δ} − 12 {δ } [ K ]{δ } − 12 {δ } [ K ]{δ } − 12 {δ } T

T

e

T

e

T

Ge

ª¬ K pe º¼ {δ }

(23) 12

where

[ M e ] = ³³ [Q] [ I ][Q] dxdy; [ Ke ] = ³³ [ B] [ D][ B] dxdy T

T

[ KGe ] = ³³ [ BG ] [G][ BG ] dxdy; T

(

)

ª¬ K pe º¼ = γ ³³ PxT Px + PyT Py dxdy

The element stiffness matrix [ K e ] and element penalty matrix ª¬ K pe º¼ are summed up together for each element. The [ K Ge ] and [ M e ] are the element geometric stiffness matrix and element mass matrix respectively. These matrices are accordingly assembled together to form the corresponding global matrices and are reported as the corresponding terms without subscripts. 3.5. Governing equations

The Lagrange equation for the conservative system yields to the equations of motion of a structure under in-plane load with the global mass matrix [ M ] , global elastic matrix [ K ] and global geometric stiffness matrix [ K G ] which may be expressed as

[ M ]{δ} + ª¬[ K ] − N ( t ) [ KG ]º¼ {δ } = 0

(24)

Eq. (24) is a general governing equation which can be reduced to have the governing equations for the eigen value problem of buckling, vibration, and dynamic stability of the plate structure.

3.5.1. Buckling

ª¬[ K ] − Ncr [ KG ]º¼ {δ } = {0}

(25)

where N cr is the critical buckling load.

3.5.2. Vibration 13

ª¬[ K ] − N [ KG ]º¼ {δ } − ω2 [ M ]{δ } = {0}

(26)

In Eq. (26), ω become the vibration frequency when the structure is subjected to in-plane load. However, when the structure is considered without in-plane load i.e., N attains zero value, the ω become the natural frequency of vibration. The frequency of vibration becomes zero when

N = N cr .

3.5.3. Dynamic stability studies The Eq. (24) can be used as the governing equation for the dynamic stability analysis of a plate structure subjected to uniform in-plane periodic load N . Also, the periodic load N can be expressed as N ( t ) = N 0 + N 1 cos Ω e t

(27)

where Ωe is the excitation frequency, N0 is the static portion of N ( t ) and N1 is the amplitude of the dynamic portion of N ( t ) . The periodic load N can also be presented in terms of the linear static buckling load Ncr as N ( t ) = α N cr + β N cr cos Ω e t

(28)

where α = N0 / Ncr and β = N1 / N cr are denoted as the static and dynamic load factors respectively. The equation of motion may be expressed by substituting Eq. (28) into Eq. (24), as follows:

[ M ]{δ} + ª¬[ K ] −α Ncr [ KG ] − β Ncr [ KG ] cos Ωet º¼ {δ } = 0

(29)

Eq. (29) is a set of Mathieu type equation, governing the instability behavior of the plate structure. The solution of this equation may be either bounded or unbounded for the given values of three parameters α , β and Ωe . The spectrum of these values of parameters has unbounded 14

solutions for some regions of the planes because of the parametrically excited resonance. This phenomenon is termed as dynamic instability and these regions are known as dynamic instability regions (DIR). The boundaries of DIRs are determined using Bolotin’s approach [31] by the periodic solutions having periods T and 2T . The practical importance of the instability region at the boundaries of period 2T has high impact [31] and thus, the solution of Eq. (29) is determined by expressing the components {δ } in trigonometric series form as follows:

{δ } =

iΩ t iΩ t º ª ¦ «¬{a } sin 2 + {b } cos 2 »¼ ∞

i

e

i

e

with period 2T where T =

i =1,3,5..

2π , or Ωe

∞ iΩ t iΩ t 1 {δ } = {b}0 + ¦ ª«{ai } sin e + {bi } cos e º» with period T . 2 2 2 ¼ i = 2,4,6.. ¬

(30)

(31)

Substituting the above equations in Eq. (29) and equating the sums of the coefficients of each sine and cosine terms to zero, a series of algebraic equations are achieved for the determination of instability regions. The principal instability region having great practical significance [31] is correspond to i = 1 and hence, the dynamic instability equation for this case leads to

ª º Ωe2 1 K α N K β N K − ± − [ M ]» {δ } = 0 cr [ G ] cr [ G ] «[ ] 2 4 ¬ ¼

(32)

Basically, Eq. (32) is a generalized eigenvalue problem of the systems for the known values of

α , β and Ncr . The two boundaries of the DIR are indicated by two conditions under a plus and minus sign. The eigenvalue problem is solved and the excitation frequencies Ωe are determined for the known values of α and β . 4. Numerical examples To verify the accuracy of the present C0 continuous 2D FE plate model, a number of numerical problems on laminated composite and sandwich plates subjected to uniform in-plane 15

edge loading are solved. The evaluated results with the compared published results based on various deformation theories are presented in the form of different tables and figures. To avoid the shear locking phenomena, the Gauss quadrature rule with reduced integration scheme (i.e., the 3 ×3 and 2 ×2 Gauss points are used to integrate the bending and the shear terms respectively) for thin laminate whereas full integration rule are adopted for thick laminate. The dynamic instability analyzes are investigated with various combinations of the edge support such as: (i) Simply supported boundary SSSS: v = w = θ y = β y = 0 at x = 0, a and u = w = θ x = β x = 0 at y = 0, b (ii) Clamped boundary CCCC: u = v = w = θ x = θ y = β x = β y = 0 at x = 0, a and y = 0, b Example 1: Symmetric cross-ply laminated composite plates: Vibration and Buckling Analyses A simply supported symmetric cross-ply laminated composite plate (0/90/90/0) is analysed using the material properties ( E11 = Open, E22 = E22 , E33 = E22 , G12 = 0.6E22 , G13 = 0.6E22 ,

G23 = 0.5E22 , υ12 = 0.25, ρ = 1 Kg/m3). The non-dimensional frequencies ω = ωb 2 / h ρ / E22 for thick plate ( a / h = 5) are obtained for different modular ratios ( E11 / E22 = 10, 20, 30, 40) as given in Table 1. The table shows that, the present FE results agreed very well with the 3D elasticity solution (Noor [65]) and also with the ZZ results of Chalak et al. [66] and Rodrigues et al. [67], mesh free results of Liew et al. [68] and HSDT results of Phan and Reddy [71]. Notably, the present results are in good agreement with the elasticity results and the results based on various deformation theories. Also, the gradual increase in modular ratio increases the fundamental frequency of the laminated plate. Again the average percentage error from exact solution [65] is around 0.47% as compared to 0.61% of Chalak et al. [66], 1.56% of Rodrigues et al. [68] and 0.59% of Liew et al. [69]. Further, the non dimensional buckling load parameters λ = λ b 2 / E22 h 3 of simply supported uni-axially loaded three-layered laminated thick plate ( a / h = 10 ) with the variation of modular ratio ( E11 / E22 = 10, 20,30, 40 ) are evaluated and tabulated in Table 2. The evaluated results are also compared with the elasticity results [70], HSDT results 16

(Reddy and Phan [71]), local higher order theory (LHT) results (Wu et al. [72]), global local higher order theory (GLHT) results (Lo et al. [73]), ZZ results (Singh et al. [74]). It can be observed that the present results are in well agreement with the exact solution [70] and the results obtained by other deformation theories specifically results obtained by ZZ theories [74]. However, the ESL theories over predict the buckling loads for laminated composite and sandwich plates as expected. The performance of the present model is quite remarkable as it predicts better result assuming only seven number of field variables whereas Singh et al. [74] obtained the ZZ results by considering eleven number of unknowns making the analysis more complex. Also, the non-dimensional buckling loads increases significantly with the increasing modular ratio. The present formulation in terms of solution accuracy in the direction of higher degree of material anisotropy is excellent which assures the flexibility of the present model. Example 2: Simply supported isotropic plate: Dynamic Stability Analysis

A simply supported square isotropic plate subjected to uni-axial in-plane uniform edge loading is analyzed in this example. The following material properties are used for this problem,

E11 = E22 , G23 = G13 = G12 = E22 / 2.5 , υ12 = 0.3, ρ = 1.0 kg/m3. The non-dimensional excitation

( )

( )

frequencies Ω e = Ω e a 2 ρ h / D for the upper ΩUe and lower ΩeL bound region are evaluated for thin plate ( a / h = 100) using the mesh sizes of 2 × 2, 4 × 4, 6 × 6, 8 × 8 and 10 × 10. The excitation frequencies are obtained with the variation of the static load factor (α ) from 0.0-0.35 and dynamic load factor ( β ) from 0.2-1.2 and are presented in Table 3. It is observed that, the excitation frequencies are converged well for mesh division 10 × 10 and hence, subsequent analyzes are carried out with this mesh division. Again, the obtained results in the form of excitation frequencies are compared with the ZZ results of Chakrabarti and Sheikh [57],

17

analytical results of Hutt and Salam [45] and FEM results of Srivastava et al. [75] based on FSDT. Well agreements of results are observed among the present theory ZZTSF and the compared ones (Chakrabarti and Sheikh [57], Hutt and Salam [45], and Srivastava et al. [75]). Again, it is clear from the table that, the upper bound non-dimensional excitation frequencies increase with the increase in dynamic load factor and decrease in static load factor whereas the lower bound non-dimensional excitation frequencies reduce with the reduction in static and dynamic load component. Further, the non-dimensional excitation frequencies are obtained with the variation in dynamic load component and the instability regions are plotted in Fig. 4 for the discussed thin isotropic plate subjected to dynamic in-plane load. The static load component is neglected for this analysis i.e., (α = 0.0) . The figured plot of instability region is also validated with the DIR plot of Hutt and Salam [45]. It is notified that the DIR based on ZZTSF is well agreed with the corresponding figure of Hutt and Salam [45]. Again, increased dynamic load factors widen the dynamic stability opening (DIO) and the DIO is the distance between the two boundaries at a particular load level.

Example 3: Four layered symmetric cross-ply laminated composite plate

A symmetric cross-ply square laminate with equal thickness plies (0/90/90/0) subjected to uni-axial periodic load under simply supported edge conditions is studied in this example. The material properties used for this computation are as follows (Reddy [76]): E11 / E22 = 40,

G12 / E22 = G13 / E22 = 0.6, G23 / E22 = 0.5, υ12 = 0.25, ρ = 1.0 kg/m3. The non-dimensional excitation frequencies Ω e = Ω e ( a 2 / h ) ρ / E22 are evaluated for both upper and lower 18

boundaries varying the static (α ) and dynamic load factor ( β ) assuming the span-thickness ratio

( a / h) as 25 and are tabulated in Table 4. The initial part of the table contains the upper and lower bound resonance frequencies which are evaluated with the presence of only the dynamic load component and ignoring the static load factor. It is observed that, the gradual increase in dynamic load component ( β = 0.0 − 1.5 ) increases the upper bound excitation frequencies and reduces the lower bound values of excitation frequencies. The final part of the table carries the evaluated upper and lower bound values of excitation frequencies with the variation of both static and dynamic load components. The static load component varies from 0.2 to 0.8 and the dynamic load component varies from 0.06 to 0.24 assuming the ratio (α / β ) to be 0.3. The table illustrates that the increase in static and dynamic load components reduces both upper and lower bound values of excitation frequencies. It shows that, the presence of the compressive static load component makes the reduction in the stiffness of the laminate. The table also includes the published results of Wang and Dawe [54] and Kao et al. [77] along with the evaluated excitation frequencies for the validation purposes. The evaluated results based on ZZTSF are well agreed with the compared results ([54], [77]) which is clearly visible. Further, the instability regions are plotted for the same case in the plane of ( Ωe / ω , β ) in the Fig. 5. The load frequency ( Ω e / ω ) i.e., ratio of excitation frequency and computed fundamental frequency of vibration is considered as the ordinate and the load amplitude

( β = N1 / N cr ) as the abscissa. The figure also includes the published DIR plot of Moorthy et al. [39] for comparison. It may be notified that the present results are well agreed with the reference (Moorthy et al. [48]). The primary instability region occurs in the vicinity of Ωe / ω = 2.0 when the static and dynamic load components are ignored. The region with vertices at Ωe = 2ω is called as the principal instability region which is of great practical importance. Also, the width of

19

instability zones increases gradually with the increase in dynamic load component. This case is defined as the standard one for all further computations in line with Moorthy et al. [39]). Moreover, the discussed laminate is analyzed again to study the effect of modular ratio

( E11 / E22 ) on

the dynamic instability assuming the span-thickness ratio to be 25. The load

frequencies ( Ω e / ω ) with respect to the dynamic load factor are evaluated for different values of modular ratio ( E11 / E22 ) and are plotted in Fig. 6. Here, the static load components are neglected and only dynamic load components are considered that varies from 0.0 to 1.0. The plotted instability region is also compared with the published figure of Moorthy et al. [48]. It is illustrated from the figure that the evaluated results are in well agreement with the established results of Moorthy et al. [48]. Also, the reducing modular ratio makes the onset of instability to occur earlier with a wider instability opening ( a / h = 25 ) . The non-dimensional excitation frequencies are evaluated corresponding to the variation in dynamic load factor assuming the static load component as 0.4 for different boundaries (SSSS, SSSC and CCCC) and are plotted in Fig. 8. It is clearly visible from the figure that the onset of instability occurs at a higher excitation frequency with wider instability opening moving from SSSS to CCCC boundaries due to the increase of restraint at the edges. Again, the influence of the static load factor on the dynamic instability of the same thick plate ( a / h = 10 ) is analyzed. The non-dimensional excitation frequencies with respect to the variation of dynamic load factor for different values of static component of load are evaluated and plotted in Fig. 7. It is notable from the figure that the onset of dynamic instability occurs earlier with wider instability region for the increasing static component of load. Example 4: Five layered un-symmetric sandwich plate An un-symmetric simply supported sandwich plate (0/90/C/0/90) subjected to uni-axial periodic load is considered for the dynamic instability analysis. The plate consists of an isotropic core layer with an identical thickened un-symmetric cross-ply stiff face sheets. The material properties used for the face sheets are E11 = 131GPa,

E22 = 10.34GPa, G12 = 6.895GPa,

υ12 = 0.22, ρ = 1627 Kg/m3 and same for the core is E11 = 6.89 ×10−3 GPa, E22 = E11 , 20

G12 = 3.45 × 10 −3 GPa, υ12 = 0.0,

ρ = 97 Kg/m3 (Kant and Swaminathan [78]). The span-

(

thickness ratio ( a / h ) as well as the core-to-face thickness ratio hc / hf

) is assumed to be 10 for

the computation. The load frequency is evaluated with respect to the dynamic load factor and the DIRs are plotted in ( Ωe / ω , β ) plane for different values of static component of load in Fig. 9. It is clearly observed from the figure that the primary instability region occurs in the vicinity of Ωe / ω = 2 for the zero values of static and dynamic load factor as expected. Also, the onset of instability regions occurs earlier with wider instability opening due to the increase of compressive static in-plane load. Further, the effects of core-to-face thickness ratio ( hc / h f ) on the dynamic instability are studied for the same thick plate ( a / h = 10 ) . The upper and lower bound non-dimensional excitation frequencies Ω e = Ω e ( b 2 / h ) ρ f / E22 f with respect to the various thickness ratios

( h / h ) ranging from 4 to 50 are obtained and are plotted in Fig. 10. The figure illustrates that c

f

the gradual increase in hc / hf ratio increases the non-dimensional excitation frequencies with narrow instability opening zone. Further, to verify the accuracy of the results, a simple comparison among the evaluated non-dimensional natural frequency ω = ω ( b 2 / h ) ρ f / E22 f and non-dimensional excitation frequency are studied for the same sandwich plate for the thickness ratios ( hc / h f ) as of 4, 20 and 100 according to Table 5. The present computation neglects the static and dynamic load factors such as (α = β = 0.0) for the comparison purposes, as the established results are insufficient to compare both natural and excitation frequencies of sandwich plate based on ESL as well as ZZ theory. The table clearly signifies that, the non-dimensional excitation frequencies are twice of the non-dimensional natural frequencies with the non-consideration of static and dynamic load components. The table also includes the published results of Kao et al. [77] based on ESL theory and the available results of Rao et al. [78] based on refined theory. It is quite well known that the ESL theories over estimate the value of natural frequencies due to the ignorance of inter-laminar continuity effect in the assumed displacement field (Kao et al. [77]) and so as it predicts erroneous results for the excitation frequencies as obvious. However, ZZTSF predicts quite 21

better result for natural frequencies as it is well agreed with the refined theory results (Rao et al. [78]) for a certain range of property variation. Accordingly, the evaluated excitation frequencies using present approach are quite accurate as the introduction of static and dynamic load components keeps the result in the same range. This is the reason behind the analysis of dynamic instability of sandwich plate based on computationally effective ZZTSF which requires only seven numbers of unknown variables. Again, the same plate is analyzed to capture the influences of the span-thickness ratio on the considering the thickness ratio of the core thickness ( hc / h f ) as 10. The non-dimensional excitation frequencies Ω e = Ω e ( b 2 / h ) ρ f / E22 f are evaluated for a wide range of span thickness ratios ( a / h = 100,50, 40, 20,10 ) such as for thin, thick and moderately thick plates and the DIRs are plotted in Fig. 11. It can be observed from the figure that, the onset of instability occurs later with wider instability region for the increasing span-thickness ratio.

Example 5: Five layered symmetric sandwich plate This example deals with a five layered uni-axially loaded square sandwich plate with symmetric cross-ply faces (0/90/C/90/0). The ratio of thickness of core to the total thickness

( hc / h ) is considered as 0.8 and each of the plies of face sheets are assumed to be of identical thickness. The material properties of the faceplates are E11 = 276 GPa, E22 = G12 = G13 = G23 = 10.34 GPa, υ12 = 0.22, ρ = 681.8 Kg/m3and those for the core are E11 = E22 = 0.5776 GPa,

G12 = G13 = 0.1079 GPa, G23 = 0.22215GPa, υ12 = 0.0025, ρ = 1000 Kg/m3. The nondimensional excitation frequencies Ω e = 100Ω e a ρ c / E11 f are evaluated for both lower and upper bound with the variation of dynamic load factor ( β ) for different boundaries i.e., SSSS, SSSC, SSCC and CCCC and are plotted in Fig. 12. The figure signifies that the onset of instability occurs later with higher values of excitation frequencies as well as with wider instability regions due to the introduction of restraint at the edges in the clamped plate. 5. Conclusion 22

The dynamic stability analysis of the laminated composite and sandwich plates based on a recently developed zigzag theory by the authors is studied to ensure the theory’s practicability. The theory satisfies the inter-laminar continuity at the layer interfaces as well as the shear-stressfree boundary conditions on the plate boundary surface obviating the need of shear correction factor. The theory assumes the realistic non-linear distribution of transverse shear stresses. An efficient C0 finite element is implemented to solve a number of numerical problems on the dynamic stability analysis of laminated composite and sandwich plates featured with various boundary conditions, material properties, thickness ratios and modular ratios. Some important observations are drawn from the whole analyses which are enlisted as given below: 1. The increasing value of dynamic load component makes the dynamic instability opening wider. 2. The onset of instability occurs earlier with wider zone of instability opening for the gradual increase in modular ratio. 3. The onset of instability regions tend to shift to lower frequencies with the increasing static component of the compressive in-plane load, perusing a destabilizing effect on the dynamic stability analysis. 4. The onset of instability occurs later with wider zone of instability due to the introduction of restraint at the edges. 5. The reduction in core thickness ratios minimizes the value of excitation frequencies with wider zone of instability opening. 6. The increasing span-thickness ratio increases the excitation frequencies gradually with wider instability region i.e., the thick plate predicts a narrow instability zone showing more stiffness. The accuracy of the present model is verified by comparing the evaluated results with the available established results. The ZZTSF is approved to be an appropriate model in terms of

23

excellent mesh convergence, well agreed results with the compared ones and accurate prediction of dynamic instability behavior.

Appendix A:

ª1 «0 « [ H] = «0 « «0 «¬0

0 0 z 0 0 p1 0 1 0 0 z 0 0 p2

0 0

0 1 0 0 z 0 0 0 0 0 0 0 0 0 0 0 0 0

p3 0 0

0 0 0

0º 0 »» p4 0 0 0 0 » where » 0 1 0 q1 0 » 0 0 1 0 q2 »¼ 0 0 0 0 0 0 0 0

nu −1

nl −1

i =1

j =1

p1 = g ( z ) + ¦ (Ω − Ωux )( z − ziu ) H ( z − ziu ) + ¦ (Ω − Ωlx )( z − z lj ) H ( − z + z lj ) nu −1

nl −1

p2 = g ( z ) + ¦ (Ω − Ω )( z − z ) H ( z − z ) + ¦ (Ω − Ωly )( z − z lj ) H ( − z + z lj ) u y

u i

u i

i =1

nu −1

j =1

nl −1

q1 = g '( z ) + ¦ (Ω − Ωuy ) H ( z − ziu ) + ¦ (Ω − Ω ly ) H ( − z + z lj ) i =1

j =1

nu −1

nl −1

i =1

j =1

q2 = g '( z ) + ¦ (Ω − Ωux ) H ( z − ziu ) + ¦ (Ω − Ωlx ) H ( − z + z lj )

ª1 «0 « «0 [ HG ] = «0 « «0 « ¬«0

0 0 0 z 0 0 0 p1 0

0

1 0 0 0 z 0 0 0 p1 0 0 1 0 0 0 z 0 0 0 p2 0 0 1 0 0 0 z 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0

0 0

0 0 0º 0 0 0»» 0 0 0» » p2 0 0» 0 1 0» » 0 0 1¼»

References

24

[1] E. Reissner, Reflections on the theory of elastic plates, Appl. Mech. Rev. 38 (1985) 14531464. [2] A. K. Noor, W. S. Burton, Refinement of higher-order laminated plate theories, Appl. Mech. Rev. 42(1) (1989) 1–13. [3] R. K. Kapania, S. Raciti, Recent advances in analysis of laminated beams and plates, part ii: vibrations and wave propagation. AIAA J. 27(7) (1989) 935–946. [4] J. N. Reddy, A review of refined theories of laminated composite plates, Shock Vib. Dig. 22(7) (1990) 3–17. [5] Mallikarjuna, T. Kant, A critical review and some results of recently developed refined theories of fiber-reinforced laminated composites and sandwiches, Compos. Struct. 23(4) (1993) 293–312. [6] P. F. Pai, A new look at shear correction factors and warping functions of anisotropic laminates, Int. J. Solids Struct. 32(16) (1995) 2295–2313. [7] A. Toledano, H. Murakami, Composite plate theory for arbitrary laminate configuration, J. Appl. Mech. ASME 54 (1987) 181-189. [8] X. Lu, D. Liu, An inter-laminar shear stress continuity theory for both thin and thick laminates, ASME J. Appl. Mech. 59 (1992) 502-509. [9] D. H. Jr Robbins, J. N. Reddy, Modeling of thick composites using a layer-wise laminate theory, Int. J. Numer. Meth. Eng. 36 (1993) 655–677. [10] M. Touratier, An efficient standard plate theory, Int. J. Eng. Sci., 29 (8) (1991), pp. 745-752 [11] H. Arya, R.P. Shimpi, N.K. Naik, A zig–zag model for laminated composite beams, Compos. Struct., 56 (2002), pp. 21-24 [12] R.P. Shimpi, A.V. Aynapure, A beam finite element based on layerwise trigonometric shear deformation theory, Compos. Struct., 53 (2001), pp. 153-162 [13]R.P. Shimpi, Y.M. Ghugal, A layerwise trigonometric shear deformation theory for two layered cross-ply laminated beams, J. Reinf. Plast. Compos., 18 (1999), pp. 1516-1542 25

[14] J. L. Mantari, A. S. Oktem, C. G. Soares, A new trigonometric layerwise shear deformation theory for the finite element analysis of laminated composite and sandwich plates, Comput. Struct. 94-95 (2012) 45–53. [15] C. M. C. Roque, A. J. M. Ferreira, R. M. N. Jorge, Modelling of composite and sandwich plates by a trigonometric layerwise deformation theory and radial basis functions, Compos.: Part B 36 (2005) 559–572. [16] P. Vidal, O. Polit, A family of sinus finite elements for the analysis of rectangular laminated beams, Compos. Struct., 84 (2008), pp. 56-72 [17] P. Gaudenzi, R. Barboni, A. Mannini, A finite element evaluation of single-layer and multi-layer theories for the analysis of laminated plates, Comput. Struct., 30 (1995), pp. 427-440 [18] P.F. Pai, A.N. Palazotto, A higher-order sandwich plate theory accounting for 3-D stresses, Int. J. Solids Struct., 38 (2001), pp. 5045-5062. [19] A. J. M. Ferreira, C. M. C. Roque, E. Carrera, M. Cinefra, O. Polit, Two higher order ZigZag theories for the accurate analysis of bending, vibration and buckling response of laminated plates by radial basis functions collocation and a unified formulation, J. Compos. Mat. 45(24) (2011) 2523–2536. [20] T.S. Plagianakos, D.A. Saravanos, High-order layerwise laminate theory for the prediction of interlaminar shear stresses in thick composite and sandwich composite plates Compos. Struct., 87 (2009), pp. 23-35 [21] M. di Sciuva, Bending, vibration and buckling of simply supported thick multilayered orthotropic plates: an evaluation of new displacement model, J. Sound Vib. 105(3) (1986) 425– 444. [22] H. Murakami, Laminated composite plate theory with improved plate theory with improved in-plane responses, J. Appl. Mech. 53 (1986) 661–6. [23] K. Bhasker, T. K. Varadan, Refinement of higher order laminated plate theories, AIAA 27(12) (1989) 1830–1831. [24] M. di Sciuva, Multilayered anisotropic plate models with continuous inter-laminar stress Comput. Struct. 22(3) (1992) 149-167. 26

[25] M. Cho, R. R. Parmerter, Efficient higher order composite plate theory for general lamination configurations, AIAA 31(7) (1993) 1299–1306. [26] L. Demasi, Refined multilayered plate element based on Murakami's zig-zag function. Compos. Struct. 70 (2005) 308–316. [27] J. D. Rodrigues, C. M. C. Roque, A. J. M. Ferreira, E. Carrera, M. Cinefra, Radial basis functions-finite differences collocation and a Unified Formulation for bending, vibration and buckling analysis of laminated plates, according to Murakami's zig-zag theory, Compos. Struct. 93(7) (2011) 1613–1620. [28] A. M. A. Neves, A. J. M. Ferreira, E. Carrera, M. Cinefra, C. M. C. Roque, R. M. N. Jorge, Static analysis of functionally graded sandwich plates according to a hyperbolic theory considering zig-zag and warping effects, Adv. Eng. Soft. 52 (2012) 30–43. [29] U. Icardi, Higher-order zig-zag model for analysis of thick composite beams with inclusion of transverse normal stress and sublaminates approximations, Comp Part B, 32 (2001), pp. 343-354. [30] U. Icardi, Layerwise mixed element with sublaminates approximation and 3D zig-zag field, for analysis of local effects in laminated and sandwich composites, Int J Numer Methods Eng, 70 (2007), pp. 94-125. [31] E. Carrera, Historical review of Zig-Zag theories for multilayered plates and shells, Appl. Mech. Rev. 56 (2003) 1-22. [32] M. di Sciuva, A third order triangular multilayered plate finite element with continuous inter-laminar stresses, Int. J. Numer. Meth. Eng. 38 (1995) 1–26. [33] M. Cho, R. R. Parmerter, Finite element for composite plate bending based on efficient higher order theory, AIAA 32(11) (1993) 2241–2245. [34] A. Chakrabarti, A. H. Sheikh, A new triangular element to model inter-laminar shear stress continuous plate theory, Int. J. Numer. Meth. Eng. 60 (2004) 1237–1257.

27

[35] M. K. Pandit, A. H. Sheikh, B. N. Singh, An improved higher order zigzag theory for the static analysis of laminated sandwich plate with soft core, Finite Ele. Anal. Des. 44(9–10) (2008) 602–610. [36] S. Singh, A. Chakrabarti, P. Bera, J. S. D. Sony, An efficient C0 FE model for the analysis of composites and sandwich laminates with general layup, Lat. Amer. J. Solid Struct. 8 (2011) 197–212. [37] H. D. Chalak, A. Chakrabarti, M. A. Iqbal, A. H. Sheikh, An improved C0 FE model for the analysis of laminated sandwich plate with soft core, Finite Ele. Anal. Des. 56 (2012) 20–31. [38] R. Sahoo, B. N. Singh, A new trigonometric zigzag theory for buckling and free vibration analysis of laminated composite and sandwich plates, Compos. Struct. 117 (2014) 316–332. [39] R. Sahoo, B. N. Singh, A new trigonometric zigzag theory for the static analysis of laminated composite and sandwich plates, Aerosp. Sci. Tech. 35 (2014) 15-28. [40] V. V. Bolotin, The dynamic stability of elastic systems, 1964, Holden-Day, San Francisco. [41] R. M. Evan-iwanowski, On the parametric response of structures, Appl. Mech. Rev. 18 (1965) 699-702. [42] R. A. Ibrahim, Structural dynamics with parameter uncertainties, Appl. Mech. Rev. 40 (1987) 309-328. [43] G. J. Simitses, Instability of dynamically loaded structures, Appl. Mech. Rev. 40 (1987) 1403-1408. [44] S. K. Sahu, P. K. Datta, Research advances in the dynamic stability behavior of plates and shells: 1987-2005, Appl. Mech. Rev. ASME 60 (2007) 65-75. [45] J. M. Hutt, A. E. Salam, Dynamic stability of plates by finite element method, ASCE J. 3 (1971) 879–899. [46] P. J. Deolasi, P. K. Datta, Parametric instability of rectangular plates subjected to localized edge compressing (compression or tension), Compos. Struct. 54 (1995) 73–82. [47] R. S. Srinivasan, P. Chellapandi, Dynamic stability of rectangular laminated composite plates, Comput. Struct. 23(7) (1987) 1053-1061.

28

[48] J. Moorthy, J. N. Reddy, R. H. Plaut, Parametric instability of laminated composite plates with transverse shear deformation, Int. J. Solids Struct. 26(7) (1990) 801-811. [49] G. Cederbaum, Dynamic stability of shear deformable laminated plates, AIAA J. 29 (1991) 2000-2005. [50] Y. W. Kwon, Finite element analysis of dynamic instability of layered composite plates using a high order bending theory, Compos. Struct. 38 (1991) 57–62. [51] A. Chattopadhay, A. G. Radu, Dynamic instability of composite laminates using a higher order theory, Compos. Struct. 77 (2000) 453-60. [52] S. K. Sahu, P. K. Datta, Dynamic instability of laminated composite rectangular plates subjected to non-uniform harmonic in-plane edge loading, J. Aerosp. Eng. Proceedings of IMeche Part G, 214 (2000) 295-312. [53] S. Wang, D. J. Dawe, Dynamic instability of composite laminated rectangular plates and prismatic plate structures, Comput. Meth. Appl. Mech. Eng. 191 (2002) 1791-1826. [54] V. Balamurugan, M. Ganapathi, T. K. Vardan, Nonlinear dynamic stability analysis of composite plates using finite element method, Comput. Struct. 60(1) (1996) 125-130. [55] M. Ganapathi, P. Boisse, D. Solaut, Nonlinear dynamic stability analysis of composite laminates under periodic in-plane compressive loads, Int. J. Numer. Meth. Eng. 46 (1999) 943956. [56] B. P. Patel, M. Ganapathi, K. R. Prasad, V. Balamurugan, Dynamic instability of layered anisotropic composite plates on elastic foundations, Eng. Structs. 21 (1999) 988-995. [57] A. Chakrabarti, A. H. Sheikh, Dynamic instability of sandwich plates using an efficient finite element model, Thin wall-Struct. 44 (2006) 57-68. [58] J. Fazilati, H. R. Ovesy, Dynamic Instability Analysis of Composite Laminated Thin-walled Structures using Two versions of FSM, Compos. Struct. 92(9) (2010) 2060–2065. [59] L. S. Ramachandra, S. K. Panda, Dynamic instability of composite plates subjected to nonuniform in-plane loads, J. Sound Vib. 331(1) (2012) 53–65. [60] M. H. Noh, S. Y. Lee, Dynamic instability of delaminated composite skew plates subjected to combined static and dynamic loads based on HSDT, Compos.: Part B 58 (2014) 113–121. [61] M. Aydogdu, A new shear deformation theory for laminated composite plates, Compos. Struct. 89 (2009) 94-101. 29

[62] N. Grover, B. N. Singh, D. K. Maiti, New non-polynomial shear-deformation theories for the structural behavior of laminated-composite and sandwich plates, AIAA J. 51(8) (2013) 1861– 1871. [63] M. Talha, B. N. Singh, Static response and free vibration analysis of FGM plates using higher order shear deformation theory, Appl. Math. Model. 34 (2010) 3991–4011. [64] R. D. Cook, D. S. Malkus, M. E. Plesha, R. J. Witt, Concepts and applications of finite element analysis, 2005 Wiley, New York. [65] Noor AK. Free vibrations of multilayered composite plates. AIAA 1973; 11: 1038– 1109. [66]Chalak HD, Chakrabarti A, Iqbal MA, Sheikh AH. Free vibration analysis of laminated soft core sandwich plates. J Vib Acoust 2013; 135: 20–31. [67]Liew KM, Huang YQ, Reddy JN. Vibration analysis of symmetrically laminated plates based on FSDT using the moving least squares differential quadrature method. J Comput Methods Appl Mech Eng 2003; 192: 2203–2222. [68]Rodrigues JD, Roque CMC, Ferreira AJM, Carrera E, Cinefra M. Radial basis functions-finite differences collocation and a unified formulation for bending, vibration and buckling analysis of laminated plates, according to Murakami's zig-zag theory. Compos Struct 2011; 93(7): 1613–1620 [69]Noor AK. Stability of multilayered composite plates. Fiber Sci Technol 1975; 8(2): 81– 89. [70] Reddy JN, Phan ND. Stability and vibration of isotropic, orthotropic and laminated plates according to a higher order shear deformation theory. J Sound Vib 1985; 98(2): 157– 170. [71] Lo SH, Zhen W, Sze KY, Wanji C. An improved in-plane displacement model for the stability analysis of laminated composites with general lamination configuration. Compos Struct 2011; 93: 1584–1594.

30

[72]Wu Z, Chen W, Ren X. An accurate higher-order theory and C0 finite element for free vibration analysis of laminated composite and sandwich plates. Compos Struct 2010; 92: 1299–1307. [73] Singh S, Chakrabarti A. Buckling analysis of laminated composite plates using an efficient C0 FE model. Lat Am J Solids Struct 2012; 9: 353–365. [74] A. K. L. Srivastava, P. K. Datta, A. H. Sheikh, Dynamic instability of stiffened plates subjected to non-uniform harmonic edge loading, J. Sound Vib. 262 (2013) 1171-89. [75] J. N. Reddy, A simple higher order shear deformation theory for laminated composite plates, J. Appl. Mech. 51(4) (1984) 745-753. [76] J. Y. Kao, C. S. Chen, W. R. Chen, Parametric vibration response of foam-filled sandwich plates under periodic loads, Mech. Compos. Mater. 48(5) (2012) 765-782. [77] T. Kant, K. Swaminathan, Analytical solutions for free vibration of laminated composite and sandwich plates based on a higher-order refined theory, Compos. Struct. 53 (2001) 73-85. [78] M. K. Rao, K. Scherbatiuk, Y. M. Desai, A. H. Shah, Natural vibrations of laminated and sandwich plates, J. Eng. Mech. 130(11) (2004) 1268–1278.

Table Captions 31

Table 1 Non-dimensional fundamental frequency for simply supported laminated composite plate (0/90/90/0) ( a / h = 5) Table 2 Non-dimensional buckling load parameter of simply supported laminated composite plate (0/90/0) with different modular ratio ( a / h = 10 ) Table 3 Convergence and validation study of non-dimensional excitation frequency of a simplysupported isotropic plate

(

Table 4 Non-dimensional excitation frequency Ωe = Ωe a2 / h

)

ρ / E of a simply-supported

laminated composite plate (0/90/90/0) (rad/s) Table 5 Comparison of non-dimensional natural frequencies and excitation frequencies of a simply-supported un-symmetric sandwich plate (0/90/C/0/90) with (α = β = 0.0)

Figure Captions Fig. 1. Schematic diagram of a laminate Fig. 2. General lamination scheme and displacement configuration Fig. 3. Eight noded isoparametric serendipity element Fig.4. Instability region of a simply-supported isotropic plate considering (α = 0.0) with

( a / h = 100) Fig. 5. Comparison of dynamic instability region of a simply-supported laminated composite plate (0/90/90/0) for (α = 0.0) with ( a / h = 25 ) Fig. 6 Effect of modular ratio on principal instability region of a simply-supported laminated composite plate (0/90/90/0) for (α = 0.0) with ( a / h = 25 ) Fig. 7 Effect of boundary conditions on principal instability region of a laminated composite plate (0/90/90/0) for (α = 0.4) with ( a / h = 25 )

32

Fig. 8 Effect of static load factor on principal instability region of a simply-supported laminated composite plate (0/90/90/0) for ( a / h = 10 ) Fig. 9 Effect of static load factor on principal instability region of a simply-supported unsymmetric sandwich plate (0/90/C/0/90) for ( a / h = 10 ) Fig. 10 Effect of thickness ratio on principal instability region of a simply-supported unsymmetric sandwich plate [0/90/C/0/90] for (α = 0.4 ) with ( a / h = 10 ) Fig. 11 Effect of span-thickness ratio on principal instability region of a simply-supported unsymmetric sandwich plate [0/90/C/0/90] for (α = 0.4 ) with ( hc / h f = 10 ) Fig. 12 Effect of boundary conditions on principal instability region of a symmetric sandwich plate (0/90/C/90/0) for (α = 0.4 ) with ( a / h = 10 )

Table 1 Non-dimensional fundamental frequency for simply supported laminated composite plate (0/90/90/0) ( a / h = 5)

E11 / E22

References 10

20

30

40

Present

8.2693 9.5143

10.2470

10.7480

Noor (1973)

8.2100 9.5600

10.2720

10.7520

Chalak et al. (2013)

8.3456 9.5703

10.2976

10.7984

Rodrigues et al. (2011)

8.4142 9.6629

10.4013

10.9054

Liew et al. (2003)

8.2924 9.5613

10.3200

10.8490

33

Table 2 Non-dimensional buckling load parameter of simply supported laminated composite plate (0/90/0) with different modular ratio ( a / h = 10 ) References

E11 / E22 10

20

30

40

ZZTSF

9.7043

14.7160 18.6755 21.8990

Noor (1975)

9.7621

15.0191 19.3040 22.8807

Reddy and Phan (1985)

9.9406

15.2980 19.6740 23.3400

Lo et al. (2011)

9.8328

14.8801 18.8545 22.0809

Wu et al. (2010)

9.7961

15.0757 19.3761 22.9643

Singh et al. (2012)

9.6259

14.6458 18.6158 21.8527

Table 3 Convergence and validation study of non-dimensional excitation frequency of a simplysupported isotropic plate

a/ h

β

α =0.0

References U e

100

α =0.2

Ω

Ω

0.2 Present(2 × 2) Present(4 × 4) Present(6 × 6) Present(8 × 8) Present(10 × 10) Chakrabarti and Sheikh (2006)

70.2262 41.758 41.4082 41.393 41.3909 41.384

0.4 Present(2 × 2) Present(4 × 4) Present(6 × 6) Present(8 × 8) Present(10 × 10) Chakrabarti and Sheikh (2006) Hutt and Salam (1971) Srivastava et al. (2003)

L e

U e

L e

α =0.35 Ω Ω eL U e

Ω

Ω

63.5985 37.7715 37.4551 37.4413 37.4394

63.5985 37.7715 37.4551 37.4413 37.4394

56.1593 33.3114 33.0323 33.02 33.0185

58.1118 34.4805 34.1917 34.1791 34.1774

49.8288 29.5274 29.28 29.2692 29.2677

37.433

37.433

33.01

34.169

29.261

73.3062 59.9986 66.9979 52.0273 61.8258 43.6148 35.6113 39.8147 30.8403 36.7073 43.2495 35.313 39.4812 30.582 36.3998 43.2335 35.3 39.4666 30.5707 36.3864 43.2313 35.2982 39.4646 30.5692 36.3846

45.1023 26.7085 26.4848 26.475 26.4737

43.224 43 43.16

35.292 39.458 30.564 36.379 26.469 35.32 35.37 -

34

1.2 Present(2 × 2) Present(4 × 4) Present(6 × 6) Present(8 × 8) Present(10 × 10) Chakrabarti and Sheikh (2006) Hutt and Salam (1971) Srivastava et al. (2003)

84.461 50.3621 49.9402 49.9218 49.9193 49.911 49.52 49.54

42.5375 25.181 24.9701 24.9609 24.9596

79.0906 47.1094 46.7148 46.6976 46.6952

30.1211 17.8056 17.6565 17.65 17.6491

74.7964 15.0772 44.5142 8.9028 44.1413 8.8282 44.125 8.825 44.1228 8.8245

24.956 46.687 17.646 44.116 8.823 25.06 24.02 -

(

Table 4 Non-dimensional excitation frequency Ωe = Ωe a2 / h

)

ρ / E of a simply-supported

laminated composite plate (0/90/90/0) (rad/s)

a/ h α

25

0.0 0.0 0.0 0.0 0.0 0.0 0.2 0.4 0.6

β

0.0 0.3 0.6 0.9 1.2 1.5 0.06 0.12 0.18

Excitation Frequency Ωe

Present 144.4817 154.9393 164.7345 173.979 182.7565 191.1313 131.6291 117.3775 101.1372

ΩUe Wang and Kao et al. Dawe (2002) (2012) 144.57 144.36 155.03 155.64 164.83 165.12 174.08 174.43 182.87 183.21 191.25 191.75 131.71 132.12 117.45 117.96 101.20 101.84

Present 144.4817 133.2055 120.882 107.1505 91.3782 72.2408 126.7821 106.1719 80.444

Ω eL Wang and Kao et al. Dawe (2002) (2012) 144.57 144.36 133.29 133.79 120.95 121.45 107.21 107.63 91.430 91.860 72.280 72.620 126.86 127.26 106.24 106.82 80.4795 81.100 35

0.8

0.24

81.7312

81.780

82.310

40.8656

40.890

41.320

Table 5 Comparison of non-dimensional natural frequencies and excitation frequencies of a simply-supported un-symmetric sandwich plate (0/90/C/0/90) with (α = β = 0.0) hc / h f

References 4 Kao et al. (2012) Present Rao et al. (2004)

Ωe ω 14.0140 28.028 1.9427 3.8854 1.9084 -

20 Kao et al. (2012) Present Rao et al. (2004)

12.9952 25.9903 2.1461 4.2923 2.1307 -

100 Kao et al. (2012) Present Rao et al. (2004)

8.5068 17.0135 2.7871 5.5743 2.7875 -

Fig. 1 Schematic diagram of a laminate

36

Fig. 2 General lamination scheme and displacement configuration

Non-dimensional excitation frequency

Fig. 3 Eight noded isoparametric serendipity element 50

Hutt and Salam(1971) ZZTSF

45

40

35

30

25 0

0.2

0.4

0.6

0.8

1

Dynamic load factor 37

Fig.4 Instability region of a simply-supported isotropic plate considering (α = 0.0) with

( a / h = 100) 2.7

ZZTSF Moorthy et al.(1990)

2.5

ȍe/Ȧ

2.3 2.1 1.9 1.7 1.5 1.3 0

0.2

0.4

0.6 ߚ

0.8

1

1.2

Fig. 5 Comparison of dynamic instability region of a simply-supported laminated composite plate (0/90/90/0) for (α = 0.0) with ( a / h = 25 ) 3

Moorthy et al.(1990)E11/E22=40 Present_E11/E22=40 Moorthy et al.(1990)E11/E22=30 Present_E11/E22=30 Moorthy et al.(1990)E11/E22=20 Present_E11/E22=20 Moorthy et al.(1990)E11/E22=10 Present_E11/E22=10

2.5

ȍe/Ȧ

2 1.5 1 0.5 0

0.2

0.4

0.6

ߚ

0.8

1

1.2

Fig. 6 Effect of modular ratio on principal instability region of a simply-supported laminated composite plate (0/90/90/0) for (α = 0.0) with ( a / h = 25 )

38

Non-dimensional excitation frequency

80 70 60 50 40 30 20 SSSS SSSC CCCC

10 0 0

0.2

0.4

0.6

0.8

1

Dynamic load factor Fig. 7 Effect of boundary conditions on principal instability region of a laminated composite plate (0/90/90/0) for (α = 0.4) with ( a / h = 25 )

Non-dimensional excitation frequency

40 35 30 25 20 ɲ=0.0 ɲ=0.2 ɲ=0.4

15 10 0

0.3

0.6

Dynamic load factor

0.9

1.2

Fig. 8 Effect of static load factor on principal instability region of a simply-supported laminated composite plate (0/90/90/0) for ( a / h = 10 )

39

2.4 2.2

ёe/ʘ

2 1.8 1.6 1.4

ɲ=0.0 ɲ=0.2 ɲ=0.4 ɲ=0.6

1.2 1 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

ɴ

Non-dimensional excitation frequency

Fig. 9 Effect of static load factor on principal instability region of a simply-supported unsymmetric sandwich plate (0/90/C/0/90) for ( a / h = 10 ) 5.8 5.3 4.8 4.3 3.8 3.3 2.8 hc/hf=4 hc/hf=20 hc/hf=30 hc/hf=50

2.3 1.8 1.3 0

0.2

0.4

0.6

0.8

1

Dynamic load factor Fig. 10 Effect of thickness ratio on principal instability region of a simply-supported unsymmetric sandwich plate [0/90/C/0/90] for (α = 0.4 ) with ( a / h = 10 )

40

Non-dimensional excitation frequency

30

a/h=10 a/h=20 a/h=40 a/h=50 a/h=100

25 20 15 10 5 0 0

0.2

0.4

0.6

0.8

1

Dynamic load factor Fig. 11 Effect of span-thickness ratio on principal instability region of a simply-supported unsymmetric sandwich plate [0/90/C/0/90] for (α = 0.4 ) with ( hc / h f = 10 )

Non-dimensional excitation frequency

23.5 22.5 21.5 20.5 19.5 18.5 17.5 16.5 15.5

SSSS SSSC SSCC CCCC

14.5 13.5 0

0.2

0.4

0.6

0.8

1

Dynamic load factor

Fig. 12 Effect of boundary conditions on principal instability region of a symmetric sandwich plate (0/90/C/90/0) for (α = 0.4 ) with ( a / h = 10 )

41