On the Free Vibration Analysis of Laminated Composite and Sandwich Plates: A Layerwise Finite Element Formulation

In this paper, a new higher-order layerwise finite element model, developed earlier by the present authors for the static analysis of laminated composite and sandwich plates, is extended to study the free vibration behavior of multilayer sandwich plates. In the present layerwise model, a first-order displacement field is assumed for the face sheets, whereas a higher-order displacement field is assumed for the core. Thanks for enforcing the continuity of the interlaminar displacement, the number of variables is independent of the number of layers. In order to reduce the computation effort, a simply four-noded C0 continuous isoparametric element is developed based on the proposed model. In order to study the free vibration, a consistent mass matrix is adopted in the present formulation. Several examples of laminated composite and sandwich plate with different material combinations, aspect ratios, boundary conditions, number of layers, geometry and ply orientations are considered for the analysis. The performance and reliability of the proposed formulation are demonstrated by comparing the author’s results with those obtained using the three-dimensional elasticity theory, analytical solutions and other advanced finite element models. From the obtained results, it can be concluded that the proposed finite element model is simple and accurate in solving the free vibration problems of laminated composite and sandwich plates.


INTRODUCTION
Due to their low weight, high stiffness and high strength properties, the composites sandwich structures are widely used in various industrial areas e.g. civil constructions, marine industry, automobile and aerospace applications. A sandwich is a three layered construction, where a low weight thick core layer (e.g., rigid polyurethane foam) of adequate transverse shear rigidity, is sandwiched between two thin laminated composite face layers of higher rigidity (Pal and Niyogi 2009). Despite the many advantages of sandwich structures, their behavior becomes very complex due to the large variation of rigidity and material properties between the core and the face sheets. Therefore, the accuracy of the results for sandwich structures largely depends on the computational model adopted (Pandey and Pradyumna 2015).
In the literature, several two-dimensional theories and approaches have been used to study the behavior of composite sandwich structures. Starting by the simple classical laminated plate theory (CLPT), based on the Kirchhoff's assumptions, which does not includes the effect of the transverse shear deformation, the first-order shear deformation theory (FSDT), where the effect of the transverse shear deformation is considered (Reissner 1975, Whitney and Pagano 1970, Mindlin 1951, Yang et al. 1966, but this theory gives a state of constant shear stresses through the plate thickness, and the higher-order shear deformation theories (HSDT) where a better representation of transverse shear effect can be obtained (Lo et al. 1977, Manjunatha and Kant 1993, Reddy 1984, Lee and Kim 2013. Regarding the approaches used to model the behavior of composite structures, we distinguish the equivalent single layers (ESL) approach where all the laminate layers are referred to the same degrees-of-freedom (DOFs). The main advantages of ESL models are their inherent simplicity and their low computational cost, due to the small number of dependent variables. However, ESL approach fails to capture precisely the local behavior of sandwich structures. This drawback in ESL was circumvented by the Zig-Zag (ZZ) and layerwise (LW) approaches in which the variables are linked to specific layers (Belarbi and Tati 2015, Belarbi et al. 2016, Ćetković and Vuksanović 2009, Chakrabarti and Sheikh 2005, Kapuria and Nath 2013, Khalili et al. 2014, Khandelwal et al. 2013, Marjanović and Vuksanović 2014, Maturi et al. 2014, Sahoo and Singh 2014, Singh et al. 2011, Thai et al. 2016. For more details, the reader may refer to (Carrera 2002, Ha 1990, Khandan et al. 2012, Sayyad and Ghugal 2015. In the last decades, the finite element method (FEM) has become established as a powerful method and as the most widely used method to analyze the complex behavior of composite sandwich structures (e.g., bending, vibration, buckling). This is due to the limitations in the analytical methods which are applicable only for certain geometry and boundary conditions (Kant and Swaminathan 2001, Mantari and Ore 2015, Maturi et al. 2014, Noor 1973, Plagianakos and Papadopoulos 2015, Srinivas and Rao 1970. Khatua and Cheung (1973) were one of the first to use the FEM in the analysis of this type of structures. They developed two triangular and rectangular elements to study the bending and free vibration of sandwich plates.
In the recent years, many researchers have investigated the dynamic response of laminated composite and sandwich plates using finite element models based on Zig-Zag theory. Chakrabarti and Sheikh (2004) developed a C 1 continuous six-noded triangular plate element with 48 DOFs for dynamic analyses of laminate-faced sandwich plate using higher-order zig-zag theory (HZZT). Afterwards, Kulkarni and Kapuria (2008) extended the application of a newly improved discrete Kirchoff quadrilateral element, based on third order zigzag theory to vibration analysis of composite and sandwich plates. Wanji (2006, 2010) carried out free vibration analyses of laminated composite and sandwich plates using an eight-noded quadrilateral element based on a global-local higher order shear deformation theories (GLHSDT). Chalak et al. (2013) developed a nine node finite element by taking out the nodal field variables in such a manner to overcome the problem of continuity requirement of the derivatives of transverse displacements for the free vibration analysis of laminated sandwich plate. An efficient nine-noded quadratic element with 99 DOFs is developed by Khandelwal et al. (2013) for accurately predicting natural frequencies of soft core sandwich plate. The formulation of this element is based on combined theory, HZZT and least square error (LSE) method. Therefore, the zig-zag plate theory presents good performance but it has a problem in its finite element implementation as it requires C 1 continuity of transverse displacement at the nodes which involves finite element implantation difficulties. Also, it requires high-order derivatives for displacement when obtaining transverse shear stresses from equilibrium equations (Pandey and Pradyumna 2015).
Recently, various authors have adopted the layerwise approach to assume separate displacement field expansions within each material layer. Lee and Fan (1996) described a new layerwise model using the FSDT for the face sheets whereas the displacement field at the core is expressed in terms of the two face sheets displacements. They used a nine-nodded isoparametric finite element for bending and free vibration of sandwich plates. On the other hand, a new three-dimensional (3D) layerwise finite element model with 240 DOFs has been developed by Nabarrete et al. (2003) for dynamic analyses of sandwich plates with laminated face sheets. They used the FSDT for the face sheets, whereas for the core a cubic and quadratic function for the in-plane and transverse, displacements, was adopted. In the same year, Desai et al. (2003) developed an eighteen-node layerwise mixed brick element with 108 DOFs for the free vibration analysis of multi-layered thick composite plates. Later, an eight nodes quadrilateral element having 136 DOFs was developed by Araújo et al. (2010) for the analysis of sandwich laminated plates with a viscoelastic core and laminated anisotropic face layers. The construction of this element is based on layerwise approach where the HSDT is used to model the core layer and the face sheets are modeled according to a FSDT. Elmalich and Rabinovitch (2012) have undertaken an analysis on the dynamics of sandwich plates, using a C0 four-node rectangular element. The formulation of this element is based on the use of a new layerwise model, where the FSDT is used for the face sheets and the HSDT is used for the core. In 2015, Pandey and Pradyumna (2015) presented a new higher-order layerwise plate formulation for static and free vibration analyses of laminated composite plates. A high-order displacement field is used for the middle layer and a firstorder displacement field for top and bottom layers. The authors used an eight-noded isoparametric element containing 104 DOFs to model the plate. The performance of these layerwise models is good but it requires high computational effort as the number of variables dramatically increases with the number of layers.
According to the presented literature review on the sandwich models, we found that many authors used finite element models having large number of nodes and/or DOFs, especially those based on the layerwise approach. Therefore, the present work aims to propose a new C 0 layerwise model competitor to the majority of aforementioned finite element models, having a reduced number of nodes and DOFs. This new model is used for the calculation of natural frequencies of laminated composite and sandwich plates. Thanks for enforcing the continuity of the interlaminar displacement, the number of variables is independent of the number of layers. The numerical results obtained by developed model are compared favorably with those obtained via analytical solution and numerical results obtained by other models. The results obtained from this investigation will be useful for a more understanding of the bending and free vibration behavior of sandwich laminates plates.

