Three-Dimensional Free Vibration Analysis of Thermally Loaded FGM Sandwich Plates

Using the finite element code ABAQUS and the user-defined material utilities UMAT and UMATHT, a solid brick graded finite element is developed for three-dimensional (3D) modeling of free vibrations of thermally loaded functionally gradient material (FGM) sandwich plates. The mechanical and thermal material properties of the FGM sandwich plates are assumed to vary gradually in the thickness direction, according to a power-law fraction distribution. Benchmark problems are firstly considered to assess the performance and accuracy of the proposed 3D graded finite element. Comparisons with the reference solutions revealed high efficiency and good capabilities of the developed element for the 3D simulations of thermomechanical and vibration responses of FGM sandwich plates. Some parametric studies are carried out for the frequency analysis by varying the volume fraction profile and the temperature distribution across the plate thickness.


Introduction
Sandwich panels are usually used instead of traditional structural elements made of metals and alloys, when increased strength and stiffness with little resultant weight are required for engineering applications [1][2][3]. Although sandwich panels provide outstanding structural features, this structural design has also drawbacks. A typical sandwich panel configuration has a high mismatch of material and geometrical properties between the face sheet and the core [4]. Due to this, a variation of the interfacial stresses induced by thermal or/and mechanical loads is significant at the face sheet-to-core interface [5][6][7]. Therefore, the performance and reliability of such tri-layer composites are eventually defined by the quality of the bonded interface [8,9]. When debonding arises between the skin and core material layers, sandwich panels significantly lose their load bearing capacity [10,11]. The modal dynamic characteristics of such panels damaged by debonding are changed [12][13][14][15] and their overall dynamic responses are modified [16][17][18][19] as well. Moreover, the debonding may cause eventual failure of the sandwich panels [20][21][22] under dynamic loads.
In regards to reducing or avoiding the debonding issue, functionally graded materials (FGMs), with mechanical and thermal properties that are smoothly distributed over the volume, have potential for use as basic layers in sandwich panels. Usually, such sandwich configuration is achieved by gradually changing the volume fraction of the FGM constituents across the sandwich plate thickness from the bottom face sheet to the top face sheet [23]. This removes the interface stress concentration and allows controlling deformation, dynamic responses, and etc, by customizing the gradation profile [24][25][26]. Another characteristic of FGM sandwich panels, with metal and ceramic material, is that they can be used at elevated service temperatures [27,28]. This fact, however, gives rise to a new demand for providing safe operation of FGM sandwich panels in thermal environments. Therefore, an accurate description of the thermomechanical behavior of FGM sandwich plates becomes mandatory, mainly, to prevent their thermal failure, as considered for FGM coatings in [29][30][31][32][33], amongst the most recent studies. At the same time, the analysis of free vibrations of FGM sandwich panels subjected to temperature, in the design stage, is important for estimating their overall performance.
To date, a considerable amount of research on the behavior of FGM plates and shells under temperature loading is found in the literature. Various analytical (or semi-analytical) and numerical methods have been developed for this, as recently reviewed by [34] reporting that there is intense research activity with respect to modeling thermally loaded FGM composites because reliable and practical analysis procedures, in terms of accuracy and computational efforts, are in high demand. In this regard, the development of two-dimensional (2D) models is motivated by computational efficiency, e.g., [35][36][37][38][39]. On the other hand, such models partially lose the accuracy of predictions due to simplifications caused by prescribing the behavior of shear deformations across the panel thickness. So-called quasi-three-dimensional theories, which adopt assumptions for both shear and normal deformations in the thickness direction, have recently been proposed as an improvement to 2D theories for more reliable analyses of FGM sandwich panels, for example, in [40][41][42]. Nevertheless, only three-dimensional (3D) models, which are not spoiled by any additional prescriptions inherent to 2D or quasi-3D theories, are able to address unique aspects stemming from the complex dynamic response of FGM sandwich panels [43]. However, there are two main obstacles to the use of 3D models for predictions of FGM plates. First, analytical exact 3D solutions for FGM sandwich plates are only available for simple cases of boundary conditions and geometries [44][45][46][47][48]. Secondly, although the finite element method (FEM), which is one of the leading computational tools, can solve problems associated with the complex geometry and different boundary conditions, the efficiency of conventional 3D finite element models is not suitable for real-scale FGM sandwich structures. The main reason for this is the layered approach for modeling material gradient with conventional finite elements because the conventional finite element has constant elastic properties over its domain. To overcome this obstacle inherent in models with conventional finite elements, the graded finite element, which assigns the material gradation profile at the element level, must be used instead. There are two techniques to elaborate a graded finite element. A nodal approximation of element material properties with interpolation functions identical to those utilized for the displacement field is used in the first technique, for example, as done in [49,50]. The second technique is based on sampling the material properties directly at the integration points of the element [51].
Not to mention fully coupled thermomechanical problems of FGM sandwich plates involving nonlinear material behaviors, the modal dynamics of FGM sandwich plates in thermal environment is strongly dependent on the distribution of thermal stresses across the plate thickness [52][53][54][55]. On the other hand, the thermomechanical behavior of such plates is defined by the material gradation profiles. Hence, a prerequisite for high-fidelity modeling is a precise thermoelastic analysis of such plates with accurately prescribed FGM properties. It is well-recognized that the FEM is a powerful means for solving multiphysical thermomechanical problems. Moreover, the method has been implemented in a series of commercially available codes, for instance, ABAQUS [56] which is popular among researchers and engineers. However, the package does not provide finite elements with a spatial variation of material properties. To accomplish this, programming of externally prescribed subroutines within the package environment is required. Some existing works suggest the use of the ABAQUS UEL subroutine which implements a material gradation in the FGM plates via a user-developed finite element, for example, [57]. Although this approach gives high flexibility in modeling, it requires the knowledge of an experienced user and extensive benchmarks for the performance of the element before simulations. Another approach for the implementation of varying material properties has been reported for ABAQUS 2D plane strain elements in [30,58,59], where the strain-stress state of FGM pavement and the thermomechanical behavior of the FGM plate have both been analyzed. The necessary material properties of studied functionally graded materials have been distributed at the Gauss points by coding appropriate material user-defined subroutines such as UMAT and UMATHT. In addition, this modeling technique has been extended to 2D plate/shell and 3D models of FGM plates in [60,61], respectively, however, only the static bending analysis has been simulated there. Recently, three-dimensional finite element models of FGM sandwich plates for dynamic modal analyses have been developed in [62,63].
The aim of this study is to propose an efficient approach for implementing a 3D graded finite element into ABAQUS code to perform a computationally accurate free vibration analysis of thermally loaded FGM sandwich plates. A novel graded finite element has been developed based on the 3D brick graded finite element proposed for the free vibration analysis of 3D FGM sandwich plates in [62,63]. In our study, we extend the functionality of the element to use it for modeling FGM sandwich plates under thermal loading. First, we consider a thermomechanical analysis to compute through-the-thickness distributions of displacements and stresses in the sandwich plates. Then, with a known thermally induced stress state, the free vibration analysis is carried out. The 3D brick graded element has been developed by coding a combination of the subroutines such as UMAT, UMATHT, and USDFLD similar to the 2D thermomechanical finite element analysis presented in [59]. The performance of the proposed 3D graded element has been demonstrated by the 3D modelling of heat transfer and free vibrations of FGM sandwich plates subjected to thermal loading. The accuracy of the graded element has been validated by comparison with results available in the literature for FGM sandwich plates. Parametric studies have also been carried out to determine the effect of varying volume fraction profiles and the temperature on natural frequencies and associated mode shapes. We believe that the results of this research can be used as a benchmarks for 2D solutions and results obtained by other numerical methods.

