Analytical Solutions Based on Fourier Cosine Series for the Free Vibrations of Functionally Graded Material Rectangular Mindlin Plates

This study aimed to develop series analytical solutions based on the Mindlin plate theory for the free vibrations of functionally graded material (FGM) rectangular plates. The material properties of FGM rectangular plates are assumed to vary along their thickness, and the volume fractions of the plate constituents are defined by a simple power-law function. The series solutions consist of the Fourier cosine series and auxiliary functions of polynomials. The series solutions were established by satisfying governing equations and boundary conditions in the expanded space of the Fourier cosine series. The proposed solutions were validated through comprehensive convergence studies on the first six vibration frequencies of square plates under four combinations of boundary conditions and through comparison of the obtained convergent results with those in the literature. The convergence studies indicated that the solutions obtained for different modes could converge from the upper or lower bounds to the exact values or in an oscillatory manner. The present solutions were further employed to determine the first six vibration frequencies of FGM rectangular plates with various aspect ratios, thickness-to-width ratios, distributions of material properties and combinations of boundary conditions.


Introduction
Functionally graded materials (FGMs) were first produced in the mid-1980s [1]. An FGM is composed of varying mixtures of different materials, such as ceramics and metals. The material properties of FGMs smoothly and continuously vary, in contrast to conventional laminated composite materials. Consequently, FGMs do not comprise stress singularities formed due to discontinuities in the material properties. FGMs can be designed to possess the high heat resistance and corrosion resistance of ceramics as well as the high mechanical strength of metals. Over the previous three decades, FGMs have been extensively explored in various fields including aerospace, energy, electronics, optics, biomedicine, and mechanical engineering.
Plates are employed in a wide range of mechanical and structural system components in civil, mechanical and aeronautical engineering. The behaviors of FGM plates have attracted research attention. Different reviews [2][3][4][5][6] have provided exhaustive summaries of the studies published on the free vibrations and buckling of FGM plates according to various plate theories and the three-dimensional elasticity theory.
Numerous studies have investigated the free vibrations of FGM rectangular plates, and most of these studies have employed various numerical methods. For example, on the basis of the classical plate theory, Abrate [7] and Zhang and Zhou [8] reported that an FGM plate behaves similar to a homogeneous plate if a suitable reference plane is adopted, while Yang and Shen [9] employed a one-dimensional differential quadrature approximation and the Galerkin procedure to determine the frequencies of initially stressed plates. According to the first-order shear deformation plate theory, Zhao et al. [10] applied the element-free kp-Ritz method to analyze the vibrations of square and skew plates under different combinations of boundary conditions. Fu et al. [11] employed the Ritz method with admissible functions consisting of double Fourier cosine and several closed-form auxiliary functions to study the vibrations of orthotropic FGM plates with general boundary restraints. Ferreira et al. [12] investigated the vibrations of square plates employing the first-order and third-order shear deformation plate theories and the collocation method with multiquadric radial basis functions. Huang et al. [13] adopted the third-order shear deformation plate theory and Ritz method for analyzing the vibrations of rectangular plates with and without side cracks. Hong [14] investigated the thermal vibrations of plates via the generalized differential quadrature method and the third-order shear deformation plate theory. Using the higher-order shear deformation plate theories, Qian et al. [15] applied the Petrov-Galerkin meshless method and Roque et al. [16] adopted the multiquadric radial basis function method to find the vibration frequencies of thick plates. Using the Ritz method and three-dimensional elasticity theory, Uymaz and Aydogdu [17] studied the vibrations of plates with various combinations of boundary conditions, and Cui et al. [18] performed vibration analysis of an FGM sandwich rectangular plate resting on an elastic foundation using admissible trigonometric functions. Huang and his coworkers [19,20] proposed a set of admissible functions, which can accurately describe the behaviors of a crack, for examining the vibrations of cracked FGM rectangular plates and also showed the natural frequencies of intact plates. Burlayenko et al. [21] employed the commercial finite element package ABAQUS to analyze the vibrations of thermally loaded FGM sandwich plates.
Only a few studies have been devoted to analytical solutions for the free vibrations of FGM rectangular plates based on various plate theories. The solutions in these studies consider rectangular plates with two opposite edges or four simply supported edges (faces). Using the first-order shear deformation plate theory, Hosseini-Hashemi et al. [22,23] introduced new potential and auxiliary functions to construct exact closed-form solutions for the vibrations of rectangular plates having two opposite edges simply supported with and without considering the in-plane displacement components, respectively, while Ghashochi-Bargh and Razavi [24] proposed analytical solutions for the vibrations of orthotropic FGM rectangular plates without considering the in-plane displacement components. Hosseini-Hashemi et al. [25] extended their studies by using a third-order shear deformation theory. Considering simply supported conditions on four edge surfaces, Matsunaga [26] and Sekkal et al. [27] developed solutions for FGM rectangular plates and sandwich plates, respectively, according to the higher-order shear deformation plate theories. Based on three-dimensional elasticity theory, Vel and Batra [28] used the power series method to construct solutions for the vibrations of FGM rectangular plates, while Reddy and Cheng [29] employed an asymptotic approach along with a transfer matrix. Huo et al. [30] employed the recursive matrix method to develop the solutions for the vibrations of FGM sandwich plates.
According to plate theories, there are 21 distinct combinations of boundary conditions (i.e., free, simply supported and clamped) for rectangular plates. The literature review found that except for the six cases in which two opposite edges are simply supported, no analytical solution for the vibrations of FGM rectangular plates with various combinations of boundary conditions exists. The present study aims to fill a gap in the literature and proposes analytical solutions based on the Mindlin plate theory for the vibrations of FGM rectangular plates with 21 combinations of boundary conditions. The proposed solutions are established using the Fourier cosine series with polynomial supplementary functions, which eliminate the validity requirement for the term-wise differentiation of the Fourier sine series of a function to accurately represent the differential of the function [31]. The validity of the present solutions is confirmed through comprehensive convergence studies for plates with four combinations of boundary conditions and by comparing the obtained vibration frequencies with those published in the literature. The material properties along the thickness of an FGM plate are estimated using the power-law or the Mori-Tanaka scheme. The solutions are further applied to determine the vibration frequencies of Al/Al 2 O 3 FGM plates with nine combinations of boundary conditions. The obtained analytical results can serve as a benchmark for the evaluation of other solutions obtained through various approximate or numerical approaches.