MATHEMATICAL MODEL
Sandwich plate is a structure composed of three principal layers as shown in Figure.1, two face sheets (top-bottom) of thicknesses (ht), (hb) respectively, and a central layer named core of thickness (hc) which is thicker than the previous ones. Total thickness (h) of the plate is the sum of these thicknesses. The plane (x, y) coordinate system coincides with mid-plane plate.

Displacement Field
In the present model, the core layer is modeled using the HSDT. Hence, the displacement field is written as a third-order Taylor series expansion of the in-plane displacements in the thickness coordinate, and as a constant one for the transverse displacement: where 0 u , 0 v and 0 w are respectively, in-plane and transverse displacement components at the midplane of the sandwich plate. For the two face sheets, the FSDT is adopted. The compatibility conditions as well as the displacement continuity at the interface (top face sheet-core-bottom face sheet), leads to the following improved displacement fields (Fig. 2): where t x  and t y  are the rotations of the top face-sheet cross section about the y and x axis, respectively, and the displacement of the core for (z=hc/2) is given by: The substitution of Eq. (3) in Eq.
(2) led finally to the following expressions: According to Figure 2, the displacement components of the bottom face-sheet can be written as:  are the rotations of the bottom face-sheet cross section about the y and x axis respectively, and the displacement of the core for (z=-hc/2) is given by: Substituting equation (6) into equation (5), leads to the following expression:

Strain-Displacement Relations
The strain-displacement relations derived from the displacement model of Eqs. (1), (4) and (7) are given as follows: For the core layer, For the top face sheet, For the bottom face sheet,