Problem Formulation
For the sake of completeness, the thermomechanical problem for a continuum made of a functionally graded material is briefly summarized in this section. Throughout the section we adopt the usual notations used in most books on continuum mechanics, which can be referred to for more details.

Thermomechanical Problem
Let us consider the FGM sandwich plate as a 3D deformable medium occupying the domain Ω ∈ [0, a] × [0, b] × − h 2 , + h 2 bonded by the surface, ∂Ω ⊂ Ω, at an instant of time, t ∈ [0, t end ] The plate is defined at a given temperature T 0 in the unstressed reference configuration with respect to a rectangular Cartesian co-ordinate system, x i = (x, y, z), with the z-axis aligned along the plate thickness and with the plane, z = 0, coinciding with the mid-plane of the sandwich plate. In addition, the planes z = ±h/2 refer to the bottom ∂Ω − and the top ∂Ω + plate surfaces, respectively, where ∂Ω\∂Ω − ∪ ∂Ω + = − h 2 , + h 2 , as shown in Figure 1a. In the Lagrangian description, the equations of mechanical motion and thermal equilibrium at each spatial point, x, of the domain, Ω, at a time instant, t, in the absent of body forces and internal heat sources can be presented as [64]: where, σ and q are the Cauchy stress tensor and the heat flux vector, respectively; u is a displacement field and θ stands for a temperature field associated with a change of the instantaneous temperature T(x, t) above the reference temperature T 0 at time t ∈ [0, t end ]; and ρ(x), c(x), and β(x) denote the mass density, the specific heat, and the stress-temperature modulus, respectively, which are functions of a spatial position x. Assuming small displacements and deformations, the infinitesimal strain tensor as a sum of elastic mechanical "el" and thermal "th" parts is expressed by where, α(x) is the coefficient of thermal expansion and I is the identity tensor. Consider the FGM is a linear isotropic material that complies with the classical law of thermal conductivity. Then, the thermoelastic constitutive equations of the FGM sandwich plate have a form: where, the Lamé constants λ(x) and µ(x) of the elasticity tensor D(x), the modulus β(x) = α(3λ + 2µ) and the conductivity κ(x) are pointwise functions of location. Two types of the boundary conditions must be specified on the plate surface ∂Ω at any instant in time. The mechanical boundary conditions are prescribed by displacements u on the boundary ∂Ω u and tractiont on the boundary ∂Ω t , where, ∂Ω u ∪ ∂Ω t = ∂Ω. In a similar manner, the thermal boundary conditions are defined by a prescribed temperature T and heat flux q 1 and/or an exposure to an ambient temperature through convection so that q 2 =ĥ(x)(T − T ∞ ) on the plate surfaces ∂Ω θ and ∂Ω q = ∂Ω q 1 ∪ ∂Ω q 2 , respectively. Here,ĥ(x) is the heat film transfer coefficient and T ∞ is the temperature of the surrounding medium [64].
Using the principle of virtual work and collecting (1) to (3) with appropriate boundary conditions, the system of mechanical and energy equations can be rewritten in the weak form as follows: for all kinematically admissible virtual displacement δu and temperature δθ fields.