Material Models
Depicted in Figure 1 is a rectangular FGM plate with a length of a, width of b and thickness of h. The material properties of FGM plates are assumed to vary along their thickness (z) according to the power-law or Mori-Tanaka scheme, which are two popular material models used in the literature. The FGM plates under consideration are made of aluminum (Al) and ceramic (zirconia (ZrO 2 ) or alumina (Al 2 O 3 )), the material properties of which are given in Table 1.
Materials 2020, 13, x FOR PEER REVIEW 3 of 20 are estimated using the power-law or the Mori-Tanaka scheme. The solutions are further applied to determine the vibration frequencies of Al/Al2O3 FGM plates with nine combinations of boundary conditions. The obtained analytical results can serve as a benchmark for the evaluation of other solutions obtained through various approximate or numerical approaches.

Material Models
Depicted in Figure 1 is a rectangular FGM plate with a length of a, width of b and thickness of h. The material properties of FGM plates are assumed to vary along their thickness (z) according to the power-law or Mori-Tanaka scheme, which are two popular material models used in the literature. The FGM plates under consideration are made of aluminum (Al) and ceramic (zirconia (ZrO2) or alumina (Al2O3)), the material properties of which are given in Table 1.
along the thickness of an FGM plate are given as follows: ; Pb and Pt denote the material properties at the bottom face (z = −h/2) and top face (z = h/2), respectively; ΔP is the difference between Pb and Pt , and m is the material property gradient index that governs the material variation profile in the thickness direction. Equation (1) indicates that if Pb = Pt or m = 0, P(z) is constant. FGM plates consisting of Al and ZrO2 (or Al2O3) exhibit a constant Poisson's ratio because Al and ZrO2 (or Al2O3) have the same Poisson's ratio. Figure   2 illustrates the distributions of ( )   A power-law distribution of material properties is often assumed for FGMs. The material properties (i.e., Young's modulus (E = E(z)), Poisson's ratio (ν(z)), and mass density (ρ = ρ(z))) along the thickness of an FGM plate are given as follows: where V(z) = ( z h + 1 2 ) m ; P b and P t denote the material properties at the bottom face (z = −h/2) and top face (z = h/2), respectively; ∆P is the difference between P b and P t , and m is the material property gradient index that governs the material variation profile in the thickness direction. Equation (1) indicates that if P b = P t or m = 0, P(z) is constant. FGM plates consisting of Al and ZrO 2 (or Al 2 O 3 ) exhibit a constant Poisson's ratio because Al and ZrO 2 (or Al 2 O 3 ) have the same Poisson's ratio.  The Mori-Tanaka scheme is also frequently used to describe the material properties of FGMs. The effective mass density along the thickness of an FGM plate is given by where subscripts 1 and 2 indicate materials 1 and 2, respectively, and 1 where After the effective moduli K and G are estimated, the effective Young's modulus and Poisson's ratio are obtained using the following equation: In the Mori-Tanaka scheme, the Poisson's ratio is a function of z even when materials 1 and 2 have the same Poisson's ratio. Equation (2) can be converted to the form of Equation (1), so the density distribution based on the Mori-Tanaka scheme is the same as that described by Equation (1). The distributions of E(z) and ρ( ) z along the thickness of Al/Al2O3 plates with m = 0.5, 2, and 5 are illustrated in Figure 2 and denoted by "M-T."

