Free vibrations and static analysis of functionally graded sandwich plates with three-dimensional finite elements

A three-dimensional modelling of free vibrations and static response of functionally graded material (FGM) sandwich plates is presented. Natural frequencies and associated mode shapes as well as displacements and stresses are determined by using the finite element method within the ABAQUS™ code. The three-dimensional (3-D) brick graded finite element is programmed and incorporated into the code via the user-defined material subroutine UMAT. The results of modal and static analyses are demonstrated for square metal-ceramic functionally graded simply supported plates with a power-law through-the-thickness variation of the volume fraction of the ceramic constituent. The through-the-thickness distribution of effective material properties at a point are defined based on the Mori-Tanaka scheme. First, exact values of displacements, stresses and natural frequencies available for FGM sandwich plates in the literature are used to verify the performance and estimate the accuracy of the developed 3-D graded finite element. Then, parametric studies are carried out for the frequency analysis by varying the volume fraction profile and value of the ceramic volume fraction.


Introduction
Sandwich panels are promising thin-walled structural elements combining high performance with lightweight. This is due to their excellent strength-toweight and stiffness-to-weight ratios, acoustic and thermal insulation, protection against impacts and a possibility of tailoring properties for optimizing structural responses [1]. Conventional sandwich plates are tri-layer structures consisting of two thin, strength face sheets (skins) and a thick, lightweight core [2]. However, when dissimilar materials are joined, a high mismatch in material and geometrical properties between them leads to an inter-laminar stress concentration at the face sheet-to-core interface and as a result sandwich composites often suffer from premature failure caused by debonding the face sheet from the core. This issue has already become one of major topics extensively studied over the past several decades involving static, dynamic and fracture problems, e.g. recent papers [3][4][5][6][7][8][9][10][11][12]. Another important aspect is that usual sandwich materials have no high service temperature capabilities. These disadvantages of the sandwich structure can be reduced or eliminated by using a functionally graded material (FGM) as constitutive layers.
FGMs are a new type of composite materials, which were initially developed as thermal barrier coatings to minimize thermal and residual stresses as well as to prevent large deflections under high temperature in aerospace structures and fusion reactors [13]. The original idea was to use heat resistant ceramic on a high-temperature side and, then, to gradually reduce its volume fraction to metal on the other side. Nowadays, this FGM concept has received a wide spreading in automotive, optoelectronic, biomedical and sport industries. Combining two or more material phases often with incompatible properties of individual materials, the specific features desirable for certain engineering applications are accomplished. Thus, with spatially varying volume fraction of one of the constituents from the bottom to top face sheets of the sandwich material, the bond strength enhancement, removing interface stress concentration, providing multi-functionality and ability to control deformation, dynamic response, etc. can be achieved [14]. Since the use of FGMs as a branch of advanced composite materials, in particular functionally graded sandwich plates, increases, there is a high demand in predictions of their responses at a design stage to ensure their reliable exploitation. In this respect, many studies have been performed to estimate the thermal-mechanical behaviour of FGM plates including fracture, e.g. [15][16][17][18][19][20] among of the latest. On the other hand, the analysis of free vibrations and static response of FGM plates is of primary importance for their performance assessment.
A considerable amount of research on free vibrations and static analysis studies of FGM plates and shells can be found in the literature. Various analytical (or semi-analytical) and numerical methods have been developed for this, as recently reviewed in [21]. The main idea of many of them is to incorporate some new methodologies into known methods with appropriate minimal modifications, if needed. In most of analytical approaches, through-the-thickness assumptions of two-dimensional (2-D) plate/shell theories such as classic (CPT), first order shear deformation (FSDT) and high order shear deformation (HSDT) theories in conjunction with a given variation of elastic constants along the thickness direction are used to obtain the solution, e.g. [22][23][24][25][26][27] among of the many others. As concluded in [28], the material gradient along the thickness direction leads to complex deformations of the core that require the use at least a HSDT theory. An improvement of the 2-D FGM theories so-called quasi-three dimensional (3-D) theories accounting for the effect of both shear and normal deformations in the thickness direction to analyse FGM sandwich plates have been proposed in [29][30][31][32][33]. Analytical 3-D solutions for FGM sandwich plates have also been developed, e.g. in [34][35][36][37][38][39]. Although such approaches are sophisticated and require much efforts, they accurately describe complicating effects, do not call for additional assumptions and are usually used as a benchmarks for 2-D solutions and numerical results. Besides, the complexity of analytic solutions for FGM plates which is associated with their geometry different from rectangular, can be successfully overcome by using semi-analytic approaches. For instance, a numerical-analytical approach based on R-functions theory and the Ritz method has been worked out for studying geometrically nonlinear vibrations of functionally graded shallow shells of complex planform in [40]. An approach implementing the differential transform method for free vibration analysis of nonuniform cross-section of functionally graded beams has been presented, e.g. in [41,42].
Nevertheless, the analysis of more general FGM plate shapes as well as loading and boundary conditions requires the use of numerical tools such as the finite element method (FEM). Conventional finite elements are typically formulated assuming constant elastic properties over the whole element. The modelling of FGM structures with such elements is related to the presentation of a graded region by a number of homogeneous strips. This approach requires a fine enough mesh discretization to accurately capture the gradient in materials properties as well as the gradients in calculated strain and stress fields. As a result, layered models are quite cumbersome in preparation of the input data and lead to excessive computational costs, especially in the case of 3-D modelling. Another approach adopting conventional finite elements for modelling FGM plates is based on the assumption that a FGM plate behaves like a homogenous with corrected location of the neutral axis, e.g. [43]. However, this technique is suitable only for 2-D modelling. Alternatively, to improve the efficiency of the FEM for modelling FGM structures, finite elements that include the gradient in material properties at the element level so-called graded elements have been proposed in literature. There are two approaches to model the material gradient within the finite element. In the first one, element material properties are approximated by using the same interpolation functions used for the displacement field as done e.g. in [44,45], whereas in the second one, the material parameters are sampled directly at the integration points of the element, thereby the inhomogeneity is distributed, e.g. [46].
In recent years, commercially available finite element code ABAQUS TM [47] has become popular among engineers and researchers due to its versatility and high accuracy. However, the package presents a difficulty in implementation of a spatial variation of material properties into FE models. To solve this issue, some existing works suggest to use the ABAQUS' UEL subroutine implementing a material gradation in FGM plates via a user-developed finite element, e.g. [48]. While this approach gives a high flexibility in modelling, it requires the knowledge of an experienced user and extensive benchmarks of the element performance before simulations. Alternatively, to assign the gradation of properties in ABAQUS, a user-defined material subroutine UMAT can be used. This option allows exploiting conventional finite elements from ABAQUS' library that do not need to be tested, but which are endowed with necessary properties at the Gauss points. Such approach for implementation of varying material properties has been reported for ABAQUS' 2-D plane strain elements in [49][50][51], where the strain-stress state of FGM pavement and the thermo-mechanical behaviour of FGM plate have been analysed, respectively. Also, this modelling technique has been extended on 2-D plate/shell and 3-D models of FGM plates in [52,53], respectively, however only the static bending analysis has been simulated there.
Despite a substantial volume of works and developed finite elements for the modelling of FGM structures, the literature search reveals that 3-D finite element models are still high required for predictions of dynamic and static responses of FGM sandwich plates. This is primarily motivated by the fact that 3-D solutions are direct outcomes of the analysis and no any assumptions on through-the-thickness deformations like in 2-D models are needed. Thus, the development of 3-D finite elements for high-fidelity dynamic and static predictions of FGM sandwich plates is of importance. In the present work, we propose a simple formulation of the 3-D brick graded finite element within the ABAQUS code for its application to modal dynamic and static analyses of FGM sandwich plates. For this purpose, the userdefined material subroutine UMAT is programmed to assign a smooth variation of elastic properties in a conventional 3-D brick finite element and an averaging procedure is proposed for obtaining the mass density equivalent to its actual gradation profile within the element. This modelling technique is simpler than the UEL option and is efficient compared with analytical solutions because the FE scheme can be applied to any complicated configurations and boundary conditions in sandwich plates. The performance of the 3-D graded element has been demonstrated by the 3-D modelling of free vibrations and static response of simply supported rectangular metal-ceramic FGM sandwich plates. The material is assumed to be linearly isotropic with mechanical properties varying across the plate thickness in accordance with a powerlaw. The accuracy of the graded element has been validated by comparisons with the results available in the literature for FGM sandwich plates. Parametric studies have also been carried out to find out the effect of varying volume fraction profiles and ceramic volume fractions on natural frequencies and associated mode shapes.