Properties of FGM
We assume that the sandwich plate is made of a two-phase metal-ceramic functionally graded material and, also, without loss of generality, a smooth variation of material thermomechanical properties across the plate thickness only, i.e., in the z-direction, see Figure 1b. Herewith, it is deemed that the face sheets (or skins) are homogeneous, i.e., pure metal on one side and pure ceramic on the other one, with a small or negligible small thickness as compared with the thick metal-ceramic FGM core. In addition, we suggest that the gradation profile of the ceramic volume fraction from the bottom to top sandwich plate skins is known and is determined by a power-law function in the form: where, V − c and V + c are the volume fraction of ceramic on the bottom and top surfaces, respectively. The case of V − c = 0 and V + c = 1 refers to the gradation profile from pure metal on z = −h/2 to pure ceramic on z = +h/2. It follows from (5) that the FGM plate is ceramic-rich when the parameter p < 1, and metal-rich when the parameter p > 1. Figure 1c shows the volume fraction variation of the ceramic phase along the plate thickness depending on the values of the power-law index p.
The effective mass density, and the thermal and mechanical properties at a point of the FGM are specified based on the "rule of mixture" as follows: where, P(z) represents either the mass density or any of the thermomechanical parameters; and the subscripts "m" and "c" are the metallic and ceramic phases whose volume fractions are such that In turn, each of the parameters may depend on the temperature in the form: where, P 0 stands for material parameters at the reference temperature T 0 and P −1 , P 1 , P 2 , and P 3 are constants specifying the temperature dependence of the material at the instantaneous temperature T. Finally, since the gradient in properties occurs only along the plate thickness direction, then the material tensor in (3) in the Voigt notation is as follows:

Method of Solution
A displacement-based FEM framework is used for solving the problem formulated in Section 2.

Finite Element Discretization
In the context of FEM, the actual continuous model of the sandwich plate is idealized by an assemblage of arbitrary non-overlapping finite elements, Ω = U N e=1 Ω e , interconnected at nodal points. In each base element the displacement vector, u (e) , and a scalar function of temperature, θ (e) , are approximated by suitable interpolation functions such that Here, the summation over all nodal points of the base element is intended, also, N = [N I (x)] and N = [N P (x)] are matrices of the shape functions N I and N P for the displacements and the temperature, respectively, associated with certain nodes I and P. The vectors U (e) and Θ (e) are the nodal unknown displacements and temperature at those nodes. The coupling between the mechanical and thermal problems is assumed to be due to temperature only, i.e., there is no feedback on the energy expression through the displacement field. This assumption is reasonable if a thermomechanical model is used that does not involve internal variables, such as plastic strains for computing the energy dissipation rate [56]. Substituting the displacement and temperature approximations (9) into the variational equalities (4) and accounting for the material laws in (3), we arrive at the system of semidiscrete equations of a one-way thermomechanical problem at the element level as follows: ..
Forms of the element matrices involved in (10) are found, for example, in [59]. The assembly operation (•) = A N e=1 (•) (e) over all the finite elements leads to the global system of semidiscrete equations for the thermoelastic problem in the form: where, M, K u , K θ and C are the usual global mass, stiffness, conductivity, and capacity matrices, respectively; K uθ is the coupling thermoelasticity matrix; U and Θ are the global vectors of nodal displacements and nodal temperature, the dots over them correspond to the time derivatives of these vectors; and F u and F θ are global vectors of the mechanical and thermal forces.
In the case of an uncoupled formulation, the temperature is given as an external load and it is not a primary variable in the mechanical analysis. More specifically, the thermal virtual work leads to the initial stress matrix known as a geometric stiffness matrix K G , which contains the terms due to the temperature loading on the leading diagonal. Thus, first, the temperature field is computed at the given thermal and displacement boundary conditions. Then, a temperature profile, known for solving the mechanical problem with the same displacement boundary conditions, is used to find the displacement and stress fields. The nonlinear finite element equations of the thermomechanical problem are solved by an iterative method, where the nonlinear terms of linearized equations are evaluated as a known solution from the preceding iteration. The Newton-Raphson iterative method is used in ABAQUS [56]. Finally, the frequency analysis, which accounts for the initial deformed state associated with the temperature-induced stresses is carried out by extracting natural frequencies from the eigenvalue-type equation: where, ω is an undamped circular frequency and φ is a vector associated with mode shape at a specific frequency ω.