Governing Equations and Boundary Conditions
In the Mindlin plate theory [32], the displacement components of a plate are expressed as follows: The Mori-Tanaka scheme is also frequently used to describe the material properties of FGMs. The effective mass density along the thickness of an FGM plate is given by where subscripts 1 and 2 indicate materials 1 and 2, respectively, and V t 1 and V b 1 are the volume fractions of material 1 on the top and bottom surfaces of the plate, respectively. The effective local bulk modulus K and the shear modulus G are given by 6(K 1 +2G 1 ) . After the effective moduli K and G are estimated, the effective Young's modulus and Poisson's ratio are obtained using the following equation: In the Mori-Tanaka scheme, the Poisson's ratio is a function of z even when materials 1 and 2 have the same Poisson's ratio. Equation (2) can be converted to the form of Equation (1), so the density distribution based on the Mori-Tanaka scheme is the same as that described by Equation (1). The distributions of E(z) and ρ(z) along the thickness of Al/Al 2 O 3 plates with m = 0.5, 2, and 5 are illustrated in Figure 2 and denoted by "M-T".

Governing Equations and Boundary Conditions
In the Mindlin plate theory [32], the displacement components of a plate are expressed as follows: where u, v, and w are the displacement components in the x-, y-, and z-directions, respectively; u 0 , v 0 , and w 0 are the displacements on the mid-plane, and ψ x and ψ y are the rotations of the mid-plane normal in the xand y-directions, respectively. The stress resultants are defined as follows: where the subscript β represents x or y, and σ ij represents the stress components. The equations of motion are given as follows: . where ρ(z)z l dz (l = 0, 1 and 2). The subscript comma denotes the partial derivative with respect to the coordinates defined by the variable after the comma. Substituting linear strain-displacement (Equation (10)) and stress-strain relationships (Equation (11)) into Equation (8) yields the expressions for stress resultants in terms of displacement related components (Equation (12)). where The parameter κ is the transverse shear correction coefficient and is taken as 5/6 in the following analyses.
Each edge of a rectangular plate is simply supported (S), clamped (C) or free (F). For the edge with y = constant, the S, C and F boundary conditions are defined as follows: Similar definitions of boundary conditions are also applicable for the edge with x = constant.