Constitutive Relationships
In this work, the two face sheets (top and bottom) are considered as laminated composite. Hence, the stress-strain relations for th k layer in the global coordinate system are expressed as: The core is considered as an orthotropic composite material and its loads resultants are obtained by integration of the stresses through the thickness direction of laminated plate.
where , N M denote membrane effort, bending moment, respectively, and , N M denote higher order membrane and moment resultants, respectively. V is the shear resultant; S and R are the higher order shear resultant.
It is informative to relate the loads resultants of the core defined in Eq. (12) to the total strains in Eq. (8). From Eqs. (8) and (12), we obtain: According to the FSDT, the constitutive equations for the two face sheets are given by: where the elements of reduced stiffness matrices of the face sheets are given as follows: a. Top face sheet,

FINITE ELEMENT FORMULATION
In the present study, a four-node C 0 continuous quadrilateral element, named QSFT52 (Quadrilateral Sandwich First Third with 52-DOFs), with thirteen DOFs per node has been developed. Each node contains: two rotational DOF for each face sheet, six rotational DOF for the core, while the three translations DOF are common for sandwich layers ( Figure.3). The displacements vectors at any point of coordinates (x, y) of the plate are given by: where i  is the nodal unknown vector corresponding to node i (i = 1, 2 3 4), Ni is the shape function associated with the node i. The generalized strain vector for three layers can be expressed in terms of nodal displacements vector as follows: where the matrices  relate the strains to nodal displacements.

GOVERNING DIFFERENTIAL EQUATION
In this work, Hamilton's principle is applied in order to formulate governing free vibration problem, which is given as: where t is the time, T is the kinetic energy of the system and U is the potential energy of the system. The first variation of kinetic energy of the three layers sandwich plate can be expressed as: where ui, vi and wi are the displacement in x,y and z directions, respectively, of the three-layered sandwich (i = t, c, b), i  and i V are the density of the material and volume, of each component, respectively, and ( .. ) is a second derivative with respect to time. The first variation of the potential energy of the sandwich plate is the summation of contribution from the two face sheets and from the core as: In the present analysis, the work done by external forces and the damping are neglected. Hence, Eq. (20) leads to the following dynamic equilibrium equation of a system.
where   e M and   e K denote the total element mass matrix and the total element stiffness matrix respectively, which are computed using the Gauss numerical integration. The total element stiffness matrix is the summation of contribution from the two face sheets and from the core as: where the element stiffness matrix of the core For the two face sheets, the element stiffness matrix can be written as: a. Top face sheet: coupling membrane-bending coupling bending-membrane coupling membrane-bending coupling bending-membrane The total element mass matrix, for the three-layer sandwich plate, can be written as  are the consistent mass matrices of the top face sheet, core and the bottom face sheet, respectively, containing inertia terms. Now, after evaluating the stiffness and mass matrices for all elements, the governing equations for free vibration analysis can be stated in the form of generalized eigenvalue problem.

     
where,  denote the natural frequency,   K is he global stiffness matrix,   M is the global mass matrix,    are the vectors defining the mode shapes.

NUMERICAL RESULTS AND DISCUSSIONS
In this section, several examples on the free vibration analysis of laminated composite and sandwich plates will be analyzed to demonstrate the performance and the versatility of the developed finite element model. The MATLAB programming language is used to solve the eigenvalue problem. The obtained numerical results are compared with the analytical solutions and others finite elements numerical results found in the literature. Table 1 shows the boundary conditions, for which the numerical results have been obtained, where CCCC, SSSS, CSCS and CFCF respectively indicate: fully clamped, fully simply supported, two opposite edges clamped and other two simply supported, two opposite edges clamped and other two free. Table 2 shows the material models (MM) considered for different numerical evaluation.