Three-Dimensional Graded Element
The thermomechanical analysis and the modal frequency extraction procedure for FGM sandwich plates are carried out with the ABAQUS/Standard code using three-dimensional models. Since conventional 3D finite elements, which are available in the ABAQUS finite element library, are not able to model a variation of the thermal and mechanical properties within the element volume, a 3D graded finite element incorporating gradients of material properties is developed.
As mentioned in the Introduction, a 3D graded finite element has been developed for performing the modal frequency analysis of FGM sandwich plates in [62,63]. To incorporate a variation in the elastic properties of the heterogeneous material into the finite element, the material user-defined subroutine UMAT that establishes the tangent element stiffness matrix was programmed, while the average mass density value was adopted to determine the element mass matrix. In our study, the performance of the 3D graded finite element is extended to provide thermal loading and temperature-dependent material properties for carrying out the thermomechanical analysis of the FGM sandwich plates. Such an analysis requires computing and storing the internal thermal energy and the heat flux which comply with the energy balance equation (1) and Fourier's law of heat conduction (3), respectively, as well as modifying the mechanical behavior by accounting for thermal strains (2). The implementation of a spatial variation of the thermal and modified mechanical properties in the selected direction in the 3D graded finite element has been done following the procedure outlined in [59] for a 2D graded finite element. With this approach, a master 3D temperature-displacement finite element, either eight-node linear C3D8 or twenty-node quadratic C3D20 brick isoparametric element with either reduced or fully integration scheme, which are available in ABAQUS (Figure 2), has been supplemented by a combination of user-defined subroutines such as UMAT, UMATHT, and USDFLD, [56]. In addition, it should be mentioned that ABAQUS interpolates the temperature field using only the first-order approximation regardless of the order of approximation of the displacement field, as shown in Figure 2b,c. The material subroutine UMAT was programmed to define a through-the-thickness variation of the Lamé constants, λ(x) and µ(x), and the stress-temperature modulus, β(x), associated with the thermal expansion coefficient, α(x). The distributions across the thickness of the material thermal properties were incorporated into the element by coding the specific heat, c(x), the thermal conductivity, κ(x), and the film heat transfer,ĥ(x), coefficients in the UMATHT subroutine. Finally, the through-the-thickness variation of mass density was incorporated into the element using the USDFLD subroutine. Moreover, if the material parameters were deemed to exhibit a temperature dependence, appropriate relationships (7) with given coefficients for each material parameter were also programmed as functions of the temperature in the mentioned subroutines. In doing so, the instantaneous temperature, T, was known as it is a variable passed for information at each time increment in the subroutines, and therefore, the property was able to be computed at any current temperature value.
By running the ABAQUS code, the element matrices, presented in (10), which involve variations of the material thermal and mechanical properties coded in accordance with a certain relation for FGM constituents, have been generated by calling the corresponding user-defined subroutines. Thus, the material properties that account for the given material gradation profiles have been assigned directly at the Gauss integration points of the element (see Figure 2b,c). In such a way, any arbitrary material gradient, for example, a power-law distribution of the ceramic phase in the thickness direction of the plate can be prescribed. More details of the implementation of graded elements into the ABAQUS code are found in [30,59] and the ABAQUS manual [56]. In addition, it is important to note that the mass density, averaged over the FGM sandwich plate volume [62,63], was used in the free vibration analysis instead of its spatial representation in the case of thermomechanical analysis.

Comparison Study
The performance of the 3D graded element described above for solving thermomechanical and free vibration problems has been verified by comparing calculated numerical solutions with results available in the literature.

Cube Problem
As a first validation problem, a unit FGM cube (L = 1) that has been subjected to prescribed temperatures on two opposite sides and insulated in all the other sides is considered, Figure 3a. This problem has been solved analytically and numerically using the boundary element method in [65]. The top surface of the cube z = 1 is maintained at the temperature T L = 100 • C, while the temperature of the bottom surface at z = 0 is 0 • C. The reference temperature is also assumed to be 0 • C. The variations of thermal conductivity and specific heat along the z-axis are defined by the expressions: and The analytical solution for the temperature in the transient thermal analysis is known [65] as where, the coefficients, B n , are given in the form: In the finite element simulations, the cube is discretized with 4 × 4 × 4 twenty-node quadratic brick graded elements, as illustrated by Figure 3b. A transient heat transfer analysis is performed with ABAQUS, calling the material subroutines UMATHT and USDFLD. The temperature profile along the z−axis is plotted at different times for the given exponential material variations and compared with the analytical solutions, as shown in Figure 3c. It is evident from the plot that the numerical and analytical results are in excellent agreement.