Series Solutions
To establish the Fourier cosine series solutions for the vibrations of plates, let And mn cosα m x cosβ n y + where α m = mπ/a, β n = nπ/b, and ξ l (x) and η l (y) are supplementary functions. Tolstov [31] showed the following theorem on the differentiation of the Fourier series of a function: is expanded as When f(x) is expanded as λ n a n sin λ n x The theorem indicates that the Fourier cosine series can be differentiated term-by-term, while such an operation can be applied to the Fourier sine series only if f (0) = f (L) = 0. To remedy such shortcoming of the sine series, Li [33] proposed to add some supplementary functions into the cosine series in Equation (26) and to determine f (x) and f (iv) (x) via term-by-term differential.
According to Li [33,34], the supplementary functions are typically determined by satisfying the following conditions: If polynomial functions are used, Equation (27) leads to Substituting Equations (18)-(23) and (28)  ln . For instance, when U 0 (a, y) = 0, which is one of the fixed boundary conditions at x = a, the following equation is obtained: In establishing Equation (29), η l (y) is expressed in its Fourier cosine series as follows: where c ln = b 0 η l (y) cosβ n ydy/ b 0 (cosβ n y) 2 dy. Equation (29) includes N + 1 functions of cos β n y, and each coefficient of cos β n y must equal zero in order to satisfy the equation. Consequently, Equation (29) provides N + 1 linear algebraic equations for the coefficients A where P = (B mn ) T (l = 1, 2; m= 0, 1, · · · M; n = 0, 1, · · · N).
In establishing Equation (33), the following functions are expressed in terms of the Fourier cosine series to factor out cos α m x cosβ n y:   17)). Such set of equations can be further expressed in the following matrix form: Equation (32) gives Substituting Equation (36) into Equation (35) yields which forms an eigenvalue problem.

