Free Vibrational Analysis of a Functionally Graded Five-Layer Sandwich Plate Resting on a Winkler Elastic Foundation in a Thermal Environment

Abstract

The effect of adhesive layers bonding to the core of functionally graded (FG) surface layers is investigated using the free vibration of a five-layer sandwich composite plate resting on a Winkler elastic foundation in a thermal environment. It is assumed that all layers are experiencing a steady-state temperature ΔT. The layer-wise theory is used to derive the governing equations with the help of Hamilton’s principle. The Navier solution is employed to obtain the closed-form solutions. The numerical results obtained using the present theory are compared with three-dimensional finite elements implemented by ABAQUS software. The results show that the proposed theory is not only accurate but also efficient in predicting the natural frequencies of sandwich plates resting on Winkler foundations.

1. Introduction

Functionally graded materials (FGMs) are special composites made of two or more materials with properties that vary spatially. Due to their superior characteristics, FGMs are used for specific applications in space planes, space structures, civil structures, nuclear reactors, turbine rotors, flywheels, gears, and thermal barrier systems. Due to the high demand for these advanced materials, research on understanding their mechanical behavior has become most popular in recent years [1]. It is well known that functionally graded materials are capable of resisting high-temperature environments or extremely large temperature gradients and therefore are more suitable for use in aerospace structures and nuclear plants [2].

For plated structures using FG materials and resting on elastic foundations, scientists have proposed various kinds of models [3]. The simplest model for the elastic foundation is the Winkler model. In this model, the foundation is modeled as a series of separated springs without coupling effects [4]. To better understand the behavior of such structures, several studies were performed to analyze the behavior of these materials used in sandwich plates. Benferhat et al. [4] developed an analytical solution to study the free vibration behavior of a simply supported functionally graded plate resting on an elastic foundation by taking into account the effect of transverse shear deformations. Cetkovic and Vuksanovic [5] studied bending, free vibrations and buckling of laminated composites and sandwich plates using a layer-wise displacement model. Hafizah et al. [6] used higher-order shear deformation plate theory (HSDT) to study the free vibration of antisymmetric angle-ply laminated plates with variable thickness. Zenkour and Radwan [7] presented a four-variable refined plate theory to study the free vibration of laminated composites and soft-core sandwich plates resting on Winkler–Pasternak foundations. Moita et al. [8] developed a finite element model for vibration analysis of pure functionally graded material structures, as well as for passive damped sandwich structures, with a soft viscoelastic core between the FGM layers. The FGM layers were modeled using the classical plate theory while Reddy’s third-order shear deformation theory was used to model the core.

Liu et al. [9] studied the free vibration analysis of functionally graded sandwich plates by using a refined higher-order sandwich panel theory. Vel and Batra [10] presented a three-dimensional exact solution for free and forced vibrations of simply supported functionally graded rectangular plates. Suitable displacement functions that identically satisfied the boundary conditions were used to reduce the governing equations to a set of coupled ordinary differential equations, which were solved using the power series method. Mantari et al. [11] presented a free vibration analysis of functionally graded plates (FGPs) resting on a Pasternak elastic foundation. The displacement field was based on a novel non-polynomial higher-order shear deformation theory (HSDT). Dozio [12] dealt with the formulation of advanced two-dimensional Ritz-based models for the accurate prediction of natural frequencies of thin and thick sandwich plates with cores made of functionally graded material (FGM). Singh and Harsha [13] studied the free and forced nonlinear vibration characteristics of a functionally graded material (FGM) sandwich plate resting on a Pasternak elastic foundation. The formulation was based on non-polynomial higher-order shear deformation theory with inverse hyperbolic shape function. A new modified sigmoid law was presented to compute the effective material properties of sandwich FGM plates. Sobhy [14] studied the vibration and buckling behavior of exponentially graded material (EGM) sandwich plates resting on Pasternak elastic foundations under various boundary conditions. The EGM sandwich plate was assumed to be made of a fully ceramic core sandwiched between metal/ceramic EGM coats.

Zhai et al. [15] analyzed the free vibration of two kinds of five-layered composite sandwich plates with two-layered viscoelastic cores based on the first-order shear deformation theory. Raisi et al. [16] studied stress distribution in a five-layer sandwich plate with FG face sheets subjected to a uniform transverse load using the layer-wise method and FSDT. Kulkarni and Kapuria [17] developed a new improved discrete Kirchhoff quadrilateral element based on the third-order zigzag theory for the static analysis and free vibration response of composite and sandwich plates. The free vibration of functionally graded material (FGM) rectangular sandwich plates resting on simply supported and clamped edges was studied by Li et al. [18] using the Ritz method.

The static deformations and free vibration of shear flexible isotropic and laminated composite plates were studied by Ferreira et al. [19]. They used the layer-wise method in their analysis. Reddy and Cheng [20] obtained a three-dimensional solution of a smart functionally graded plate consisting of a plate made of functionally graded material and actuators made of active material.

Tinh Quoc Bui et al. [21] studied the high frequency of plates based on the Reissner–Mindlin plate theory. High frequencies of plates are numerically explored through numerical examples for both thick and thin plates with different boundary conditions. A novel nonlocal shear deformation theory is established by Le Kha Hoa et al. [22] to investigate functionally graded nanoplates. The numerical solutions of simply supported rectangular functionally graded material nanoplates are carried out by applying the Navier procedure. Wenbin Zhou et al. [23] studied thermomechanical analysis in threads of porous metal–ceramic functionally graded composite joints by ABAQUS codes. Adelina Miteva et al. [24] reviewed some aerospace applications of functionally graded materials and their paper is motivated by the huge interest in the rapidly developing field of material science, namely, functionally graded (or functionally gradient) materials (FGMs). The extensive survey of the literature reveals that the performed analyses so far have neglected the presence of an adhesive layer in their formulation. Hence, the aim of this work is to extend the static analysis performed by Raisi et.al [16] on a five-layer sandwich composite plate and determine the free vibrational behavior of this structure resting on a Winkler elastic foundation when exposed to a thermal environment.

2. Problem Formulation

Figure 1 illustrates a sandwich plate resting on a Winkler elastic foundation. The face sheets shown in this figure have functional properties with exterior surfaces made of pure ceramic. The variations in mechanical properties of the face sheets in the z-direction (thickness direction) for the top and bottom layers ($z^{\left(k\right)} ,k=1.5$) are defined based on the origin of the global coordinate system located at the center of the core (along one of the left edges of the sandwich plate, as shown in Figure 1). Moreover, the effective material properties for each layer, namely Young’s modulus E, coefficient of thermal expansion $\alpha$, and density ρ are expressed according to Equation (1).

$$P\left(z\right)=\left(P_c-P_m\right)V_c+P_m$$

(1)

where $P_c$ and $P_m$ refer to the material properties of the ceramic and metal constituents. $V_c$ denotes the volume fraction of the ceramic phase within the face sheet. The volume fractions of the FGMs are assumed to obey a power-law function in the thickness direction, as in Equation (2).