Analysis of FGM Plates
As a second example, an aluminum-zirconia functionally graded square plate with sides a = b = 0.2 m and thickness h = 0.01 m, as studied in [35] is considered. The plate is assumed to be simply supported on all the edges and is exposed to a temperature field such that the ceramic-rich top surface is held at 300 • C and the metal-rich bottom surface is held at 20 • C. A stress-free state is assumed to be at 0 • C. The material thermomechanical parameters of the FGM plate are listed in Table 1. For the purpose of comparison, the values of the volume fraction exponent, p, in (5) have been accepted as 0.0 (pure ceramic), 0.2, 0.5, 1.0, 2.0 and ∞ (pure metal). Figure 4a shows the variation of the temperature through-the-thickness profiles of the FGM plate depending on the exponent values. By comparing the temperature profiles shown in Figure 4a with those illustrated in [35] (p. 680), one can conclude that the plots are nearly identical. With these calculations it was found that the calculated results converge to the reference data in [35] with an increase of the number of graded elements in the thickness direction. The best agreement between both the solutions was achieved using eight graded elements across the plate thickness. This resulted in time-consuming computations in the case of cubic elements used in the mesh. In order to speed up the computations, brick graded elements, with an aspect ratio 4:4:1, have been used instead. The elements provided more than 10 times faster computations with the same through-the-thickness profiles for variations of temperature and quite acceptable results for temperature-induced deflections and stresses across the thickness, as illustrated in Figure 4b-f, respectively, as compared with those in [44]. Therefore, such elements are used in the calculations to follow.  Figure 4b. This is related to the fact that the deflection depends on the product of the temperature and the thermal expansion coefficient. The latter is larger in the metal-rich region, while the temperature is higher in the ceramic-rich portion. As a result, the thermal strains are not uniform over the plate thickness and the responses of the FGM plates are not intermediate to those of the pure metal and ceramic plates. The temperature-induced through-the-thickness distributions of central longitudinal and transverse normal stresses, and a transverse shear stress at the center of plate edge and an in-plane shear stress at the plate corner, shown in Figure 4c-f, respectively, demonstrate that the longitudinal normal stresses in the FGM plates exhibit nonlinear profiles in contrast to linear ones in the pure metal and ceramic plates, whereas, the transverse normal stress and the shear stresses of all the plates have rather similar profiles and, in the case of FGM plates, the stress distributions crucially depend on the power-law index p.
In order to evaluate the accuracy of the developed graded element in the free vibration analysis of thermally loaded FGM plates, the natural frequencies of a fully clamped (CCCC) square FGM plate with the thickness-side ratio h/a = 0.1 have been computed and compared with those available in [46]. Steel-silicon SUS304/Si 3 Ni 4 functionally graded plates were considered. The properties of constituents of the FGM are given in Table 1, while the temperature-dependent constants of the constituents can be found in [46] (p. 737, Table 1). The FGM plates were assumed to be subjected to different temperatures equal to 300 K, 600 K and 800 K which are uniformly distributed across the plate thickness. The natural frequencies of the FGM plates extracted from the finite element analysis have been nondimensionalized as follows: where, I m = hρ m and D m = E m h 3 / 12 1 − ν 2 m are expressed using the appropriate values of the stainless steel at the reference temperature T 0 = 300 K. Table 2 shows a good agreement between the computed nondimensional frequencies (ω FEM ) and the results (ω Re f . ) reported in [46].