Boundary conditions Abbreviations Restrained edges
Simply supported SSSS
The results of the comparison show the performances and convergence of the present formulation.   In this problem, a simply supported square sandwich plate having two laminated stiff layers is investigated. The thickness of each laminate layer is 0.05h, whereas the thickness of the core is 0.8h. The mechanical properties MM3 and MM4 of table 2 are adopted, respectively, for laminated face sheets and core. The non-dimensional natural frequencies, for the first six modes, are presented in table 4 using a mesh size of 12×12. In the present analysis, different thickness ratios (a/h = 6.67, 10 and 20) are considered. The obtained results are compared with the 3D-elasticity solution given by Kulkarni and Kapuria (2008), analytical results of Wang et al. (2000) using p-Ritz method and some existing 2 c c a E h     finite element results based on HZZT Sheikh 2004, Kulkarni andKapuria 2008  Moreover, the same plate is analyzed by considering two different boundary conditions, CCCC and SCSC. The non-dimensional natural frequencies, for the first six modes, are reported in table 5 for different thickness ratios (a/h = 5, 10 and 20). The first six mode shapes obtained for SSSS, CFCF and CFFF square laminated sandwich plate with a/h =10 are shown in Figures 4, 5 and 6. It can be observed that, in comparison with the FEM solution based on HZZT (Khandelwal et al. 2013, Chalak et al. 2013, Kulkarni and Kapuria 2008, Chakrabarti and Sheikh 2004, the present element gives more accurate results than the other models.

Skew Laminated Plates
In order to evaluate the performance of the developed element for the study of free vibration response of irregular plates, a five layer symmetric cross-ply skew laminated plates (90/0/90/0/90) with simply supported edges is considered. The geometry of the skew plates is shown in Figure 7. The material properties MM5 of Table 2 is used for this analysis. The skew angle α is varied from 0°, 15°, 30°, 45° and 60°. The non-dimensional natural frequencies for the first four modes are reported in Table 6, considering the thickness ratios (a/h) as 10. A mesh size of 12×12 is considered for the analysis. The first six flexural mode shapes obtained for α = 45° are shown in Figure. 8. The comparison was made with the analytical solutions of Wang (1997) using B-spline Rayleigh-Ritz method, the solution of Ferreira et al. (2005) based on Radial Basic Function (RBF), as well as with the finite element models of Nguyen-Van (2009) and Garg et al. (2006). The results of the comparison show the effectiveness of the present element in the analysis of this type of structures.
are obtained for thickness ratio a/h = 5. The non-dimensional results of natural frequencies are reported in Table 4 using a 12×12 mesh.
It is clear, from the table 7, that the results obtained from developed element are in excellent agreement when compared with those obtained from the 3D-elasticity solution given by Noor (1973), the FEM-Q9 and Q4 solution based on LW (Marjanović and Vuksanović 2014), the FEM-Q8 solution based on GLHSDT (Zhen et al. 2010), the FEM-Q9 and Q4 solution based on HSDT (Nayak et al. 2002) and other analytical results (Owen and Li 1987, Vuksanović 2000, Matsunaga 2000. From Figure 9, it can be seen that the values of natural frequencies of laminated plate increase with increasing in E11/E22 modular ratios, whatever the number of layer.

References FE Models
No. of layers 3 10 20 30 Present element Noor (1973) Owen and Li (1987) Vuksanović (2000) Marjanović and Vuksanović (2014) Marjanović and Vuksanović ( Noor (1973) Owen and Li (1987) Vuksanović (2000) Marjanović and Vuksanović (2014) Marjanović and Vuksanović (2014) Nayak et al. (2002) Matsunaga (    To study the effect of the core thickness ratio (hc/hf) and aspect ratio (a/b) on the fundamental frequencies, a simply supported square sandwich plate with unsymmetric laminated face sheets and isotropic core is considered. The mechanical properties MM7 and MM8 of table 2 are adopted, respectively, for the laminated face sheets and the core. The thickness ratio (a/h) is taken to be 10. A comparison has been made with 3D-elasticity solution of Rao et al. (2004) to assess the suitability of the present formulation. Figure 10 shows the effect of the core thickness on the fundamental frequency of vibration. It is seen that the values of non-dimensional natural frequencies, increase with the increasing in hc/hf ratio. Further, the same sandwich plate was analyzed for different aspect ratios (a/b) keeping the same ratios a/h = 10 and hc/hf = 10. Figure 11 shows the effect of aspect ratio on the fundamental frequency. It is found that the variation of the fundamental frequency decrease with increase in aspect ratio. It is concluded that, from Figure 10 and 11, the present results are in very close agreement with the 3D-elasticity solution (Rao et al. 2004).

CONCLUSION
A new higher-order layerwise finite element model was proposed for free vibration analysis of laminated composite and sandwich plates. The developed model is based on a proper combination of higher-order and first-order, shear deformation theories. These combined theories satisfy interlaminar displacement continuity. Although the model is a layerwise one, the number of variables is independent of the number of layers. Thus, the plate theory enjoys the advantage of a single-layer plate theory, even though it is based on the concept of a layerwise plate approach. Based on this model, a fournoded C 0 continuous isoparametric element is formulated. The performance and the efficiency of the newly developed FE model are demonstrated by several numerical examples on free vibration analysis of laminated composite, symmetric/unsymmetric sandwich and skew plates, with varying material combinations, aspect ratios, number of layers, geometry and boundary conditions. The results obtained by our model were compared with those obtained by the analytical results and other finite element models found in literature. The comparison showed that the element has an excellent accuracy and a broad range of applicability. It is important to mention here, that the proposed FE formulation is simple and accurate in solving the free vibration problems of laminated composite and sandwich plates.