Convergence Studies and Comparisons
The boundary conditions for a rectangular plate at the edges x = 0, y = 0, x = a and y = b are specified by four letters in a respective series. For example, CSFF boundary conditions mean a clamped boundary condition at x = 0, a simply supported boundary condition at y = 0 and a free boundary condition at x = a and y = b. To validate the proposed solutions, comprehensive convergence studies were carried out for the nondimensional vibration frequencies Ω (=ω(a 2 /h) ρ c /E c , where the subscript "c" indicates ceramic material properties) of the first six modes of square plates with h/b = 0.1 and with SSSS, SCSC, CFFF and FFFF boundary conditions. The obtained results were compared with the published results from the literature. In the following, the solution terms M and N in Equations (19)-(23) are set equal, and M (= M + 1) mainly equals 5, 10, 15, 25, 30 and 35 in the convergence studies. An Al/Al 2 O 3 FGM, whose material properties are described by Equation (1) (the power-law model), is mainly considered. Table 2 presents the convergence studies for Al/Al 2 O 3 and Al/ZrO 2 FGM (m = 1) plates with SSSS boundary conditions. Notably, the material properties of the Al/ZrO 2 FGM are determined by the Mori-Tanaka scheme. In the following tables, the mode denoted by "*" is the in-plane displacement-dominated mode. The results determined from simple exact closed-form solutions (see Appendix A) are given, and the superscripts "(m, n)" denote the wave numbers in the x-direction and y-direction, respectively. Some published results based on the Mindlin plate theory are also given in Table 2 for comparison. The published results presented in the table include (1) the results provided by Hosseini-Hashemi et al. [23], who proposed exact analytical solutions for FGM rectangular plates having two simply supported opposite edges; (2) the numerical results of Zhao et al. [10], who obtained solutions by using an element-free kp-Ritz method with shape functions constructed based on the kernel particle concept; (3) the numerical results of Ferreira et al. [12], who used the global collocation method with multi-quadric radial basis functions.  [12]; "*" denotes the in-plane displacement-dominated model; the superscripts "( )" denote the wave numbers in the x-direction and y-direction.
The present results converge from the upper-bounds of solutions as the number of solution terms increases. The results obtained using M = 15 show the agreement of three significant figures with the results determined from the simple exact closed-form solutions given in Appendix A, and the differences are less than 0.1%. It is interesting to observe that although the solutions of Hosseini-Hashemi et al. [23] provided accurate results for out-of-plane displacement-dominated modes, they did not provide the results for a wave number of zero in the x-direction or y-direction, which corresponded to the in-plane displacement-dominated modes. The first six modal shapes are depicted in Figure 3. In this figure, the contours of out-of-plane displacement (W 0 ) (represented by solid lines) and the nodal lines (represented by dashed lines) are depicted for the out-of-plane flexural modes. Moreover, the in-plane modal deformations are displayed for the in-plane displacement-dominated modes. Compared with the results obtained from the exact closed-form solutions, the results of Zhao et al. [10] exhibited a difference of approximately 1.6% and the results of Ferreira et al. [12] exhibited a difference between 0.9% and 3.3%. Consequently, the present results obtained using M = 15 are more accurate than those of Zhao et al. [10] and Ferreira et al. [12].     [23] for at least three significant figures. Again, the solutions of Hosseini-Hashemi et al. [23] neglected the in-plane mode for a wave number of zero in the x-direction.     [35] and Du et al. [36], who investigated the vibrations of homogeneous rectangular plates. Liew et al. [34] determined the frequencies of out-of-plane modes via the conventional Ritz method with polynomial admissible functions, while Du et al. [36] established series solutions for the in-plane vibrations of rectangular plates. The present results for out-of-plane displacement-dominated modes obtained using M ≥ 25 for the FGM (m = 0 and 1) square plates are consistent with the results of Hosseini-Hashemi et al. [23] for at least three significant figures. Again, the solutions of Hosseini-Hashemi et al. [23] neglected the in-plane mode for a wave number of zero in the x-direction. Table 4 presents the convergence of the non-dimensional natural frequencies Ω of cantilevered FGM (m = 0 and 5) square plates with various numbers of solution terms. Interestingly, the results indicate different convergence trends from those observed from Tables 1 and 2. The values of Ω for the homogeneous plate (m = 0) reveal convergence from the upper-bounds of solutions for the first and third modes, convergence from the lower-bounds of solutions for the second, fourth and sixth modes, and oscillatory convergence for the fifth mode. These findings are not applied to the results for the FGM plate with m = 5. For example, oscillatory convergence is observed for the first and fifth modes. The present results obtained using M ≥ 25 exhibit excellent agreement with those of Liew et al. [35] with the differences less than 0.1%. The differences between the present convergent results and those of Zhao et al. [10] can be larger than 1%.  Similar to Table 4, Table 5 considers the plates with FFFF boundary conditions. Notably, using supplementary functions given in Equation (28) yields singular B p in Equation (32), and its inverse cannot be found for Equation (36). To overcome such numerical difficulties, in addition to the conditions presented in Equation (27), the following conditions are proposed to establish the polynomial supplementary functions: Satisfying Equations (27) and (38) yields (39) Table 5 lists the results obtained using M = 5, 15, 25, 35, 40 and 45. Notably, six rigid body modes with zero frequencies are not considered in the table. The convergence of the numerical results is slower than that of the results presented in Tables 2-4. When the results of the homogenous plate (m = 0) are under consideration, oscillatory convergence is found for the third to sixth modes, while convergence from lower-bounds and upper-bounds is observed for the first and second modes, respectively. The results obtained using m ≥ 35 are consistent with those of Liew et al. [35] with the differences less than 0.1%. The convergence behaviors of the results for the FGM plate with m =5 and FFFF boundary conditions are different from those for the homogeneous plate. The convergence of natural frequencies for different modes is monotonic from the upper-or lower-bounds. The results obtained using M ≥ 35 are in good agreement with those obtained from the Ritz method based on three-dimensional elasticity [20] with the differences less than 0.6%.
To simply demonstrate the convergence rates of the present solutions for other combinations of boundary conditions, Table 6 shows the average relative differences in the Ω values of the first six modes obtained using M = N = 15 and M = N = 35 for Al/Al 2 O 3 FGM square plates with h/a = 0.1 and m = 5 at 16 combinations of boundary conditions. All the differences are less than 0.1%, which indicates that the solutions derived in this study provide accurate results even when M = N = 15 is used. Differences larger than 0.08% occur in the results for SFSF and SFFF boundary conditions, while the differences are less than 0.03% when CSSF, CSCF, CCSF and CCCF boundary conditions are under consideration.