Parametric Study
After establishing the correctness of the developed 3D graded finite element, parametric studies are performed to investigate the effects of gradation profiles in thermo-elastic properties and temperature distributions on the free vibrations of FGM sandwich plates.
First, we consider the free vibration analysis of SUS304/Si 3 Ni 4 functionally graded square sandwich plates with the skins' thickness negligible as compared with the core thickness. The geometry and material properties of the sandwich plates were identical to the analysis for the FGM plate in Section 4.2. It is assumed that across the thickness, the sandwich plates may be subjected to either a uniform temperature field, T b = T t = T, or a temperature profile associated with a steady-state heat transfer due to differently prescribed temperatures on the bottom, T b , and top, T t , plate surfaces. Two types of boundary conditions, i.e., all edges simply supported (SSSS) and all edges clamped are used in the calculations. Five different material gradations defined by the power-law index p = 0.2, 0.5, 1, 5 and 10 as well as pure ceramic (p = 0) and metal ( p → ∞) homogeneous plates are examined. In Tables 3  and 4, the first ten nondimensional natural frequencies ω of simply supported and clamped FGM plates under the uniform temperature of T=600 K are presented, respectively. In addition, in Tables 3  and 4, the frequencies calculated with the 3D graded elements are compared with those available in [46] for the studied FGM plates. Similar to the previous study, the first ten nondimensional natural frequencies of simply supported and clamped FGM plates subjected to the temperature profiles following from the solution of the thermomechanical analysis under steady-state conditions with the prescribed temperature on the top (ceramic) surface T t = 600 K and at the reference temperature on the bottom (metal) surface T b = T 0 = 300 K are collected in Tables 5 and 6, respectively. The contour plot of the temperature distribution within the plate (a quarter of plate is removed from the presentation to illustrate the temperature distribution inside the plate) and the variations of temperature across the thickness (at the central section of the plate) depending on the power-law index p, which have been predicted by the thermomechanical analysis for the SSSS and CCCC plates, are illustrated in Figure 5a,b. It is evident from the latter plot that the effect of p is not significant for this material and the nonlinear temperature profiles are very close to the linear temperature distribution. Table 5. The first ten nondimensional frequencies of SSSS SUS304/Si 3 Ni 4 functionally graded square plates thermally loaded by the nonlinear temperature rise as shown in Figure 5.  Table 6. The first ten nondimensional frequencies of CCCC SUS304/Si 3 Ni 4 functionally graded square plates thermally loaded by the nonlinear temperature rise as shown in Figure 5. It is worth noting that there is very good agreement between the present results and the referenced solutions, as seen in Table 3 to Table 6. This demonstrates the accuracy and effectiveness of the 3D graded finite element developed in the present work. In addition, for the sake of clear demonstration of the effect of the material gradation profile on the natural frequencies, some frequencies from Tables 5  and 6 for the simply supported and clamped FGM plates subjected to the nonlinear temperature rise are plotted as functions of the volume fraction exponent p in Figure 6a,b, respectively. It is obvious from the plots that the frequencies decrease with an increasing in the power-law index, i.e., growing the percentage of ceramic fraction in the top thickness of the FGM plates. In doing so, the higher frequencies show more intensive descending trends.  Next, sandwich plates with thickness of skins h f = 0.1h and thickness of core h c = 0.8h, which are referred to as 1-8-1 sandwich configurations, are considered for the free vibration analysis. Three different boundary conditions, such as fully simply supported, fully clamped, and two edges simply supported and two edges clamped (SCSC) are examined. It is assumed that the sandwich plates have homogenous skins such that the top and bottom skins are pure ceramic and pure metal, respectively, and the core is SUS304/Si 3 Ni 4 functionally graded material identical to that in the previous study with the volume fraction exponent p = 0.2, 0.5, 2 and 10. The thermomechanical constants of the material constituents are listed in Table 1. The sandwich plates are assumed to be subjected to the prescribed temperature on the top surface T t , and the bottom surface is at the reference temperature, i.e., T b = T 0 . For the free vibration analyses, first, temperature gradients across the thickness of the plates are computed using the thermomechanical analysis under steady-state conditions. Three different temperatures T t = 300, 400 and 600 K are applied to the ceramic surface, while the metal surface is kept constant and equal to the reference temperature T 0 = 300 K.
The nondimensional natural frequencies ω = ω a 2 /h ρ m /E m of the FGM sandwich plates for SSSS, CCCC and SCSC boundary conditions, subjected to different temperature profiles following from the solution of the thermo-mechanical analysis and for different material gradients defined by the power-law index are tabulated in Tables 7-9, respectively. In Tables 7-9, for the sake of controlling the accuracy of simulations, the first two frequencies have been compared with those presented in [39]. By inspecting Table 7 to Table 9, one can observe that the first two frequencies obtained from the present finite element model involving the 3D graded element are in a good agreement with the finite element results reported in [39] for all the thermal and displacement boundary conditions considered in the calculations. In addition, to show the effect of temperature on the natural frequencies of the SSSS sandwich plates with different material gradients in SUS304/Si 3 Ni 4 cores, several nondimensional frequencies from Table 7 are presented as a two-dimensional function of the temperature and the power-law index in Figure 7. It is clearly seen from the plots that all the frequencies decrease with increasing temperature for each value of p, in other words, as expected, the FGM sandwich plates become more compliant because of the decrease in material stiffness at higher temperatures. Here, the variation of the volume fraction exponent affects the natural frequencies to a greater extent than the temperature in the considered range. It is also important to mention that the natural frequencies of the CCCC and SCSC sandwich plates exhibit similar responses with rising temperature, and by increasing the power-law index. Table 9. The first nine nondimensional frequencies of two edges fully supported and two edges clamped (SCSC) SUS304/Si 3 Ni 4 functionally graded square 1-8-1 sandwich plates thermally loaded by the nonlinear temperature rise.