Problem statement
For the sake of completeness, the standard mechanical initial boundary value problem is briefly summarized in the section. Throughout this section we adopt the notations usual in most of the books on Continuum Mechanics, to which we address the reader for more details.

Equations of motion
Let us consider a FGM sandwich plate as a 3-D deformable medium occupying the domain X 2 ½0; a Â ½0; b Â ½À h 2 ; þ h 2 bonded by the surface oX & X at an instant of time t 2 ½0; T. The plate is defined 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. Also, the planes z ¼ AEh=2 refer to the bottom oX À and the top oX þ plate surfaces, respectively, where Fig. 1a.
The motion of the plate is described by a displacement field u at time t 2 ½0; T. Assume that during the motion the plate experiences only infinitesimal deformations with a strain measure defined by a tensor e ¼ 1 2 ru þ ðruÞ T À Á , and the stress state associated with the deformations is described by the Cauchy stress tensor r. Also, the plate is subjected to prescribed displacements u on the boundary oX u and prescribed surface traction t on the boundary oX t , the boundaries are such that oX t [ oX u ¼ oX and oX t \ oX u ¼ ø. Neglecting the body forces, the deformations of the sandwich plate during its motion have to obey the following system of equations [54]: where qðxÞ is the density of material, and the superscript dot means a time derivative, i.e. _ u and € u stand for a velocity and an acceleration, respectively, n is an outward unit normal to the boundary surface oX.
Using the dynamic principle of virtual work, the system of Eq. (1) can be rewritten in weak form as follows: for all virtual (kinematically admissible) displacement fields du.

Constitutive equations and FGM
Let the plate be made of a functionally graded material, which, without loss of generality, is a mixture of two material phases such as metal and ceramic. The composite material is assumed to be linearly isotropic with smoothly varying mechanical properties in the z-direction only, Fig and pure ceramic on the top side with negligible small thickness in comparison with the thick metallicceramic FGM core. The constitutive equation of the linear isotropic FGM material is defined by the generalized Hooke's law: where the Lamé constants kðxÞ and lðxÞ of the elasticity tensor C are point-wise functions of location.
Since the gradient in properties occurs only along the plate thickness direction, then the material tensor in (3) in the Voigt notation reads as [54]: Further, following the homogenization technique, the two-constituent microscopically inhomogeneous metallic-ceramic composite material is to be replaced by an equivalent macroscopically homogeneous FGM with effective material properties. The effective mechanical parameters of the FGM are usually evaluated based on the volume fraction distribution and the approximate shape of the dispersed phase [34]. It is assumed that the grading of the volume fraction of ceramic phase from the bottom to top plate surfaces is determined by a power-law function as follows: 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.
There exist different approaches to estimate the effective properties of such equivalent FGM. The Mori-Tanaka homogenization method is the most popular among others in literature. In this scheme, the effective Lamé constants are computed from corresponding effective bulk modulus K and shear modulus l as follows: where f m ¼ l m ð9K m þ 8l m Þ=6ðK m þ 2l m Þ; the subscripts 'm' and 'c' stand for the metallic and ceramic phases, respectively, and V m þ V c ¼ 1. From (5) it follows that the FGM 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 for various values of the power-law index p. Finally, the effective mass density at a point of the FGM is calculated using the 'rule of mixture' in the form:

Finite element solution procedure
A displacement-based FEM framework is used for solving the problem formulated above. Accordingly to this method, the actual continuous model of the sandwich plate is idealized by an assemblage of arbitrary non-overlapping finite elements X ¼ S N e¼1 X ðeÞ interconnected at nodal points. In each a base element the components of the displacement field are approximated by suitable interpolation functions such that a vector of the displacement field u at an arbitrary point of the element can be written in the matrix form as follows: Here, the summation over all nodal points of the base element is intended. Also, N ¼ ½N I ðxÞ is a matrix of the shape functions N I for the displacements associated with a certain node I and U ðeÞ is a vector of the nodal unknowns.
With the adopted displacement approximation (8), the strain-displacement relations can be rewritten through the nodal unknowns in the matrix notation as the following: e ¼ BU, where B ¼ ½B I is the matrix of gradients of the shape functions, i.e. B I ¼ N I ;x . Inserting all these matrix relations into the variational equality (2) and taking into account a linear elastic material definition (3), we arrive at the system of semidiscrete equations of motion at the element level as: where M ðeÞ , K ðeÞ and F ðeÞ are the mass matrix, the stiffness matrix and the force vector, defined by the expressions [55]: Thereafter, the use of the assembly operation ðÞ ¼ A N e¼1 ðÞ over all the finite elements leads us to the global system of finite element semi-discrete equations of motion: It should be noted that damping is not included in the dynamic system of Eq. (11). The governing system of Eq. (11) can be employed to study the statics and free vibrations by ignoring the time and excluding the appropriate terms. For the static (or quasi-static) analysis, the problem reduces to the system of algebraic equations: In the case of free vibration analysis, the governing system of equations should be stated in the form of the eigenvalue problem as follows: where x is an undamped circular frequency and / is a vector of mode shape associated with found frequency x.
3 Three-dimensional graded finite element The modal analysis for extraction of natural frequencies and associated mode shapes and the static analysis of FGM sandwich plates are carried out with the ABAQUS/Standard code using three-dimensional models. Conventional 3-D finite elements available in the ABAQUS' finite element library do not enable to model variation of mechanical properties within the element volume. To eliminate this inconvenience, we have developed within the package environment a graded 3-D finite element incorporating gradients of material properties. For this purpose, either eight-node linear C3D8 or twenty-node quadratic C3D20 brick isoparametric finite elements can be used as a base element, Fig. 2. Since the quadratic elements usually yield better results at less expense than their linear counterparts [55], C3D20 elements will be used in the succeeding calculations. Isoparametric interpolation within the element volume is defined in terms of the natural coordinates ðn; g; fÞ each spanning the range from À1 to þ1, as shown in Fig. 2a. The node numbering convention used for these elements with either full or reduced integration is also demonstrated in Fig. 2b, c. The shape functions of both the 8-node and 20-node brick elements can be presented as follows: i ; for i ¼ 9; 11; 13; 15; i ; for i ¼ 10; 12; 14; 16; The spatial variation of the material parameters have been achieved by coding the user-defined material subroutine UMAT. The routine evaluates the increments of strain and stress vectors and the material Jacobian to provide an implementation of the incremental form of a defined mechanical constitutive law into the code solver. At the end, it updates the stress and strain vectors, and other solution dependent variables over the element volume at the Gauss points, [47]. Since the material constants depend on the position within the element, the calculations of the constitutive material relation will assign the material parameters directly at different integration points of the element. Finally, by means of numerical integration of the tangential stiffness matrix of the finite element as: the actual gradient of material properties is completely specified within the element. In this expression, i, j and k are the Gauss points, jJ ijk j is the determinant of the Jacobian matrix and w i are the Gaussian weights. The matrix C given in (15) contains the material constants k and l depending on the position as shown in (4).
To complete the development of the graded finite element, the mass density, which is also a function of position for the FGM plate, should be incorporated into the element formulation. However, there is a complication in terms of formulating the element mass matrix with sampling the density values at the Gauss points similar to (15). The frequency extraction procedure is a linear perturbation analysis in ABAQUS. The perturbations are assumed to be about the base state of a model and ignore any inelastic effects and any parametric dependencies of the model on temperature or other field and solution dependent variables [47]. Moreover, the code does not provide a direct access to the computation of elements of the mass matrix. The latter is formed with the conventional homogeneous elements, where the element volume is simply multiplied by a given constant mass density. Herein, we suggest a simple ''engineering solution'' to account for changes in inertial properties due to the use of FGMs. A constant value of the mass density required by ABAQUS can be obtained by averaging a proper density distribution over the whole volume of the FGM plate, V as follows: Although the density in the form of an integrated value averaged over the volume (16), does not reflect its actual gradual change, this formula can reasonably approximate the inertial property of a FGM sandwich plate, as shown in all the frequency analyses below.

Numerical results
In this section numerical results for studying the static behaviour and free vibrations of two-constituent metallic-ceramic functionally graded sandwich plates are presented. In each the analysis presented below, first, convergence studies with mesh refinement have been conducted. An optimal mesh was selected from the conditions that differences between the results, namely, a displacement of the central point in the case of static analysis and a fundamental frequency in the case of eigenfrequency problem, obtained with current and finer meshes are no more than 5%, but from the other hand, the calculations should be the best compromise between computational cost and solution accuracy.

Static analysis
The accuracy of the 3-D graded finite element programmed via UMAT in the ABAQUS environment is verified, first, performing the static analysis of a FGM sandwich plate. A simply supported Al/SiC square sandwich plate subjected to a normal pressure given by q 0 sin px a sin py b on the top surface was considered. We accepted the parameters of (5) as V À c ¼ 0, V þ c ¼ 1 and p ¼ 2, the plane dimensions and the thickness-to-length ratio as a ¼ b ¼ 100 mm and h a ¼ 0:2, respectively, and q 0 ¼ 1 MPa. The material constants used in the calculations are listed in Table 1. The numerically calculated results are compared with those obtained analytically in [34].
A 3-D finite element model of the FGM sandwich plate considered in the present work is illustrated in Fig. 3a. The meshes, optimal from the standpoint of the convergence analysis mentioned above, for both the developed 3-D FGM elements with six elements through-the-thickness and the conventional homogeneous 3-D elements with ten elements through-thethickness are shown in Fig. 3b, c, respectively. It should be noted that the latter model requires a tedious and cumbersome procedure for preparing and, then, introducing the input data different for each the layer of elements in the mesh, while the former one uses a simpler pre-processing, i.e. the input data are assigned on the whole finite element model at once.
The physical quantities compared with the solutions reported in [34] are presented in a dimensionless form as follows: The results of comparison between the analytic solutions in [34] and the results obtained from the both graded and layered models in the static analysis are summarized in Table 2, where D; % denotes relative errors between the reference data and the results of the model with FGM elements and the layered model with homogeneous elements, relatively. A good agreement with the analytic solutions can be seen for both the finite element models. However, the graded model is superior to the layered one because it is more accurate and efficient, i.e. gives better results with less number of the elements and is more convenient for using. The minor discrepancies between the results of the model with FGM elements and the reference values may be addressed to the problem of non-equivalency between the ways of applying the boundary conditions to plate edges in 3-D analytic and FEM models, as discussed elsewhere in literature, e.g. in [56,57]. The distributions of the dimensionless deflection and longitudinal r xx and transverse normal r zz stresses through-the-thickness of the FGM sandwich plate, computed at a point ð a 2 ; b 2 ; zÞ using both the FGM elements in the continuous model and the homogeneous elements in the layered model, in comparison with some known analytic values given in [34] are shown in Fig. 4. One can see that the numerical models have very close results to the analytic solutions and they are able to accurately reproduce the throughthe-thickness behaviours. Although both the numerical models are slightly stiffer than the analytical one due to the discretized nature of the models, they give the exact result for the transverse normal stress and the results which are very close to the exact solution for the longitudinal stress. In doing so, the layered model is a bit stiffer compared with the graded elements and exhibits a stepwise change of the longitudinal stress (see Fig. 4d). Such a step-type variation with constant longitudinal stress in each element is inherently related to homogeneous elements and cannot be avoided in modelling with layered models. In contrast to this, the graded model provides a continuous stress variation through the plate thickness.
To show the efficiency of the developed graded 3-D finite element for studying FGM sandwich plates, the influence of the gradation profile on the static response of the simply supported square FGM sandwich plate is further analysed. Figure 5   it is clearly seen that the deflections of the FGM sandwich plates are not symmetrical with respect to the plate midplane and this non-symmetry is larger when more aluminum phase is in the graded material. The stresses have also nonlinear patterns across the thickness. These conclusions are in a compliance with the results in [29,36]. With increasing the exponent p the deflections of the bottom, middle and top layers are reduced smoothly and same manner to the values close to deflections of the metallic plate (Fig. 5d), whereas the longitudinal stress on those layers varies in a more complicated way (Fig. 5e). The tensile longitudinal stress r xx at the bottom layer first a bit reduces with increasing the metal phase and, then, its value smoothly rises with growing the index p and, finally, tends to the magnitude of stress in the metallic plate. In contrast to this, the compressive longitudinal stress r xx at the top layer first remarkably increases with increasing the metal phase and, then, it gradually decreases with increasing the exponent p nearly to the stress of the metallic plate. The transverse normal stress r zz at bottom and top layers almost does not depend on the gradation profile, and its behaviour at the midplane is similar to the longitudinal stress at this layer, i.e., first, increases and, then, slowly tends to the stress of the metal plate, as seen in Fig. 5e, f.

Free vibration analysis
The correctness of the developed 3-D graded finite element in conjunction with the proposed procedure of the mass density averaging for finite element modelling of free vibrations has been verified by comparing the fundamental frequency, displacements and stresses of a simply supported square Al=ZrO 2 sandwich plate with the results available in [36]. The material properties of the FGM plate studied are listed in Table 1. The dimensions and the parameters of the power-law given by (5) for this FGM sandwich plate are taken the same as the previous static analysis. Table 3 shows the comparisons between the reference values and the results obtained with the 3-D FGM elements model and the layered model with the conventional homogeneous 3-D elements.
It is known that besides standard flexural and inplane predominant modes resulted from 2-D models of FGM sandwich panels, the 3-D solution predicts additionally thickness-stretching predominant modes  Table 3 The comparison of the results of free vibration analysis between the present calculations and the analytic solution in [36] x  [35]. Thus, to validate the ability of the proposed graded element for producing such response, the deflection and stresses at several points of the thickness have also been considered for comparison with analytical values in [36]. All the compared quantities have been nondimensionalized by the formulae: A glance at the values in Table 3 reveals that the computed results are in satisfactory agreement with those presented in [36]. The finite element solutions of the both models with FGM and homogeneous elements are almost the same, and they slightly overestimate the analytical reference values. The sources of these discrepancies could be, first, the normalization formulae involve a deflection which is not exact value itself, but being calculated. Second, not actual distribution of the mass density, but its average value is implemented into the model. Finally, because of an incompatibility between the ways of applying boundary conditions in the FEM and analytical techniques. It should be noted that for the sake of simplicity, needed to carry out a large number of computations, the simply supported conditions in the present research have been realised using the nodal displacement constraints applied over each a whole edge of the sandwich plate. As expected, the 3-D finite element frequency analysis has been able to reveal the existence of peculiar mode shapes associated with the natural frequencies of the thick sandwich plate. That is apart from the usual bending or flexural mode shapes, one has been observed the following additional modes: the pumping or thickness-stretching mode, where the deflections are predominant and are symmetric with respect to the mid-plane; the in-plane or extensional modes, where one of longitudinal displacements prevails; and finally so-called bi-inplane modes or extensional modes, where two longitudinal displacements act simultaneously and either symmetrically or anti-symmetrically. Such mode shapes have been mentioned, but not exhibited in [35]. These types of the mode shapes and the differences between them for metallic, Al=ZrO 2 sandwich and ceramic plates are illustrated in Fig. 6. One can clearly see that those peculiar mode shapes are different between the two homogeneous plates and they, in turn, differ from the functionally graded plate. The latter has the configurations of the mode shapes which are intermediate between the metallic and ceramic plates. Herewith, the influence of the material gradation is more evident on thickness-stretching and bi-inplane extensional mode shapes than on in-plane extensional ones.
The ability of the developed 3-D graded element to properly capture the gradation profile across the plate thickness for the frequency analysis has also been examined. The fundamental frequency of the simply supported Al=ZrO 2 square thick sandwich plate considered above is calculated for different values of the power index p. The results of calculations are compared with those known in [29]. The values of the frequencies normalized by x ¼ xh ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi q m =E m p are exhibited in Table 4, where the numerical results obtained by using both the FGM elements and the homogeneous elements in the layered model are compared with the analytic solutions.
It is obviously that the numerical results show a good agreement with the analytic solutions. This confirms the accurate performance of the developed 3-D graded finite element. Moreover, it should be mentioned that due to the convergence requirement, the element size of the layered model was smaller than the size of the graded element as 1 mm and 2.5 mm, respectively. As a result, the total CPU time of calculations with homogeneous elements was about 3 times longer than the computations with the graded elements. This demonstrates the efficiency of the 3-D graded element in the predictions. For this FGM sandwich plate, the first sixteen dimensionless natural frequencies depending on the power index p are computed and summarized in Table 5.
For the sake of better illustration of the power-index dependence of the natural frequencies of the Al=ZrO 2 sandwich plate, the evolutions of four frequencies with increasing the metal phase are depicted in Fig. 7. One can see that, first, quantitatively the natural frequencies of the FGM sandwich plate are between those of the same metallic and ceramic plates. Second, the power-index or an assumption about the material phase distribution through-the-thickness of the FGM sandwich plate affects its natural frequencies such that with increasing p the frequencies tend to values close to those of the metallic plate. However, the rate of this trend is different for the lower frequencies, as illustrated in Fig. 7a, b and Fig. 6 In-plane, pumping and bi-inplane mode shapes of square thick: a-c metallic plate; d-f Al=ZrO 2 sandwich plate; and g-i ceramic plate shown in Fig. 7c, d. The former ones are slightly changed at the small p and, then, slowly are reduced to the frequencies of the metallic plate with growing the exponent to infinity. In contrast to this, the latter ones decrease quickly enough with increasing the metal phase and tend to the frequencies of the metal plate.  Fig. 7 The variations of the dimensionless natural frequencies with increasing the power-index p: a mode (1,1); b mode (1,2); c inplane mode no. 9; and d pumping mode no. 14 c Fig. 8 The variations of the dimensionless natural frequencies with increasing the power-index p and the volume fraction V c : a mode (1,1); c mode (2,2); e in-plane mode no. 9; and g pumping mode no. 14 (a) Finally, the influence of both the ceramic volume fraction on the top surface and the metal phase in the sandwich core on the modal dynamics of FGM sandwich plates have been studied. The same Al=ZrO 2 sandwich plate as considered above was modeled. Different values of V þ c such as 0.0, 0.3, 0.5, 0.7 and 1.0 at a zero-value of V À c were used in the simulations. The variations of four natural frequencies depending on the both parameters of V þ c and p are shown in Fig. 8 as 2-D and 3-D plots. As is evident from these plots, the natural frequencies of the simply supported Al=ZrO 2 square thick sandwich plate regardless the position in the frequency spectrum become less sensitive to increasing the metallic phase in the gradation profile with decreasing the volume fraction of ceramic on the top surface. Thus, the concentration of the ceramic constituent is an important factor for the dynamic design of FGM sandwich plates. More detailed studies of the role of different constituents on the modal dynamics of FGM sandwich plates are still needed, but they are out of the scope of this work.

Conclusions
The static response and free vibrations of simply supported metal-ceramic FGM square thick sandwich plates have been analyzed by using the three-dimensional finite element modelling with the ABAQUS code. The effective material properties at a plate point were assumed to be defined in accordance with the Mori-Tanaka approach. A power-law variation of the volume fraction of the ceramic phase through-thethickness of the sandwich plate was adopted. The 3-D brick graded finite element incorporating the variation of material properties in direction of the plate thickness was programmed and, then, was implemented into ABAQUS via the user-defined material subroutine UMAT (the code can be downloaded from http:// polonez.pollub.pl/deliverables/). A simple procedure was proposed to estimate an average constant value of the mass density within the graded finite element. The use of this element in predictions turned out to be accurate and very efficient since the graded element unlike the layered model with conventional homogeneous elements neither requires extensive procedure of input data preparation nor fine enough mesh, but it allowed one using the whole power of the ABAQUS code.
With the 3-D graded finite element developed and implemented into ABAQUS, the calculated displacements and stresses in the static analysis and the computed natural frequencies in the free vibration analysis of the FGM sandwich plates have been found to match well the analytic solutions reported in literature. The minor discrepancies between the results of the finite element models and the reference values could be addressed to the problem of non-equivalency of boundary conditions applied to the plate edges between analytic and numerical 3-D modelling procedures. Also, the predictions showed that the deflections and the values of natural frequencies of the FGM sandwich plates are between those for pure ceramic and pure metallic plates. However, the stresses have a nonlinear pattern across the sandwich plate thickness and their variations with increasing the metallic phase have more sophisticated character, especially for the longitudinal stress. The gradients in the material properties affect the natural frequencies of a simply supported square FGM sandwich plate. The values of the natural frequencies with growing the metallic phase tend to those of a pure metallic plate. However, the lower frequencies are less sensitive to the powerlaw index than the higher ones. It was also found that the mode shapes associated with the frequencies of the thick sandwich plates exhibit not only flexural dominant deformations, but thickness-stretching dominant and extensional dominant deformations as well. Yet, it was found out that the volume fraction of ceramic on the top surface is an important factor for dynamic design of FGM sandwich plates. The predictions demonstrated that the less is the ceramic concentration on the top surface, the less sensitive are the frequencies to increasing the metallic phase in the sandwich core. Finally, it needs to mention that although the present results are demonstrated only for the simply supported metallic-ceramic FGM square sandwich plate, the graded 3-D element developed in this work can be used for the 3-D modelling of sandwich plates with any other geometry, gradation profile and boundary conditions. Thus, the results presented in this paper may provide a benchmark for studying statics and modal dynamics of the FGM sandwich plate by other methods and models.