$$V_c=\{\begin{array}{c}{\left(\frac{1}{2}-\frac{z^{\left(1\right)}}{h_1}\right)}^n k=1 -h_1/2\le z^{\left(1\right)}\le h_1/2 \\ {\left(\frac{1}{2}+\frac{z^{\left(5\right)}}{h_5}\right)}^n k=5 -h_5/2\le z^{\left(5\right)}\le h_5/2 \end{array}$$

(2)

where “n” is the power-law index. Note that when n = 0, one obtains a fully homogeneous ceramic face sheet.

The layer-wise theory (LT) adopted in this study uses the first-order shear deformation theory (FSDT) in each layer and the imposition of displacement continuity at the layers’ interfaces. According to the first-order shear deformation theory, the planar displacement components along the thickness of each layer change according to the first-order polynomials [25]. These polynomials are written in terms of the z-coordinate described locally for each layer along the plate thickness. Assuming a continuous displacement across the layers, the displacement fields in each layer are deduced and presented in Appendix A. Here, it was assumed that u and v are the x and y components of the total displacement measured at the core mid-plane (z = 0) where, $u^{\left(k\right)}$, $v^{\left(k\right)}$, and $w^{\left(k\right)}$ indicate the displacement components associated with the kth layer. The layers are numbered in consecutive order from bottom-up. Furthermore, it was assumed that the out-of-plane displacement $w^{\left(k\right)}$ is only a function of x, y and t (but not z) and there is a perfect bond between any two neighboring layers. Based on the linear elasticity, deformations for the kth layer are given by Equation (3) [25].

$$\begin{array}{l}{\epsilon }_{xx}^{\left(k\right)}={\frac{\partial u}{\partial x}}^{\left(k\right)}, {\epsilon }_{yy}^{\left(k\right)}={\frac{\partial \mathrm{v}}{\partial y}}^{\left(k\right)} , {\epsilon }_{zz}^{\left(k\right)}=0 \\ {\gamma }_{yz}^{\left(k\right)}={\frac{\partial \mathrm{v}}{\partial z}}^{\left(k\right)}+{\frac{\partial w}{\partial y}}^{\left(k\right)}={\frac{\partial \mathrm{v}}{\partial z}}^{\left(k\right)}+\frac{dw}{dy}\\ {\gamma }_{xz}^{\left(k\right)}={\frac{\partial u}{\partial z}}^{\left(k\right)}+{\frac{\partial \mathrm{w}}{\partial x}}^{\left(k\right)}={\frac{\partial u}{\partial z}}^{\left(k\right)}+\frac{dw}{dx} \\ {\gamma }_{xy}^{\left(k\right)}={\frac{\partial u}{\partial y}}^{\left(k\right)}+{\frac{\partial \mathrm{v}}{\partial x}}^{\left(k\right)}\end{array}$$

(3)

Substituting for displacement components (from Appendix A) in Equation (3), the strain vectors in each layer of the five-layer composite sandwich plate can be written as:

$$\left[\begin{array}{c}{\epsilon }_{xx}^{p\left(k\right)}\\ {\epsilon }_{yy}^{p\left(k\right)}\\ {\gamma }_{xy}^{p\left(k\right)}\end{array}\right]=\left[\begin{array}{c}{\epsilon }_{xx}^{0\left(k\right)}\\ {\epsilon }_{yy}^{0\left(k\right)}\\ {\gamma }_{xy}^{0\left(k\right)}\end{array}\right]+z^{\left(k\right)}\left[\begin{array}{c}{\epsilon }_{xx}^{1\left(k\right)}\\ {\epsilon }_{yy}^{1\left(k\right)}\\ {\gamma }_{xy}^{1\left(k\right)}\end{array}\right]-\left[\begin{array}{c}\alpha {\left(z\right)}^k\Delta \mathrm{T}\\ \alpha {\left(z\right)}^k\Delta \mathrm{T}\\ 0\end{array}\right]$$

(4a)

$$\left[\begin{array}{c}{\gamma }_{yz}^{\left(k\right)}\\ {\gamma }_{xz}^{\left(k\right)}\end{array}\right]=\left[\begin{array}{c}{\gamma }_{yz}^{0\left(k\right)}\\ {\gamma }_{xz}^{0\left(k\right)}\end{array}\right]$$

(4b)

In Equation (4a,b), the superscript “p” and “0” correspond to total strains and mid-plane strains for the elastic and isotropic FGMs, respectively. The constitutive relations can be written as in Equation (5) where (${\sigma }_{xx}^{\left(k\right)}$,${\sigma }_{yy}^{\left(k\right)},{\tau }_{yz}^{\left(k\right)},{\tau }_{xz}^{\left(k\right)},{\tau }_{xy}^{\left(k\right)}$) and (${\epsilon }_{xx}^{\left(k\right)},{\epsilon }_{yy}^{\left(k\right)},{\gamma }_{yz}^{\left(k\right)},{\gamma }_{xz}^{\left(k\right)},{\gamma }_{xy}^{\left(k\right)}$) are the stress and strain components, respectively. Using the material properties defined in Equation (1), the stiffness coefficients $(Q_{ij}^{\left(k\right)})$ can be expressed as in Equation (6).

$$\left[\begin{array}{c}{\sigma }_{xx}^{\left(k\right)}\\ {\sigma }_{yy}^{\left(k\right)}\\ {\tau }_{yz}^{\left(k\right)}\\ {\tau }_{xz}^{\left(k\right)}\\ {\tau }_{xy}^{\left(k\right)}\end{array}\right]=\left[Q_{ij}^{\left(k\right)}\right]\left(\left[\begin{array}{c}{\epsilon }_{xx}^{\left(k\right)}\\ {\epsilon }_{yy}^{\left(k\right)}\\ {\gamma }_{yz}^{\left(k\right)}\\ {\gamma }_{xz}^{\left(k\right)}\\ {\gamma }_{xy}^{\left(k\right)}\end{array}\right]-\left[\begin{array}{c}1\\ 1\\ 0\\ 0\\ 0\end{array}\right]\alpha {\left(z\right)}^{\left(k\right)}\mathsf{\Delta }\mathrm{T}\right)$$

(5)

$$\left[Q_{ij}^{\left(k\right)}\right]=\left[\begin{array}{ccccc}{Q}_{11}^{\left(k\right)}& Q_{12}^{\left(k\right)}& 0& 0& 0\\ Q_{12}^{\left(k\right)}& Q_{22}^{\left(k\right)}& 0& 0& 0\\ 0& 0& Q_{44}^{\left(k\right)}& 0& 0\\ 0& 0& 0& Q_{55}^{\left(k\right)}& 0\\ 0& 0& 0& 0& Q_{66}^{\left(k\right)}\end{array}\right]$$

(6)

where:

$$Q_{11}^{\left(k\right)}=Q_{22}^{\left(k\right)}=\frac{E^{\left(k\right)}\left(z\right)}{1-{\left(v^{\left(k\right)}\right)}^2} , Q_{12}^{\left(k\right)}=v^{\left(k\right)}{Q}_{11}^{\left(k\right)}, Q_{44}^{\left(k\right)}=Q_{55}^{\left(k\right)}=Q_{66}^{\left(k\right)}=G^{\left(k\right)}\left(z\right)$$

(7)

In Equation (5), $\mathsf{\Delta }\mathrm{T}$ is the temperature change. Note that the effective properties of the plate vary along the thickness direction according to Equation (1) and thus the elastic coefficients $Q_{ij}^{\left(k\right)}$ are a function of $z$. Hamilton’s principle is used herein to derive the equations of motion to determine the displacement field and the constitutive equations. This principle is stated in analytical form as:

$$\delta {\int }_{t_1}^{t_2}\left(U+V-T+U_f\right)dt=0$$

(8)

where U is the strain energy; T is the kinetic energy of the FG sandwich plate; $U_f$ is the strain energy of the foundation; and V is the work of external forces (for the free vibration, V = 0). It can be shown that the virtual strain energy $\delta U$ is given by:

$$\begin{gathered} \delta U=\frac{1}{2}{\int }_{\Omega }\sum _{k=1}^5{\int }_{-\frac{h_k}{2}}^{\frac{h_k}{2}}({\sigma }_{xx}^{\left(k\right)}\delta {\epsilon }_{xx}^{p\left(k\right)}+{\sigma }_{yy}^{\left(k\right)}\delta {\epsilon }_{yy}^{p\left(k\right)}+{\tau }_{xy}^{\left(k\right)}\delta {\gamma }_{xy}^{p\left(k\right)}+{\tau }_{xz}^{\left(k\right)}\delta {\gamma }_{xz}^{p\left(k\right)} \\ +{\tau }_{yz}^{\left(k\right)}\delta {\gamma }_{yz}^{p\left(k\right)}) dx dy d{z}_k ,k=1, 2, 3 ,4, 5 \end{gathered}$$

(9)

where $\Omega$ is the material volume. For the Winkler foundation model, the virtual strain energy $\delta U_f$ is given by:

$$\delta U_f={\int }_{\Omega }{f}_e\delta wd\Omega$$

(10a)

where $f_e$ is the density of the reaction force of the foundation.

For the Winkler foundation model

$$f_e=k_{w}w$$

(10b)

where $k_w$ is the parameter of the Winkler foundation.

The variation in the kinetic energy is:

$$\begin{array}{ll}{\int }_{t_1}^{t_2}\delta \mathrm{T}\text{d}t={\int }_{t_1}^{t_2}& \{{\int }_0^a{\int }_0^b{\int }_{-\frac{h_1}{2}}^{\frac{h_1}{2}}{\rho }_1\left(z\right)({\dot{u}}_1\delta {\dot{u}}_1+{\dot{v}}_1\delta {\dot{v}}_1+{\dot{w}}_1\delta {\dot{w}}_1)dxdydz\\ & +{\int }_0^a{\int }_0^b{\int }_{-\frac{h_2}{2}}^{\frac{h_2}{2}}{\rho }_2({\dot{u}}_2\delta {\dot{u}}_2+{\dot{v}}_2\delta {\dot{v}}_2+{\dot{w}}_2\delta {\dot{w}}_2)dxdydz\\ & +{\int }_0^a{\int }_0^b{\int }_{-\frac{h_3}{2}}^{\frac{h_3}{2}}{\rho }_3({\dot{u}}_3\delta {\dot{u}}_3+{\dot{v}}_3\delta {\dot{v}}_3+{\dot{w}}_3\delta {\dot{w}}_3)dxdydz\\ & +{\int }_0^a{\int }_0^b{\int }_{-\frac{h_4}{2}}^{\frac{h_4}{2}}{\rho }_4({\dot{u}}_4\delta {\dot{u}}_4+{\dot{v}}_4\delta {\dot{v}}_4+{\dot{w}}_4\delta {\dot{w}}_4)dxdydz\\ & +{\int }_0^a{\int }_0^b{\int }_{-\frac{h_5}{2}}^{\frac{h_5}{2}}{\rho }_5\left(z\right)({\dot{u}}_5\delta {\dot{u}}_5+{\dot{v}}_5\delta {\dot{v}}_5+{\dot{w}}_5\delta {\dot{w}}_5)dxdydz\}dt\end{array}$$

(11)

where $I_{\mathrm{n}}^{\left(k\right)} \left(k=1,2,3,4,5\right)$ are the inertia constants defined by:

$$I_{\mathrm{n}}^{\left(k\right)}={{\displaystyle \int }}_{-\frac{h_{\mathrm{i}}}{2}}^{\frac{h_i}{2}}{\rho }_z{z}_k^n d{z}_k$$

(12)

The stress resultants are defined as:

$$\begin{gathered} \left\{\begin{array}{c}{N}_{xx}^{\left(k\right)}\\ N_{yy}^{\left(k\right)}\\ N_{xy}^{\left(k\right)}\end{array}\right\}={\int }_{-\frac{h_k}{2}}^{\frac{h_k}{2}}\left\{\begin{array}{c}{\sigma }_{xx}^{\left(k\right)}\\ {\sigma }_{yy}^{\left(k\right)}\\ {\tau }_{xy}^{\left(k\right)}\end{array}\right\} d{z}_{\left(k\right)} \\ \left\{\begin{array}{c}{M}_{xx}^{\left(k\right)}\\ M_{yy}^{\left(k\right)}\\ M_{xy}^{\left(k\right)}\end{array}\right\}={\int }_{-\frac{h_k}{2}}^{\frac{h_k}{2}}\left\{\begin{array}{c}{\sigma }_{xx}^{\left(k\right)}\\ {\sigma }_{yy}^{\left(k\right)}\\ {\tau }_{xy}^{\left(k\right)}\end{array}\right\}{z}_{\left(k\right)} d{z}_{\left(k\right)} \\ \left\{\begin{array}{c}{\mathrm{Q}}_{xx}^{\left(k\right)}\\ {\mathrm{Q}}_{yy}^{\left(k\right)}\end{array}\right\}={\int }_{-\frac{h_k}{2}}^{\frac{h_k}{2}} K\left\{\begin{array}{c}{\tau }_{xz}^{\left(k\right)}\\ {\tau }_{yz}^{\left(k\right)}\end{array}\right\} d{z}_{\left(k\right)} \end{gathered}$$

(13)

and

$$\begin{gathered} N_{xx}^{T\left(k\right)}=N_{yy}^{T\left(k\right)},k=1,2,3,4,5 \\ {M}_{xx}^{T\left(k\right)}=M_{yy}^{T\left(k\right)},k=1,2,3,4,5 \\ {N}_{xx}^{\mathrm{T}\left(k\right)}={{\displaystyle \int }}_{-\frac{h_k}{2}}^{\frac{h_k}{2}}\frac{E\left(z\right)}{1-v}\alpha {\left(z\right)}^{\left(k\right)}\mathsf{\Delta }\mathrm{T} d{z}_{\left(k\right) }k=1,2,3,4,5 \\ {M}_{xx}^{\mathrm{T}\left(k\right)}={{\displaystyle \int }}_{-\frac{h_k}{2}}^{\frac{h_k}{2}}\frac{E\left(z\right)}{1-v}\alpha {\left(z\right)}^{\left(k\right)}\mathsf{\Delta }\mathrm{T} z_{\left(k\right) } d{z}_{\left(k\right) }k=1,2,3,4,5 \end{gathered}$$

(14)

where $\kappa$ is the shear correction factor. In this study, a shear correction factor of 0.616 was used [26]. The stress resultants defined in Equation (13) can be related to the strains defined in Equation (4) by the following equations:

$$\left[\begin{array}{c}{N}_{xx}^{p\left(k\right)}\\ N_{yy}^{p\left(k\right)}\\ N_{xy}^{p\left(k\right)}\\ M_{xx}^{p\left(k\right)}\\ M_{yy}^{p\left(k\right)}\\ M_{xy}^{p\left(k\right)}\end{array}\right]=\left[\begin{array}{cc}{\left[\mathit{A}\right]}^{\left(k\right)}& {\left[\mathit{B}\right]}^{\left(k\right)}\\ {\left[\mathit{B}\right]}^{\left(k\right)}& {\left[\mathit{D}\right]}^{\left(k\right)}\end{array}\right]\left[\begin{array}{c}{\epsilon }_{xx}^{0\left(k\right)}\\ {\epsilon }_{yy}^{0\left(k\right)}\\ {\gamma }_{xy}^{0\left(k\right)}\\ {\epsilon }_{xx}^{1\left(k\right)}\\ {\epsilon }_{yy}^{1\left(k\right)}\\ {\gamma }_{xy}^{1\left(k\right)}\end{array}\right]-\left\{\begin{array}{c}{N}_{xx}^{T\left(k\right)}\\ N_{yy}^{T\left(k\right)}\\ 0\\ M_{xx}^{T\left(k\right)}\\ M_{yy}^{T\left(k\right)}\\ 0\end{array}\right\}$$

(15)

Here, the elements of matrices ${\mathit{A}}_{ij}^{\left(k\right)}$, ${\mathit{B}}_{ij}^{\left(k\right)}$, and ${\mathit{D}}_{ij}^{\left(k\right)}$ (k = 1,…,5), can be written as:

$$({\mathit{A}}_{ij}^{\left(k\right)},{\mathit{B}}_{ij}^{\left(k\right)},{\mathit{D}}_{ij}^{\left(k\right)})={\int }_{-\frac{h_k}{2}}^{\frac{h_k}{2}}{Q}_{ij}^{\left(k\right)}\left(1,z_{\left(k\right)},z_{\left(k\right)}{}^2\right) d{z}_{\left(k\right)} i,j=1,2,6$$

(16)

Using the generalized displacement–strain relations (Equations (3) and (4)), stress–strain relations (Equations (5) and (13)), and integrating by parts as well as the fundamental lemma of variational calculus and collecting the coefficients of $\delta u,\delta \mathrm{v},\delta w,\delta {\varnothing }_x^{\left(k\right)} \mathrm{and} \delta {\varnothing }_y^{\left(k\right)}$ in Equation (8), the equations of motion are obtained as:

$$\begin{gathered} \delta u:{\displaystyle \sum }_{k=1}^5\left(\frac{\partial N_x^{\left(k\right)}}{\partial x}+\frac{\partial N_{xy}^{\left(k\right)}}{\partial y}-\frac{\partial N_x^{T\left(k\right)}}{\partial x}\right)=\left(J_1\ddot{u}+J_2{\ddot{\varnothing }}_x^{\left(1\right)}+J_3{\ddot{\varnothing }}_x^{\left(2\right)}+J_4{\ddot{\varnothing }}_x^{\left(3\right)}+J_5{\ddot{\varnothing }}_x^{\left(4\right)}+J_6{\ddot{\varnothing }}_x^{\left(5\right)}\right) \\ \delta \mathrm{v}:{\displaystyle \sum }_{k=1}^5\left(\frac{\partial N_y^{\left(k\right)}}{\partial y}+\frac{\partial N_{xy}^{\left(k\right)}}{\partial x}-\frac{\partial N_y^{T\left(k\right)}}{\partial y}\right)=\left(J_1\ddot{v}+J_2{\ddot{\varnothing }}_y^{\left(1\right)}+J_3{\ddot{\varnothing }}_y^{\left(2\right)}+J_4{\ddot{\varnothing }}_y^{\left(3\right)}+J_5{\ddot{\varnothing }}_y^{\left(4\right)}+J_6{\ddot{\varnothing }}_y^{\left(5\right)}\right) \\ \delta w:{\displaystyle \sum }_{k=1}^5\left(-\frac{\partial Q_x^{\left(k\right)}}{\partial x}-\frac{\partial Q_y^{\left(k\right)}}{\partial y}\right)+k_{w}w=J_1\ddot{w} \\ \delta {\varnothing }_x^{\left(1\right)}:\frac{h_1}{2}\frac{\partial N_x^{\left(1\right)}}{\partial x}+\frac{h_1}{2}\frac{\partial N_{xy}^{\left(1\right)}}{\partial y}-\frac{\partial M_x^{\left(1\right)}}{\partial x}-\frac{\partial M_{xy}^{\left(1\right)}}{\partial y}+Q_x^{\left(1\right)}-\left(\frac{h_1}{2}\frac{\partial N_x^{T\left(1\right)}}{\partial x}-\frac{h_1}{2}\frac{\partial M_x^{T\left(1\right)}}{\partial x}\right) \\ =\left(J_2\ddot{u}+J_7{\ddot{\varnothing }}_x^{\left(1\right)}+J_8{\ddot{\varnothing }}_x^{\left(2\right)}+J_9{\ddot{\varnothing }}_x^{\left(3\right)}\right) \\ \delta {\varnothing }_y^{\left(1\right)}:\frac{h_1}{2}\frac{\partial N_y^{\left(1\right)}}{\partial y}+\frac{h_1}{2}\frac{\partial N_{xy}^{\left(1\right)}}{\partial x}-\frac{\partial M_y^{\left(1\right)}}{\partial y}-\frac{\partial M_{xy}^{\left(1\right)}}{\partial x}+Q_y^{\left(1\right)}-\left(\frac{h_1}{2}\frac{\partial N_y^{T\left(1\right)}}{\partial y}-\frac{h_1}{2}\frac{\partial M_y^{T\left(1\right)}}{\partial y}\right) \\ =\left(J_2\ddot{v}+J_7{\ddot{\varnothing }}_y^{\left(1\right)}+J_8{\ddot{\varnothing }}_y^{\left(2\right)}+J_9{\ddot{\varnothing }}_y^{\left(3\right)}\right) \\ \delta {\varnothing }_x^{\left(2\right)}:h_2\frac{\partial N_x^{\left(1\right)}}{\partial x}+h_2\frac{\partial N_{xy}^{\left(1\right)}}{\partial y}+\frac{h_2}{2}\frac{\partial N_x^{\left(2\right)}}{\partial x}+\frac{h_2}{2}\frac{\partial N_{xy}^{\left(2\right)}}{\partial y}-\frac{\partial M_x^{\left(2\right)}}{\partial x}-\frac{\partial M_{xy}^{\left(2\right)}}{\partial y}+Q_x^{\left(2\right)} \\ -\left(h_2\frac{\partial N_x^{T\left(1\right)}}{\partial x}+\frac{h_2}{2}\frac{\partial N_x^{T\left(2\right)}}{\partial x}-\frac{\partial M_x^{T\left(2\right)}}{\partial x}\right) \\ =\left(J_3\ddot{u}+J_8{\ddot{\varnothing }}_x^{\left(1\right)}+J_{10}{\ddot{\varnothing }}_x^{\left(2\right)}+J_{11}{\ddot{\varnothing }}_x^{\left(3\right)}\right) \\ \delta {\varnothing }_y^{\left(2\right)}:\frac{h_3}{2}\frac{\partial N_x^{\left(1\right)}}{\partial x}+\frac{h_3}{2}\frac{\partial N_{xy}^{\left(1\right)}}{\partial y}+\frac{h_3}{2}\frac{\partial N_x^{\left(2\right)}}{\partial x}+\frac{h_3}{2}\frac{\partial N_{xy}^{\left(2\right)}}{\partial y}-\frac{\partial M_x^{\left(3\right)}}{\partial x}-\frac{\partial M_{xy}^{\left(3\right)}}{\partial y}+Q_x^{\left(3\right)}-\frac{h_3}{2}\frac{\partial N_x^{\left(4\right)}}{\partial x} \\ -\frac{h_3}{2}\partial N_{xy}^{\left(4\right)}\partial y-\frac{h_3}{2}\frac{\partial N_x^{\left(5\right)}}{\partial x}-\frac{h_3}{2}\frac{\partial N_{xy}^{\left(5\right)}}{\partial y} \\ -\left(\frac{h_3}{2}\frac{\partial N_x^{T\left(1\right)}}{\partial x}+\frac{h_3}{2}\frac{\partial N_x^{T\left(2\right)}}{\partial x}-\frac{h_3}{2}\frac{\partial N_x^{T\left(4\right)}}{\partial x}-\frac{h_3}{2}\frac{\partial N_x^{T\left(5\right)}}{\partial x}-\frac{\partial M_x^{T\left(3\right)}}{\partial x}\right) \\ =\left(J_3\ddot{v}+J_8{\ddot{\varnothing }}_y^{\left(1\right)}+J_{10}{\ddot{\varnothing }}_y^{\left(2\right)}+J_{11}{\ddot{\varnothing }}_y^{\left(3\right)}\right) \\ \delta {\varnothing }_y^{\left(2\right)}:h_2\frac{\partial {\mathrm{N}}_y^{\left(1\right)}}{\partial y}+h_2\frac{\partial {\mathrm{N}}_{xy}^{\left(1\right)}}{\partial x}+\frac{h_2}{2}\frac{\partial {\mathrm{N}}_y^{\left(2\right)}}{\partial y}+\frac{h_2}{2}\frac{\partial {\mathrm{N}}_{xy}^{\left(2\right)}}{\partial x}-\frac{\partial {\mathrm{M}}_y^{\left(2\right)}}{\partial y}-\frac{\partial {\mathrm{M}}_{xy}^{\left(2\right)}}{\partial x}+{\mathrm{Q}}_y^{\left(2\right)}=f_y^{T\left(2\right)} \\ \delta {\varnothing }_x^{\left(3\right)}\frac{h_3}{2}\frac{\partial N_x^{\left(1\right)}}{\partial x}+\frac{h_3}{2}\frac{\partial N_{xy}^{\left(1\right)}}{\partial y}+\frac{h_3}{2}\frac{\partial N_x^{\left(2\right)}}{\partial x}+\frac{h_3}{2}\frac{\partial N_{xy}^{\left(2\right)}}{\partial y}-\frac{\partial M_x^{\left(3\right)}}{\partial x}-\frac{\partial M_{xy}^{\left(3\right)}}{\partial y}+Q_x^{\left(3\right)}-\frac{h_3}{2}\frac{\partial N_x^{\left(4\right)}}{\partial x} \\ -\frac{h_3}{2}\frac{\partial N_{xy}^{\left(4\right)}}{\partial y}-\frac{h_3}{2}\frac{\partial N_x^{\left(5\right)}}{\partial x}-\frac{h_3}{2}\frac{\partial N_{xy}^{\left(5\right)}}{\partial y} \\ -\left(\frac{h_3}{2}\frac{\partial N_x^{T\left(1\right)}}{\partial x}+\frac{h_3}{2}\frac{\partial N_x^{T\left(2\right)}}{\partial x}-\frac{h_3}{2}\frac{\partial N_x^{T\left(4\right)}}{\partial x}-\frac{h_3}{2}\frac{\partial N_x^{T\left(5\right)}}{\partial x}-\frac{\partial M_x^{T\left(3\right)}}{\partial x}\right) \\ =\left(J_4\ddot{u}+J_9{\ddot{\varnothing }}_x^{\left(1\right)}+J_{11}{\ddot{\varnothing }}_x^{\left(2\right)}+J_{12}{\ddot{\varnothing }}_x^{\left(3\right)}+J_{13}{\ddot{\varnothing }}_x^{\left(4\right)}+J_{14}{\ddot{\varnothing }}_x^{\left(5\right)}\right) \\ \delta {\varnothing }_y^{\left(3\right)}:\frac{h_3}{2}\frac{\partial N_y^{\left(1\right)}}{\partial y}+\frac{h_3}{2}\frac{\partial N_{xy}^{\left(1\right)}}{\partial x}+\frac{h_3}{2}\frac{\partial N_y^{\left(2\right)}}{\partial y}+\frac{h_3}{2}\frac{\partial N_{xy}^{\left(2\right)}}{\partial x}-\frac{\partial M_y^{\left(3\right)}}{\partial y}-\frac{\partial M_{xy}^{\left(3\right)}}{\partial x}+Q_y^{\left(3\right)}-\frac{h_3}{2}\frac{\partial N_y^{\left(4\right)}}{\partial y} \\ -\frac{h_3}{2}\frac{\partial N_{xy}^{\left(4\right)}}{\partial x}-\frac{h_3}{2}\frac{\partial N_y^{\left(5\right)}}{\partial y}-\frac{h_3}{2}\frac{\partial N_{xy}^{\left(5\right)}}{\partial x} \\ -\left(\frac{h_3}{2}\frac{\partial N_y^{T\left(1\right)}}{\partial y}+\frac{h_3}{2}\frac{\partial N_y^{T\left(2\right)}}{\partial y}-\frac{h_3}{2}\frac{\partial N_y^{T\left(4\right)}}{\partial y}-\frac{h_3}{2}\frac{\partial N_y^{T\left(5\right)}}{\partial y}-\frac{\partial M_y^{T\left(3\right)}}{\partial x}\right) \\ =\left(J_4\ddot{v}+J_9{\ddot{\varnothing }}_y^{\left(1\right)}+J_{11}{\ddot{\varnothing }}_y^{\left(2\right)}+J_{12}{\ddot{\varnothing }}_y^{\left(3\right)}+J_{13}{\ddot{\varnothing }}_y^{\left(4\right)}+J_{14}{\ddot{\varnothing }}_y^{\left(5\right)}\right) \\ \delta {\varnothing }_x^{\left(4\right)}:-\frac{h_2}{2}\frac{\partial N_x^{\left(4\right)}}{\partial x}-\frac{h_2}{2}\frac{\partial N_{xy}^{\left(5\right)}}{\partial y}-h_2\frac{\partial N_x^{\left(5\right)}}{\partial x}-h_2\frac{\partial N_{xy}^{\left(2\right)}}{\partial y}-\frac{\partial M_x^{\left(4\right)}}{\partial x}-\frac{\partial M_{xy}^{\left(4\right)}}{\partial y}+Q_x^{\left(4\right)} \\ -\left(-\frac{h_2}{2}\frac{\partial N_x^{T\left(4\right)}}{\partial x}-h_2\frac{\partial N_x^{T\left(5\right)}}{\partial x}-\frac{\partial M_x^{T\left(4\right)}}{\partial x}\right) \\ =\left(J_5\ddot{u}+J_{13}{\ddot{\varnothing }}_x^{\left(3\right)}+J_{15}{\ddot{\varnothing }}_x^{\left(4\right)}+J_{16}{\ddot{\varnothing }}_x^{\left(5\right)}\right) \\ \delta {\varnothing }_y^{\left(4\right)}:-\frac{h_2}{2}\frac{\partial N_y^{\left(4\right)}}{\partial y}-\frac{h_2}{2}\frac{\partial N_{xy}^{\left(4\right)}}{\partial x}-h_2\frac{\partial N_y^{\left(5\right)}}{\partial y}-h_2\frac{\partial N_{xy}^{\left(5\right)}}{\partial x}-\frac{\partial M_y^{\left(4\right)}}{\partial y}-\frac{\partial M_{xy}^{\left(4\right)}}{\partial x}+Q_y^{\left(4\right)} \\ -\left(-\frac{h_2}{2}\frac{\partial N_y^{T\left(4\right)}}{\partial y}-h_2\frac{\partial N_y^{T\left(5\right)}}{\partial y}-\frac{\partial M_y^{T\left(4\right)}}{\partial y}\right) \\ =\left(J_5\ddot{v}+J_{13}{\ddot{\varnothing }}_y^{\left(3\right)}+J_{15}{\ddot{\varnothing }}_y^{\left(4\right)}+J_{16}{\ddot{\varnothing }}_y^{\left(5\right)}\right) \\ \delta {\varnothing }_x^{\left(5\right)}:-\frac{h_1}{2}\frac{\partial N_x^{\left(5\right)}}{\partial x}-\frac{h_1}{2}\frac{\partial N_{xy}^{\left(5\right)}}{\partial y}-\frac{\partial M_x^{\left(5\right)}}{\partial x}-\frac{\partial M_{xy}^{\left(5\right)}}{\partial y}+Q_x^{\left(5\right)}-\left(\frac{h_1}{2}\frac{\partial N_x^{T\left(5\right)}}{\partial x}-\frac{h_1}{2}\frac{\partial M_x^{T\left(5\right)}}{\partial x}\right) \\ =\left(J_6\ddot{u}+J_{14}{\ddot{\varnothing }}_x^{\left(3\right)}+J_{16}{\ddot{\varnothing }}_x^{\left(4\right)}+J_{17}{\ddot{\varnothing }}_x^{\left(5\right)}\right) \\ \delta {\varnothing }_y^{\left(5\right)}:-\frac{h_1}{2}\frac{\partial N_y^{\left(5\right)}}{\partial y}-\frac{h_1}{2}\frac{\partial N_{xy}^{\left(5\right)}}{\partial x}-\frac{\partial M_y^{\left(5\right)}}{\partial y}-\frac{\partial M_{xy}^{\left(5\right)}}{\partial x}+Q_y^{\left(5\right)} \\ -\left(\frac{h_1}{2}\frac{\partial N_y^{T\left(5\right)}}{\partial y}-\frac{h_1}{2}\frac{\partial M_y^{T\left(5\right)}}{\partial y}\right) \\ =\left(J_6\ddot{v}+J_{14}{\ddot{\varnothing }}_y^{\left(3\right)}+J_{16}{\ddot{\varnothing }}_y^{\left(4\right)}+J_{17}{\ddot{\varnothing }}_y^{\left(5\right)}\right) \end{gathered}$$

(17)

where

$$\begin{gathered} f_{xx1}^{T\left(k\right)}={\displaystyle \sum }_{k=1}^5\left(\frac{\partial N_x^{T\left(k\right)}}{\partial x}\right) \\ {f}_{yy1}^{T\left(k\right)}={\displaystyle \sum }_{k=1}^5\left(\frac{\partial N_y^{T\left(k\right)}}{\partial x}\right) \\ {f}_x^{T\left(1\right)}=\left(\frac{h_1}{2}\frac{\partial N_x^{T\left(1\right)}}{\partial x}-\frac{h_1}{2}\frac{\partial M_x^{T\left(1\right)}}{\partial x}\right) \\ {f}_y^{T\left(1\right)}=\left(\frac{h_1}{2}\frac{\partial N_y^{T\left(1\right)}}{\partial y}-\frac{h_1}{2}\frac{\partial M_y^{T\left(1\right)}}{\partial y}\right) \\ {f}_x^{T\left(2\right)}=\left(h_2\frac{\partial N_x^{T\left(1\right)}}{\partial x}+\frac{h_2}{2}\frac{\partial N_x^{T\left(2\right)}}{\partial x}-\frac{\partial M_x^{T\left(2\right)}}{\partial x}\right) \\ {f}_y^{T\left(2\right)}=\left(h_2\frac{\partial N_y^{T\left(1\right)}}{\partial y}+\frac{h_2}{2}\frac{\partial N_y^{T\left(2\right)}}{\partial y}-\frac{\partial M_y^{T\left(2\right)}}{\partial y}\right) \\ {f}_x^{T\left(3\right)}=\left(\frac{h_3}{2}\frac{\partial N_x^{T\left(1\right)}}{\partial x}+\frac{h_3}{2}\frac{\partial N_x^{T\left(2\right)}}{\partial x}-\frac{h_3}{2}\frac{\partial N_x^{T\left(4\right)}}{\partial x}-\frac{h_3}{2}\frac{\partial N_x^{T\left(5\right)}}{\partial x}-\frac{\partial M_x^{T\left(3\right)}}{\partial x}\right) \\ {f}_y^{T\left(3\right)}=\left(\frac{h_3}{2}\frac{\partial N_y^{T\left(1\right)}}{\partial y}+\frac{h_3}{2}\frac{\partial N_y^{T\left(2\right)}}{\partial y}-\frac{h_3}{2}\frac{\partial N_y^{T\left(4\right)}}{\partial y}-\frac{h_3}{2}\frac{\partial N_y^{T\left(5\right)}}{\partial y}-\frac{\partial M_y^{T\left(3\right)}}{\partial x}\right) \\ {f}_x^{T\left(4\right)}=\left(-\frac{h_2}{2}\frac{\partial N_x^{T\left(4\right)}}{\partial x}-h_2\frac{\partial N_x^{T\left(5\right)}}{\partial x}-\frac{\partial M_x^{T\left(4\right)}}{\partial x}\right) \\ {f}_y^{T\left(4\right)}=\left(-\frac{h_2}{2}\frac{\partial N_y^{T\left(4\right)}}{\partial y}-h_2\frac{\partial N_y^{T\left(5\right)}}{\partial y}-\frac{\partial M_y^{T\left(4\right)}}{\partial y}\right) \\ {f}_x^{T\left(5\right)}=\left(\frac{h_1}{2}\frac{\partial N_x^{T\left(5\right)}}{\partial x}-\frac{h_1}{2}\frac{\partial M_x^{T\left(5\right)}}{\partial x}\right) \\ {f}_y^{T\left(5\right)}=\left(\frac{h_1}{2}\frac{\partial N_y^{T\left(5\right)}}{\partial y}-\frac{h_1}{2}\frac{\partial M_y^{T\left(5\right)}}{\partial y}\right) \end{gathered}$$

(18)

Expressions for $J_1$ to $J_{17}$ are given in Appendix B. For simply supported edges, the boundary conditions are given as follows:

$$\begin{gathered} u^{\left(k\right)}\left(x,0,t\right)=0, u^{\left(k\right)}\left(x,b,t\right)=0, {\mathrm{v}}^{\left(\mathrm{k}\right)}\left(0,y,t\right)=0, {\mathrm{v}}^{\left(k\right)}\left(a,y,t\right)=0 \\ {w}^{\left(k\right)}\left(x,0,t\right)=0, w^{\left(k\right)}\left(x,b,t\right)=0, w^{\left(k\right)}\left(0,y,t\right)=0, w^{\left(k\right)}\left(a,y,t\right)=0 \\ {\varnothing }_x^{\left(k\right)}\left(x,0,t\right)=0, {\varnothing }_x^{\left(k\right)}\left(x,b,t\right)=0, {\varnothing }_y^{\left(k\right)}\left(0,y,t\right)=0, {\varnothing }_y^{\left(k\right)}\left(a,y,t\right)=0 \end{gathered}$$

(19)

To obtain the analytical solution of the partial differential Equations (17), the Navier method, based on the double Fourier series is used under the specified boundary conditions given in Equation (19). The solution to displacement components satisfying the above boundary conditions can be expressed in the following Fourier series as:

$$\begin{gathered} u={\displaystyle \sum }_{m=1}^{\infty }{\displaystyle \sum }_{n=1}^{\infty }{u}_{mn}\mathit{cos}\left(\alpha x\right)\mathit{sin}\left(\beta y\right)e^{i\omega t} \\ v={\displaystyle \sum }_{m=1}^{\infty }{\displaystyle \sum }_{n=1}^{\infty }{v}_{mn}\mathit{sin}\left(\alpha x\right)\mathit{cos}\left(\beta y\right)e^{i\omega t} \\ w={\displaystyle \sum }_{m=1}^{\infty }{\displaystyle \sum }_{n=1}^{\infty }{w}_{mn}\mathit{sin}\left(\alpha x\right)\mathit{sin}\left(\beta y\right)e^{i\omega t} \\ {\varnothing }_x^{\left(k\right)}={\displaystyle \sum }_{m=1}^{\infty }{\displaystyle \sum }_{n=1}^{\infty }{\varnothing }_{x mn}^{\left(k\right)}\mathit{cos}\left(\alpha x\right)\mathit{sin}\left(\beta y\right)e^{i\omega t} \\ {\varnothing }_y^{\left(k\right)}={\displaystyle \sum }_{m=1}^{\infty }{\displaystyle \sum }_{n=1}^{\infty }{\varnothing }_{y mn}^{\left(k\right)}\mathit{sin}\left(\alpha x\right)\mathit{cos}\left(\beta y\right)e^{i\omega t} \end{gathered}$$

(20)

where $\alpha =\frac{m\pi }{a}$ (m = 1, 2, …) and $\beta =\frac{n\pi }{b}$ (n = 1, 2, …), where m and n are the wave numbers along the x- and y-directions, respectively, and $\omega$ is the natural frequency of vibration.

To obtain a solution, the thermal loads and moments must be also expanded in double sine series. For this purpose, expressing $\Delta T\left(x,y,z,t\right)$ as a double sine series we have [25]:

$$\begin{gathered} \Delta \mathrm{T}\left(x,y,z,t\right)={\displaystyle \sum }_{m=1}^{\infty }{\displaystyle \sum }_{n=1}^{\infty }{T}_{mn}\mathit{sin}\left(\alpha x\right)\mathit{sin}\left(\beta y\right)e^{i\omega t} \\ where; T_{mn}=\{\begin{array}{c}\frac{16T}{mn{\pi }^2} \mathrm{odd} m,n \\ 0 \mathrm{even} m or n\end{array} \end{gathered}$$

(21)

The load components can be written as:

$$\begin{gathered} N_{xx}^{\left(k\right)T}={\displaystyle \sum }_{m=1}^{\infty }{\displaystyle \sum }_{n=1}^{\infty }{N}_{xxmn}^{\left(k\right)}\mathit{sin}\left(\alpha x\right)\mathit{sin}\left(\beta y\right)e^{i\omega t} \\ {N}_{xxmn}^{\left(k\right)}={{\displaystyle \int }}_{-\frac{h_k}{2}}^{\frac{h_k}{2}}\frac{E{\left(z\right)}^{\left(k\right)}}{1-v{\left(z\right)}^{\left(k\right)}}\alpha {\left(z\right)}^{\left(k\right)} T_{mn}{e}^{i\omega t}d{z}_{\left(k\right) } \\ {M}_{xx}^{\left(k\right)T}={\displaystyle \sum }_{m=1}^{\infty }{\displaystyle \sum }_{n=1}^{\infty }{M}_{xxmn}^{\left(k\right)}\mathit{sin}\left(\alpha x\right)\mathit{sin}\left(\beta y\right)e^{i\omega t} \\ {M}_{xxmn}^{\left(k\right)}={{\displaystyle \int }}_{-\frac{h_k}{2}}^{\frac{h_k}{2}}\frac{E{\left(z\right)}^{\left(k\right)}}{1-v{\left(z\right)}^{\left(k\right)}}\alpha {\left(z\right)}^{\left(k\right) }{T}_{mn} e^{i\omega t}{z}_{\left(k\right) }d{z}_{\left(k\right) } \end{gathered}$$

(22)

Substituting displacement functions (20) into Equation (17), the following equations are obtained in a matrix form as:

$$\left(\left[K\right]-\left[M\right]{\omega }^2\right)\{\Delta \}=\left\{0\right\}$$

(23)

where [K] and [M] are the stiffness and mass matrices, respectively. $\{\Delta \}$ is the column vector of coefficients ($u,\mathrm{v},w,{\varnothing }_x^{\left(k\right)},{\varnothing }_y^{\left(k\right)}$). The natural frequencies are obtained from the nontrivial solution of Equation (23).

3. Numerical Results and Discussions

Due to the extensive coupling between the terms in the differential Equation (17), a numerical scheme was adopted to seek a solution using a MATLAB software program. Furthermore, to investigate the validity of the results obtained, a finite element model (see Figure 2a) was prepared and solved using ABAQUS finite element software.

To input the variable properties of the FGM cover sheets, the thickness of each sheet was divided into 30 divisions, as shown in Figure 2b. Similar models were prepared with three layers in the absence of the adhesive layer to compare the results based on the solution method adopted in this work with the published results in the literature and the versatility of the layer-wise (LW) model for the solution of such problems.

As the first step, to investigate the validity of the numerical results on the natural frequencies, the natural frequencies of a three-layer composite plate were determined and compared with those in [27] obtained based on the exact solution. The three-layer model was composed of a core encapsulated by two cover sheets with a [0/90/0] stacking sequence.

3.1. Sandwich Plate with Three-Layer Composite Plate (0/90/0) Cover Sheets

To validate the solution procedure, a three-layer square shape sandwich plate with a core and two laminated [0/90/0] composite cover sheets was selected to rest on simply supported (SS) edge supports. Plate dimensions were taken as side $a$ and thickness $h$ where the side-to-thickness ratio a/h was assumed to be equal to 10. The thickness of each ply was considered to be h/3 with the material properties given in [19] and shown below.

$$\begin{gathered} E_1=173 \mathrm{MPa} ,E_2=33.1 \mathrm{MPa} ,E_3=5.17 \mathrm{MPa} G_{13}=8.27 \mathrm{MPa} \\ {G}_{23}=3.24 \mathrm{MPa}, G_{12}=9.38 \mathrm{MPa}, v_{12}=0.036,v_{23}=0.171 , v_{13}=0.25 \end{gathered}$$

(24)

Additionally, the dimensionless frequency parameter was defined as:

$$\overline{\omega }=\frac{\omega a^2}{h}\sqrt{\frac{\rho }{E_2}}$$

(25)

where $\omega$ is the circular frequency.

Table 1. Natural frequencies of an SS square plate (0/90/0) ($\frac{a}{h}=10),\overline{\omega }=\omega h\sqrt{\frac{\rho }{E_2}}$.

Reference	$\overline{\mathit{\omega }}$
	Mode 1	Mode 2	Mode 3	Mode 4
Ferreira [19]	0.0659	0.1322	0.1762	0.2150
Exact (Srinivas et al. [27])	0.06715	0.12811	0.17217	0.20798
HSDT (Nosier et al. [28])	0.06716	0.12816	0.17225	0.20808
Layer-wise (Wang and Zhang [29])	0.06716	0.12819	0.17230	0.20811
Present layer-wise formulation (Navier solution)	0.0662	0.1268	0.1661	0.2051

shows a comparison of results obtained on $\overline{\omega }$ based on the present formulation (solution) and those given by Ferreira [19], the exact solution by Srinivas et al. [27], the higher-order formulation theory (HSDT) by Nosier et al. [28] and the layer-wise B-spline finite strip method by Wang and Zhang [29]. Close agreements between the results indicate the validity of the present formulation.

3.2. Sandwich Plate with Two-Layer Cover Sheets (0/90/core/0/90)

As a second check, we now introduce a sandwich plate with an internal core encapsulated by two cover sheets with the layer arrangement of [0/90] to construct the whole structure. The resulting (0/90/core/0/90) composite plate with dimensions $a\times b\times h$ was analyzed using the layer-wise method followed by the application of a MATLAB software program to obtain a solution. The material properties of the layers and the internal core are given below.

Face sheets [30]:

$$\begin{gathered} E_1=131 \mathrm{GPa} ,E_2=E_3=10.34 \mathrm{GPa} G_{12}=G_{23}=6.895 \mathrm{GPa} \\ {G}_{13}=6.205 \mathrm{GPa} v_{12}=v_{13}=0.22,v_{23}=0.49 ,\rho =1627 \mathrm{kg}/{\mathrm{m}}^3 \end{gathered}$$

(26)

Isotropic core [30]:

$$\begin{gathered} E_1=E_2=E_3=6.89 \mathrm{MPa} G_{12}=G_{23}=G_{13}=3.45 \mathrm{MPa} \\ {v}_{12}=v_{13}=v_{23}=0,{\rho }_{\mathrm{c}}=97\frac{\mathrm{kg}}{{\mathrm{m}}^3} \end{gathered}$$

(27)

Results from the present study are compared with the published data in [30] for various values of $\frac{a}{h}$, $\frac{a}{b}$ and $\frac{h_3}{h_1}$ ratios in

Table 2. Comparison of non-dimensional fundamental frequencies $\overline{\omega }$ of (0°/90°/core/0°/90°) sandwich plate with $\frac{h_3}{h_1}=10 \mathrm{and} \frac{a}{b}=1$.

$\frac{\mathit{a}}{\mathit{h}}$	Present Formulation	Rao	Rao
	(LW)	(Refine)	(LW)
2	0.9024	0.7141	0.7132
4	1.0972	0.9363	0.936
10	1.9376	1.848	1.848
20	3.5283	3.4791	3.4793
30	5.0706	5.0371	5.0375
30	6.487	6.4634	6.4637
50	7.7509	7.7355	7.7358
60	8.8571	8.8118	8.8495
70	9.8124	9.8112	9.8121
80	10.6304	10.6368	10.6371
90	11.3276	11.3408	11.3409
100	11.92	11.94	11.9401

,

Table 3. Comparison of non-dimensional fundamental frequencies $\overline{\omega }$ of (0°/90°/core/0°/90°) sandwich plate with $\frac{a}{h}=10 and \frac{h_3}{h_1}=10$.

$\frac{\mathit{a}}{\mathit{b}}$	Present Formulation	Rao	Rao
	(LW)	(Refine)	(LW)
0.5	5.9769	5.7326	5.7328
1	1.9376	1.8464	1.848
1.5	1.1788	1.09	1.0884
2	0.9072	0.8048	0.8049
2.5	0.7786	0.6627	0.6626
3	0.7069	0.5804	0.5792
5	0.5954	0.4494	0.4493

and

Table 4. Comparison of non-dimensional fundamental frequencies $\overline{\omega }$ of (0°/90°/core/0°/90°) sandwich plate with $\frac{a}{\mathrm{b}}=1 \mathrm{and} \frac{h_3}{h_1}=10$.

$\frac{{\mathit{h}}_3}{{\mathit{h}}_1}$	Present Formulation	Rao	Rao
	(LW)	(Refine)	(LW)
4	2.2691	1.9084	1.9081
10	1.9376	1.848	1.848
20	2.153	2.1307	2.1311
30	2.3426	2.3321	2.3322
40	2.4756	2.469	2.469
50	2.5707	2.5658	2.5662
100	2.7899	2.7875	2.7874

, respectively. Natural frequencies were normalized using the relation $\overline{\omega }=\frac{\omega a^2}{h}\sqrt{\frac{\rho }{E_2}}$. According to the results in Table 2, for $\frac{a}{h}\succ 4$ (approaching a thin plate) the layer-wise model predicts similar values for $\overline{\omega }$ as those predicted by Rao et al. in [30].

The results in [30] were obtained based on an exact solution using the propagator matrix method in conjunction with the layer-wise model. According to Table 3, as the aspect ratio a/b increases (the composite plate approaches a long plate) the differences in results based on the present formulation and those in [30] become larger.

Table 4 shows the effect of core thickness to cover sheet thickness h₃/h₁. According to these results, the frequency ratio $\overline{\omega }$ is more susceptible to the aspect ratio a/b in comparison with h₃/h₁.

3.3. Five-Layer FGM Sandwich Plate

Now, we consider a five-layer square FGM sandwich plate resting on the Winkler foundation, as illustrated in Figure 1. The simply supported plate is assumed to be symmetric with respect to its mid-layer and all layers are experiencing a steady temperature ΔT. Each face sheet is made of functional materials while two vinyl ester (VE) based structural adhesive layers are used to bond the face sheets to an elastomeric core (Ellastollan R3000).

The advantages of vinyl ester resins include high stiffness and tensile strength, good chemical resistance, low cost, process versatility, and fast curing [31]. They possess good characteristics similar to epoxy resins, as well as unsaturated polyester resins. The shape of the stress–strain curves for a vinyl ester is strongly strain-rate dependent. Similar to strain rate effects, the shape of the stress–strain curve changes dramatically with any change in temperature. In this study, the deformation behavior of a vinyl ester polymer at strain rates 0.001/s and a wide range of temperatures (from room temperature (RT) to 100 °C) was investigated. As an illustration, the effect of temperature on the stress–strain behavior of a vinyl ester polymer at a strain rate of 0.001/s is shown in Figure 3 for different temperatures. At temperatures close to Tg (glass transition temperature) and high strains, viscoplasticity can be more pronounced and therefore, its effect on deformational behavior becomes more pronounced. However, viscoplastic deformation is not considered in the standard linear solid model of the material.

Table 5. Modules of vinyl ester at a wide range of temperatures. Data are extracted from Figure 3.

Elastic Modulus of Vinyl Ester (GPa)	Temperature (°C)
3.4	RT
3.13	50
2.8	75
2.5	100

shows the modules of elasticity of vinyl ester at a wide range of temperatures (from room temperature (RT) to 100 °C) [31].

Table 5 shows the elastic modulus of vinyl ester at different temperatures. Additionally,

Table 6. Modules of Elastollan R3000 at wide range of temperatures are extracted from Figure 5 in Ref. [32].

Elastic Modulus of Elastollan R3000 (GPa) ${\mathit{E}}_{\mathit{c}\mathit{o}\mathit{r}\mathit{e}}\equiv {\mathit{E}}_3$	Temperature (°C) T
2.8	RT
1.94	50
1.75	75
1.52	100

represents the modulus of elasticity of the Elastollan R3000 core at four different temperatures (from room temperature (RT) to 100 °C) [32]. Other mechanical properties of the layers used in this analysis are given in

Table 7. The mechanical properties of the five-layer sandwich plate.

Contituent	Mechanical Properties
(Elastollan R3000) core [32]	$v=0.45 , \alpha =20\times {10}^{-6}\frac{1}{℃},\rho =1380\frac{\mathrm{kg}}{{\mathrm{m}}^3}$
$\mathrm{Face} \mathrm{sheet} (\mathrm{AL}-{\mathrm{AL}}_2{\mathrm{O}}_3$) [1]	$E_m=70 GPa,E_C=380 GPa,v_C=v_m=0.3$ ${\alpha }_m=23.6\times {10}^{-6}\frac{1}{℃} , {\alpha }_c=6.6\times {10}^{-6}\frac{1}{℃}$ ${\rho }_m=2702\frac{\mathrm{kg}}{{\mathrm{m}}^3};{\rho }_c=3960\frac{\mathrm{kg}}{{\mathrm{m}}^3}$
Vinyl ester [33]	$v=0.375 , \alpha =50.8\times {10}^{-6}\frac{1}{℃},\rho =1106\frac{\mathrm{kg}}{{\mathrm{m}}^3}$

. The dimension of the total thickness for the sandwich plate is 12 mm and $\frac{h_1}{h}=0.1, \frac{h_2}{h}=0.02, \frac{h_3}{h}=0.76,h_5=1.2 \mathrm{mm} \left(-0.6 \mathrm{mm}\le z_1,z_5\le 0.6 \mathrm{mm}\right), h_2,h_4=0.24 \mathrm{mm}\left(-0.12 \mathrm{mm}\le z_2,z_4\le 0.12 \mathrm{mm}\right), h_3=9.12 \mathrm{mm}$.

To perform the finite element analysis, the whole plate was modeled using ABAQUS software. Three-dimensional solid element C3D20R was used for meshing. The simply supported boundary conditions were imposed on all four edges. Furthermore, to properly model the face sheets and invoke the properties of the FG material, the thickness of each cover sheet was divided into 30 thin layers with different properties defined according to Equation (1). The whole plate was subjected to thermal loading (all layers are under the same steady-state temperature T) and the sandwich plate was assumed to be resting on a Winkler elastic foundation. Finite element results on the vibrational modes are presented in

Table 8. Comparison of frequencies (HZ) based on LT and those FE findings for various (a/h). $n=1 ,\frac{a}{b}=1,\frac{h_1}{h}=0.1 \frac{h_2}{h}=0.02 , \frac{h_3}{h}=0.76,k_w=0,T_o=23 ℃, K=0.616$.

$\frac{\mathit{a}}{\mathit{h}}$	Natural Frequencies (HZ)
	Mode No.	1	2	3	4	5
10	FE	3251.2	5866.9	5866.9	7703	8583.3
	LW	3127.4	8378.2	5378.2	7003.4	7926.5
	% Difference	3.8	8.3	8.3	9.1	7.7
20	FE	1190.7	2520.6	2520.6	3473.3	4000.2
	LW	1195.4	2311.2	2311.2	3127.4	3586
	% Difference	0.4	8.3	8.3	9.9	10.4
30	FE	607.7	1387.5	1387.5	1996.5	2416.3
	LW	613.8	1296.2	1296.2	1829.6	2135.4
	% Difference	1	6.5	6.5	8.4	11.6

,

Table 9. Comparison of frequencies (HZ) based on LT and the FE findings for different n. $\frac{a}{h}=20, \frac{a}{b}=1 , \frac{h_1}{h}=0.1, \frac{h_2}{h}=0.02 , \frac{h_3}{h}=0.76, h=12 \mathrm{mm}, K_w=0, K=0.616, T_o=23℃$.

$\mathit{n}$
	Mode No.	1	2	3	4	5
0.5	FE	1244.4	2580.5	2580.5	3525.7	4031.6
	LW	1236.7	2339.8	2339.8	3139.2	3588.2
	% Difference	0.6	9.3	9.3	10.9	11
1	FE	1190.7	2520.6	2520.6	3473.3	4000.2
	LW	1195.4	2311.2	2311.2	3127.4	3586
	% Difference	0.4	8.3	8.3	9.9	10.4
2	FE	1115	2419.5	2419.5	3372.9	3923.5
	LW	1130.6	2248.1	2248.1	3079.2	3547.8
	% Difference	1.4	7.1	7.1	8.7	9.6

,

Table 10. Comparison of plate frequencies (HZ) based on the LT and those FE findings for different elastic foundations. $n=1, \frac{a}{h}=20, \frac{a}{b}=1 , \frac{h_1}{h}=0.1, \frac{h_2}{h}=0.02 , \frac{h_3}{h}=0.76, h=12 \mathrm{mm}, T_o=23℃, K=0.616$.

	${\mathit{k}}_{\mathit{w}}=0 \mathbf{N}/{\mathbf{m}}^2$			${\mathit{k}}_{\mathit{w}}={10}^8 \mathbf{N}/{\mathbf{m}}^2$			${\mathit{k}}_{\mathit{w}}={10}^9 \mathbf{N}/{\mathbf{m}}^2$			${\mathit{k}}_{\mathit{w}}={10}^{10}\mathbf{N}/{\mathbf{m}}^2$
Mode No	FE	LW	% Difference	FE	LW	% Difference	FE	LW	%Difference	FE	LW	% Difference
1	1190.7	1195.4	0.4	1239.9	1244.7	0.4	1615.9	1621	0.3	3649.1	3662	0.4
2	2520.6	2311.2	8.3	2544.1	2337.6	8.1	2746.8	2557.9	6.9	4270	4162.8	2.5
3	2520.6	2311.2	8.3	2544.1	2337.6	8.1	2746.8	2557.9	6.9	4270	4162.8	2.5
4	3473.3	3127.4	9.9	3490.4	3147.4	9.8	3640.7	3314.3	9	4892.5	4666.3	4.6
5	4000.2	3586	10.4	4015	3603.6	10.2	4146.5	3750.4	9.6	5280.5	4985.6	5.6

,

Table 11. Comparison of frequencies (HZ) based on LT and those FE findings for different elastic foundations. $n=1, \frac{a}{h}=20,\frac{a}{b}=1 ,\frac{h_1}{h}=0.1 \frac{h_2}{h}=0.02 , \frac{h_3}{h}=0.76,h=12 \mathrm{mm}, K=0.616,T_o=100 ℃$.

	${\mathit{k}}_{\mathit{w}}=0 \mathbf{N}/{\mathbf{m}}^2$			${\mathit{k}}_{\mathit{w}}={10}^8 \mathbf{N}/{\mathbf{m}}^2$			${\mathit{k}}_{\mathit{w}}={10}^9 \mathbf{N}/{\mathbf{m}}^2$			${\mathit{k}}_{\mathit{w}}={10}^{10} \mathbf{N}/{\mathbf{m}}^2$
Mode No	FE	LW	% Difference	FE	LW	% Difference	FE	LW	% Difference	FE	LW	% Difference
1	1047.4	1021.4	2.5	1103	1078.5	2.2	1513.9	1497.4	1.1	3603.7	3610.1	0.2
2	2083.5	1860.3	10.7	2112	1892.3	10.4	2352.8	2158.7	8.2	4028.5	3931.2	2.4
3	2083.5	1860.3	10.7	2112	1892.3	10.4	2352.8	2158.7	8.2	4028.5	3931.2	2.4
4	2803.8	2461.6	12.2	2825	2485.9	12	3009.4	2694.3	9.1	4444.1	4249	4.4
5	3162.8	2799.5	11.5	3181.6	2820.9	11.3	3346.5	3006.2	10.2	4679.3	4453.4	4.8

,

Table 12. Comparison of frequencies (HZ) based on LT and those FE findings. $n=1,\frac{a}{h}=20 ,\frac{a}{b}=1,\frac{h_1}{h}=0.1 \frac{h_2}{h}=0.02 , \frac{h_3}{h}=0.76,h=12 \mathrm{mm},k_w=0, K=0.616$.

	${\mathit{T}}_{\mathit{o}}=23 ℃$			${\mathit{T}}_{\mathit{o}}=50 ℃$			${\mathit{T}}_{\mathit{o}}=75 ℃$			${\mathit{T}}_{\mathit{o}}=100 ℃$
Mode No	FE	LW	% Difference	FE	LW	% Difference	FE	LW	% Difference	FE	LW	% Difference
1	1190.7	1195.4	0.4	1107.7	1094.2	1.2	1082.4	1063.7	1.8	1047.4	1021.4	2.5
2	2520.6	2311.2	8.3	2260.4	2039.1	9.8	2185.1	1962.2	10.2	2083.5	1860.3	10.7
3	2520.6	2311.2	8.3	2260.4	2039.1	9.8	2185.1	1962.2	10.2	2083.5	1860.3	10.7
4	3473.3	3127.4	9.9	3069.1	2719.9	11.4	2955.4	2608	11.8	2803.8	2461.6	12.2
5	4000.2	3586	10.4	3488.5	3102	11	3348	2970.5	11.3	3162.8	2799.5	11.5

,

Table 13. Comparison of frequencies (HZ) based on LT and FE findings (increasing the adhesive thickness). $n=1,\frac{a}{h}=20 ,\frac{a}{b}=1,\frac{h_1}{h}=0.1 , \frac{h_3}{h}=0.76,h=12 \mathrm{mm},k_w=0,T_o=23 ℃, K=0.616$.

	Mode No.
$\frac{{\mathit{h}}_2}{\mathit{h}}$		1	2	3
0	FE	1164.1	2475.7	2475.7
	LW	1168.8	2270.8	2270.8
	% Difference	0.4	8.3	8.3
0.01	FE	1177.6	2498.5	2498.5
	LW	1182.4	2291.6	2291.6
	% Difference	0.4	8.3	0.3
0.02	FE	1190.7	2520.6	2520.6
	LW	1195.6	2311.2	2311.2
	% Difference	0.4	8.3	8.3
0.04	FE	1216.2	2562.8	2562.8
	LW	1221.2	2350.7	2350.7
	% Difference	0.4	8.3	8.3
0.06	FE	1240.4	2602.5	2602.5
	LW	1245.7	2387.7	2387.7
	% Difference	0.4	8.3	8.3
0.08	FE	1269.1	2637.4	2637.4
	LW	1269.3	2422.9	2422.9
	% Difference	0.016	8.1	8.1

,

Table 14. The effect of any deviation in core elastic modulus from that of Elastollan^® R 3000. $n=1 ,\frac{h_1}{h}=0.1, \frac{h_2}{h}=0.02 , \frac{h_3}{h}=0.76, h=12 \mathrm{m},k_w=0, T_o=23 ℃, K=0.616$.

$\frac{{\mathit{E}}_3}{{\mathit{E}}_{\mathit{e}}}$	Method	Mode No.
		1	2	3	4	5
0.5	FE	1004.8	1989	1989	2672.9	3016.1
	LW	999.3	1808.3	1808.3	2387.7	2713.6
	% Difference	0.5	9.1	9.1	10.7	10
1	FE	1190.7	2520.6	2520.6	3473.3	4000.2
	LW	1195.4	2311.2	2311.2	3127.4	3586
	% Difference	0.4	8.3	8.3	9.9	10.4
2	FE	1296.2	2929	2929	4178.3	4961.6
	LW	1355.1	2810.2	2810.2	3928.7	4565.6
	% Difference	4.5	4.1	4.1	6	8
4	FE	1404.3	3311.3	3311.3	4874.7	5930.8
	LW	1474.3	3243.1	3243.1	4689.6	5535
	% Difference	5	2.1	2.1	4	6.7

and

Table 15. The effect of any deviation in the adhesive elastic modulus from that of vinyl ester. $n=1 ,\frac{h_1}{h}=0.1 \frac{h_2}{h}=0.02 , \frac{h_3}{h}=0.76,h=12 \mathrm{mm},k_w=0,T_o=23 ℃, K=0.616$.

$\frac{{\mathit{E}}_2}{{\mathit{E}}_{\mathit{v}\mathit{i}\mathit{n}}}$	Method	Mode No.
		1	2	3	4	5
0.5	FE	1183.4	2496	2496	3434.2	3935.1
	LW	1185.1	2282.4	2282.4	3083.2	3532.9
	% Difference	1.4	8.6	8.6	10.2	10.2
1	FE	1190.7	2520.6	2520.6	3473.3	4000.2
	LW	1195.4	2311.2	2311.2	3127.4	3586
	% Difference	0.4	8.3	8.3	9.9	10.4
2	FE	1195.1	2533.8	2533.8	3494	4024.8
	LW	1201.4	2327.2	2327.2	3151.7	3615.1
	% Difference	0.5	8.2	8.2	9.8	10.4
4	FE	1198.4	2541.8	2541.8	3505.6	4038
	LW	1205.1	2335.6	2335.6	3163.9	3629.5
	% Difference	0.6	8.1	8.1	9.7	10.1
6	FE	1200.5	2545.5	2545.5	3510.5	4043.1
	LW	1207.1	2338.9	2338.9	3168.2	2634.5
	% Difference	0.5	8.1	8.1	9.8	10.1

and graphically shown in Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9 followed by subsequent discussions.

Table 8 shows a comparison in results for the first five natural frequencies of the square FGM sandwich plate using the present theory (LW) and finite element analysis (FE) for various side-to-thickness ratios (a/h). According to these results, any increase in (a/h) ratio decreases the natural frequencies. As can be observed from this table, not only for thin plates but also for thick plates, the fundamental frequencies are predicted as accurately by the present method with those obtained in FE.

In Figure 4 and Figure 5, the first four mode shapes of the square sandwich plate (a/h = 20) are presented.

The effect of the variation in the power-law index on $\overline{\omega }$ for modes one to five are presented in Table 9. The presented results that are based on the FE solution, as well as the semi-analytical solution, are obtained based on values of Young’s modulus (E), mass density $\left(\rho \right)$ and the coefficient of thermal expansion $\left(\alpha \right)$ obeying the power-law index distribution given in Equation (2). According to these results, for values of n $>1$ (metallic phase dominating the ceramic), the sandwich plate experiences lower frequency values, while for n $<$ 1 (a sandwich plate with more ceramic in the core), the opposite behavior is observed due to the higher stiffness of the ceramic phase.

The natural frequencies of the five-layered square sandwich plate resting on elastic foundations are presented in Table 10 and Table 11. It is also observed for plates resting on elastic foundations that the natural frequencies are increased as the foundation stiffness increases. The results from the semi-analytical solution that are also supported by the finite element analysis show that the increasing trend in $\overline{\omega }$ is much higher for the higher values of k_w. According to the results in Table 12, the effects of a rise in temperature from $23{\hspace{0.17em}}^{\circ }\mathrm{C}$ to 100 ${}^{\circ }\mathrm{C}$ appear to have an adverse effect on the frequency ratio $\overline{\omega }$ associated with the first mode of vibration, while for the other modes, the opposite behavior is observed.

The results in Table 13 show the effect of any increase in the adhesive thickness on the natural frequencies of the FGM sandwich plate in the absence of the elastic foundation at room temperature. Although the elastic modulus of the adhesive layer is much smaller than that of the face sheets and the core, any minute increase in the thickness of this layer from $\frac{h_2}{h}=0.01$ (h₂ = 0.12 mm) to $\frac{h_2}{h}=0.08$ (h₂ = 0.96 mm) introduces a 9% increase in the natural frequencies of the sandwich plate, the magnitude of which cannot be neglected.

The effect of core thickness on the natural frequencies of the sandwich plate is shown in Figure 6 and Figure 7 for the two cases with and without elastic foundation. The results in Figure 6 indicate that in the absence of an elastic foundation (k_w = 0) for values of a × b ranging from 0.24 × 0.24 m² to 3.0 × 3.0 m², the increase in the core thickness up to a value of ≈70 mm shows its direct effect on raising the natural frequencies of the plate, while beyond this value, this effect can be neglected.

However, for plates resting on an elastic foundation with a value of k_w = 10⁸ N/m², this behavior is observed to exist only for values of $a\times b\prec 1.0\times 1.0\hspace{0.17em}\hspace{0.17em}{m}^2$. For larger plates (values of $a\times b\ge 1.0\times 1.0\hspace{0.17em}\hspace{0.17em}{m}^2$), the presence of an elastic foundation with a value of k_w = 10⁸ N/m² reduces the plate frequency as the composite plate gets thicker in the core.

Table 14 and Table 15 show the effects of $\frac{E_3}{E_e}$ and $\frac{E_2}{E_{vin}}$ on natural frequencies of the sandwich plate at room temperature and the absence of an elastic foundation. In Table 14, Ef corresponds to the elastic modulus of the vinyl ester core while E3 represents any deviations in this modulus that presents a new core material with the same density. Additionally, in Table 15, E2 corresponds to the elastic modulus of the adhesive material (which may be different from that of vinyl ester), and Ee is the elastic modulus of the Elastollan^® R 3000 selected as the core material. According to the results in Table 14, increasing the elastic modulus ratio of the core from 0.5 to 4.0 increases the mode I natural frequency by 46.7%. This increase can be interpreted as the rise in the elastic modulus of the core with respect to the Elastollan^® R 3000, which makes the sandwich plate stiffer with the assumption that the core density does not change. This increase is more pronounced at higher modes of natural frequencies. Moreover, according to Table 15, as the elastic modulus of the adhesive layer advances that of the Elastollan^® R 3000, a similar effect is observed in natural frequencies. Therefore, to have better control of the vibrational frequencies of a sandwich plate, it is essential to select the right elastic moduli for the core and adhesive layer.

Figure 8 shows the variation in the fundamental natural frequency versus the elastic modulus of the core in a sandwich plate with different dimensions resting on an elastic foundation with $k_w={10}^8 \frac{N}{m^2}$. The results show that in sandwich plates with small aspect ratios (i.e., a × b = 0.3 × 0.3 m², a thicker plate compared with other shown dimensions), the elastic modulus of the core considerably increases the fundamental natural frequency. This is partly due to higher plate stiffness attributed to the higher elastic modulus of the core. For larger plate aspect ratios (thinner plates), the presence of an elastic foundation with the given value for k_w seems to be ineffective on the natural frequencies.

However, according to Figure 9, the variations in the first three modes of natural frequencies of the sandwich plate, in the presence of a Winkler elastic foundation, seem to highly depend on small values of aspect ratio a/b. However, for values of $a/b\succ 4$, this effect is negligible.

As shown in the previous investigations by other authors [34,35], the use of the layer-wise method in conjunction with the first-order deformation theory produces good results on a few stress components, while on others, modified theories must be implemented to improve the accuracy of other stress components. However, this makes the problem more complicated. For example, the study performed by Raissi et al. [35] showed that modifying the shear deformation theory can improve some of the inaccuracies in previous findings obtained on a few stress components that resulted from the application of lower-order shear deformation theories. Consequently, the consideration of higher deformation theories in conjunction with the layer-wise theory is being considered by current authors in their future investigations for possible improvement of the present results.

The difference between a semi-analytical solution and for those finite element solutions is larger for higher modes because the element shape functions provide a better basis (approximation) of the shape of low modes and a relatively poorer basis for higher modes.

4. Conclusions

In this work, the effects of an adhesive layer bonding the core of functionally graded (FG) cover sheets were investigated using the free vibration of a five-layer sandwich composite plate resting on a Winkler elastic foundation. The whole plate was assumed to experience a steady-state temperature. On application of Hamilton’s principle, the equilibrium equations were derived using the layer-wise model and numerically solved using a MATLAB software program. A finite element model was also prepared and solved to support the results of the semi-analytical solution. The results showed that in free vibration, the first-order shear deformation theory provides good accuracy in predicting the lower modes of the fundamental frequencies of the sandwich plate. However, at higher frequencies, the semi-analytical results were slightly reduced in accuracy in comparison with the finite element findings. The results indicate the presence of a higher temperature aggravates the adverse effect of the plate frequencies, especially at higher modes. In addition, the presence of an adhesive layer with a thickness of h2 = 0.96 mm introduces a 9% increase in the natural frequencies of the sandwich plate at mode I compared with a three-layer (an inside core encapsulated by two cover sheets) sandwich plate. The accuracy of this theory in analyzing the frequencies of an FGM sandwich plate increases with the presence of an elastic foundation. In addition, the results indicate that the presence of adhesive layers does not have a significant effect on the frequencies. Furthermore, by decreasing the aspect ratio, the natural frequencies of the sandwich plate are increased where the square sandwich plate has the most natural frequencies.