Conclusions
The free vibrations of FGM sandwich plates under temperature loading conditions are examined. The natural frequencies of thermally loaded FGM plates are computed using a model based on the 3D graded finite element developed within the ABAQUS code environment. The material gradient was assumed to vary in the thickness direction of the plates according to a power-law distribution of the volume fractions. The rule of mixture was used to evaluate the effective material properties of the FGM. The FGM was implemented into the conventional 3-D elements of ABAQUS code via a combination of user-defined subroutines such as UMAT, UMATHT, and USDFLD (the codes can be downloaded from http://polonez.pollub.pl/deliverables/). In the simulations of FGM sandwich plates, the thermomechanical analysis, used to obtain a temperature profile and associated with temperature-induced displacement and stress fields, couples with a frequency analysis to calculate natural frequencies and mode shapes accounting for a temperature-defined base state. The latter analysis adopted the average mass density instead of its spatial distribution adopted by the former one. The convergence analysis of the present FE model has been done to validate the accuracy of the numerical results by comparing them with the solutions previously reported in the literature and to estimate the model computational efficiency. The effects of different thermal loads, imposed at the external surfaces in steady-state conditions, and different displacement boundary conditions and material parameters, associated with a variety of volume fractions of the material constituents, on the frequencies of FGM sandwich plates are discussed in detail. The following conclusions can be drawn from the present study: • The use of the graded finite elements for analyzing FGM sandwich plates provides a more efficient modeling approach than homogeneous elements to achieve high-fidelity results.

•
The solutions of thermomechanical analysis reveal that a temperature profile calculated with the 3D model is important to predict correct thermal-induced displacement and stress distributions, which, in turn, affect the accuracy of calculated natural frequencies and mode shapes of thermally loaded FGM plates. The thermomechanical analysis also permits the mode shapes to be analyzed in terms of the temperature and stresses.

•
It is observed from the simulations that the natural frequencies decrease as the volume fraction of ceramic decreases across the thickness of the FGM plates.

•
The natural frequencies have a tendency to decrease with an increase of temperature for each of the functionally graded material profiles studied.

•
This work forms a convenient tool for subsequent dynamic analyses of FGM sandwich plates under different temperature conditions with accurate finite element solutions provided by the ABAQUS commercial code.
The developed graded finite element can also be adopted to a 3D crack sensitivity analysis associated with nondestructive testing of FGM sandwich panels and welding-adhesive joints as proposed in [66,67]. This will be a subject of our future research.

Conflicts of Interest:
The authors declare no conflict of interest.