Numerical Results
After validating the proposed analytical solutions through convergence studies, we employed the solutions to determine the first six nondimensional frequencies (      The following inferences are drawn from Tables 7-10 and Figures 4 and 5: 1. The constraint increases when a free boundary condition changes to a simply supported boundary condition. The constraint further increases in a clamped boundary condition. Higher constraint results in higher plate stiffness and larger natural frequencies. Therefore, Ω CCCC > Ω CSSF > Ω CSFF > Ω SSFF and Ω CCCC > Ω CCCF > Ω CCFF > Ω CSFF > Ω CFFF > Ω FFFF (where the subscripts indicate the boundary conditions) if the first six rigid body modes with zero frequencies are considered for plates with FFFF boundary conditions. 2. The Mori-Tanaka material model provides a larger Young's modulus than the power-law material model does; however, both models yield the same density distribution (Figure 2). The following inferences are drawn from Tables 7-10 and Figures 4 and 5: 1.
The constraint increases when a free boundary condition changes to a simply supported boundary condition. The constraint further increases in a clamped boundary condition. Higher constraint results in higher plate stiffness and larger natural frequencies. Therefore, Ω CCCC > Ω CSSF > Ω CSFF > Ω SSFF and Ω CCCC > Ω CCCF > Ω CCFF > Ω CSFF > Ω CFFF > Ω FFFF (where the subscripts indicate the boundary conditions) if the first six rigid body modes with zero frequencies are considered for plates with FFFF boundary conditions. 2.
The Mori-Tanaka material model provides a larger Young's modulus than the power-law material model does; however, both models yield the same density distribution ( Figure 2). Consequently, FGM plates following the Mori-Tanaka material model have larger natural frequencies than those following the power-law material model.

3.
An increase in m results in a decrease in Ω. Furthermore, as displayed in Figures 4 and 5, the change rate of Ω with m gradually decreases with an increase in m. Notably, an increase in m leads to a decrease in the plate stiffness and mass (Figure 2).

4.
No in-plane displacement-dominated mode exists in the first six modes for thin square plates with h/a = 0.02; however, such a mode may exist for moderately thick plates with h/a = 0.1.

5.
The nondimensional frequencies (Ω) of plates with h/a = 0.1 are less than those of plates with h/a = 0.02 because h/a is involved in the definition of Ω. When converting Ω to ω, one finds that the trend is opposite for ω because the plate rigidity increases with h/a.

Concluding Remarks
In this study, analytical solutions based on the Mindlin plate theory were developed for the vibrations of FGM rectangular plates with various combinations of boundary conditions. The solutions were established using the Fourier cosine series with polynomial supplementary functions. Fourth-order polynomial supplementary functions were adopted in the solutions for FFFF boundary conditions, and second-order polynomial supplementary functions were adopted in the solutions for the other boundary conditions. The present solutions were validated through comprehensive convergence studies as well as comparisons with published results and the exact closed-form solutions for plates with SSSS boundary conditions. When increasing the number of solution terms, the trends of convergence in vibration frequencies (i.e., monotonous convergence from upper-bounds or lower-bounds or convergence in an oscillatory manner) varied according to the vibration modes, distributions of material properties, and boundary conditions. The vibration frequencies of the first six modes obtained using M = N = 15 were in good agreement with those obtained using M = N = 35, and the average differences were less than 0.1% for FGM square plates with h/a = 0.1 and m = 5 under 17 combinations of boundary conditions, excluding FFFF boundary conditions. The average difference was about 0.6% for FFFF boundary conditions.
The present solutions were also applied to determine the vibration frequencies of Al/Al 2 O FGM plates with CCCC, FFFF, CFFF, CFSF, SSFF, CSFF, CSSF, CCFF and CCCF boundary conditions. The effects of the plate thickness, material model (Mori-Tanaka and power-law models), and power-law index, m on the vibration frequencies were investigated. With a fixed m, the Mori-Tanaka model yielded higher plate stiffness and larger vibration frequencies for plates than the power-law model. An increase in m caused a decrease in the vibration frequencies of Al/Al 2 O FGM plates. The tabulated data in this study can be used as a standard to judge the accuracy of numerical methods. The present solutions can be simply modified to determine the buckling loads of an FGM rectangular plate under uniform initial stresses and to perform linear static and dynamic analyses of an FGM rectangular plate.