A method for natural frequency calculation of the functionally graded rectangular plate with general elastic restraints

In this study, a uniﬁed solution method to obtain the natural frequencies of the functionally graded rectangular plate (FGRP) with general elastic restraints by using ultraspherical polynomials and the Ritz method is presented. It is supposed that the material distribution properties of the FGRP depend on four different parameters. Based on the ﬁrst order shear deformation plate theory, the displacement ﬁelds of the FGRP are set and boundary and continuity conditions are modeled by adopting the virtual spring technique. All displacement functions of the FGRP including boundary and continuity conditions are expressed by ultraspherical polynomials, and the natural frequencies of the FGRP are obtained by applying the Ritz method. The accuracy of the proposed method is veriﬁed by comparison with the results of the previous literature and ﬁnite element method. Finally, the effects of several parameters on the natural frequencies of the FGRP with various classical and elastic restraints are also discussed, which can be helpful to the further development of the same topic.


I. INTRODUCTION
It is an important problem to analyze accurately and quickly the vibration characteristics of plates, which are widely being used in various engineering applications. In particular, with the development of the science and manufacturing technology, various composite plates that can guarantee the required strength and stability are widely used, so it is more important to analyze the vibration characteristics of these composite plates. Since the functionally graded material (FGM), which is one of the composite materials, has high heat resistance and can be used in structures under the mechanical load, it is widely being applied in actual engineering applications.
Research on the vibration analysis of FGM plates has been developed by many researchers. Various plate theories such as classical plate theory (CPT), first-order shear deformation plate theory (FSDPT), higher-order deformation plate theory (HSDPT), exponential shear deformation plate theory (ESDPT), hyperbolic shear deformation plate theory, and trigonometric shear deformation plate theory have been reported in the study of free vibration of FGM plates.
Kumar et al. 1,2 to investigate the free vibration behavior of the thin functionally graded rectangular plate (FGRP) used the dynamic stiffness method (DSM), solved the dynamic stiffness matrix by the Wittrick-Williams algorithm, and obtained the natural frequencies.
In this study, they set the displacement field of the FGRP by using CPT. Jin et al. 3 investigated the free vibration of arbitrarily thick functionally graded rectangular plates by using the first-order shear deformation theory. Based on Mindlin's shear deformation plate theory, Hadi and Razavi 4 studied the vibration of orthotropic and functionally graded rectangular plates with clamped and simply supported boundary conditions by selecting the displacement function according to Cramer's rule. Matsunaga 5 set fundamental dynamic equations by two-dimensional (2D) higher-order theory. Sayyad and Ghugal 6 applied exponential shear deformation theory to determine the natural frequencies and critical buckling loads of orthotropic plates. Meiche et al. 7 and Mahi et al. 8 developed a new hyperbolic shear deformation theory and applied it to the bending and free vibration analysis of isotropic, functionally graded, sandwiched, and laminated composite plates. Ferreira et al. 9 proposed a trigonometric shear deformation theory for modeling symmetric composite plates for the first time. In addition, studies on the vibration analysis of FGM plates by using various plate theories can also be found in many other literature. [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26] Various analysis and calculation methods were introduced by researchers for the vibration analyses of plates based on different plate theories. Eftekhari and Jafari 27 solved three-dimensional elasticity solutions by using the differential quadrature method, and Malik and Bert 28 proposed mixed finite element differential quadrature formulation to study the free vibration of the rectangular plate. The finite element method (FEM) is one of the calculation methods to obtain the natural frequencies of composites and FG plates with various boundary conditions, 16,29 and the Ritz method is also widely used. [30][31][32] In addition, a lot of calculation methods such as the Galerkin method, 21,33 meshless method, 18,34,35 Lèvy method, 36 and Navier method 37 have been applied for the vibration analysis of the rectangular plate.
From the review of the previous literature, it can be seen that the FSDT theory is mostly used for the vibration analysis of moderately thick rectangular plates and the Ritz method is suitable for ensuring the convergence and accuracy of the solution. In order to analyze the vibration characteristics of plates by the Ritz method, orthogonal polynomials such as the Fourier polynomial and the Jacobi polynomial have been used. Both the Fourier polynomial and the Jacobi polynomial allow for the flexible selection of displacement functions. However, the Fourier method has a defect that the calculation process is slightly complicated for the introduction of an auxiliary polynomial to remove discontinuities at boundaries. Since the Jacobi polynomial is determined by two parameters, the number of parameters can be reduced and the computational efficiency can be increased by using the ultraspherical polynomial, a special case of the Jacobi polynomial.
In addition, it can be seen that most researchers conducted vibration analysis under the assumption that the material properties of the FGM plate are determined by two classical parameters.
The purpose of this study is to establish a unified solution method for the vibration analysis of the FGRP whose material distribution characteristics are determined by four parameters by using the ultraspherical polynomial. The displacement field of the FGRP is set by the FSDPT, and all displacement functions of the FGRP are expressed by using ultraspherical polynomials. The free vibration analysis results of the FGRP by the proposed method are presented according to the material properties, geometrical dimensions, and boundary conditions.

II. THEORETICAL FUNDAMENTAL FOR ANALYSIS
A. Geometric description Figure 1 shows the geometric model and coordinate system for the vibration analysis of a FGM plate with elastic restraints. As shown in Fig. 1, the length of the plate in the x-axis direction is Lx, the length of the plate in the y-axis direction is Ly, and the thickness of the plate is h in which x, y, and z are set on the middle surface of the plate, respectively. The displacement components in the x, y, and z directions are indicated by the symbols u, v, and w, respectively.
In order to simulate the arbitrary boundary condition, the virtual spring technique is introduced, the stiffnesses of the virtual spring are denoted by symbols ku, kv, kw (liner spring) and kφ, k θ (rotation spring). Therefore, it is possible to simulate the FGM plate under the arbitrary boundary condition while changing the stiffness value from 0 to infinity. In order to increase the accuracy of results, in this work, the domain decomposition method is introduced 38 in which the plate is divided into N segments along the x and y axis directions, and the divided individual segments are connected by using connective elastic springs in the directions of displacement. Figure 2 shows the connection form of the individual segments.

B. Material properties of the FGRP
Generally, the functionally graded materials are made of a mixture of a ceramic and metal (see Fig. 1), and the physical-mechanical properties of the FGM plate can be expressed as follows: 39,40 In Eq. (1), E, ρ, and μ are Young's modulus, density, and Poisson's ratio, respectively. The subscripts m and c denote the metal and ceramic, respectively. As shown in Eq. (1), the material properties of the FGM plate change gradiently in the thickness direction z along the volume fractions Vc in which the volume fractions are defined by four distribution parameters a, b, c, and p, 39,40 Here, the parameter p is called the power-law exponent. By setting parameters a = 1 and b = 0 in Eq. (2), the classical volume fractions are obtained and the corresponding volume fraction is shown in Fig. 3.
As can be seen from Fig. 3(a), for the first distribution FGM I (a = 1/b = 0/c = 0/p), the material composition is continuously varied such that the bottom surface (z/h = −0.5) of the plate is ceramic rich, whereas the top surface (z/h = 0.5) is metal rich. Figure 3

C. Energy function of the FGRP
According to the FSDT, the displacement field at any point of the χth segment in the FGRP is defined as follows: 26,[40][41][42][43]  Based on the small deformation assumption and the linear strain-displacement relationship, the strain component of the FGRP can be expressed as follows: 3,41-43 where ε Based on the generalized Hooke's law, the stress-strain relationships at any point of the FGRP are defined as follows: 40,42,43 ⎧ where σ χ xx , σ χ yy , τ χ xy , τ χ yz , τ χ xz are the stress components and Qij(z) (i, j = 1, 2, 4, 5, 6) are the material stiffness coefficients characterizing the relationship between stress and strain, The constitutive equations representing the relation between the force components and the strain components are defined as follows: [40][41][42][43] ⎧ The stiffness coefficients {Aij, Bij, Cij}(i, j = 1, 2, 6) representing the relationship between the force components and the strain components can be obtained by integrating the material stiffness coefficients with respect to the thickness of the FGRP, According to the linear elasticity theory, the structural potential energy of the χth segment in the FGRP is expressed as follows: 40,41 (10) Substituting Eqs. (5) and (8) into Eq. (10), the detailed expression of the energy function can be obtained. For the simplicity, the structure potential energy function is divided into three parts as follows: where U T , UB, and U TB represent the tensile, bending, and tensilebending coupling potential energy expressions, respectively, 40,41 The kinetic energy of the ith segment in the FGRP is defined as follows: Substituting Eq. (3) into Eq. (15), the kinetic energy function can be rewritten as follows: where inertia items I 0 , I 1 , and I 2 are expressed as The energy function by the external forces can be expressed as where f u, f v, f w and mφ, m θ are the external forces in the x, y, and z direction and external moments. As mentioned in Sec. II A, the boundary and continuity conditions of the FGRP are modeled by the virtual spring technique; therefore, the boundary and continuity conditions depend on the stiffness values of the virtual spring.
The potential energy function stored in the boundary spring can be expressed as 3,41 where k l,e0 (l = u, v, w, φ, θ, e = x, y) and k l,e1 represent the stiffness values of the boundary spring at four edges of the FGRP, respectively. Similarly, the potential energy stored in the connective springs connecting two adjacent segments of the FGRP can be ARTICLE scitation.org/journal/adv expressed as 3,41 where α, α + 1 and β, β + 1 represents the two adjacent segments in the x and y directions, respectively.

D. Unified solution procedure
The convergence and accuracy of the analysis results depend on how to select the displacement. The purpose of this study is to analyze the vibration characteristics of the FGRP by selecting a suitable allowable displacement function. To this end, the authors selected all displacement functions including boundary and continuous conditions as ultraspherical polynomials.
An ultraspherical polynomial is a special case of the Jacobi orthogonal polynomial, which ensures good convergence and high accuracy of calculation. [44][45][46][47] The ultraspherical polynomials P The orthogonality condition is The ultraspherical polynomials P (λ) m (ξ) may also be generated by using the recurrence relation, 48,49 (m + 2λ)P (λ) where m = 1, 2, 3, . . .. Thus, the allowable displacement function of the FGRP is more generalized by the ultraspherical polynomials and can be expressed as where Um, Vm, Wm, Φm, and Θm are the unknown coefficients of the ultraspherical polynomials and M, N are the maxima m-order and norder, respectively. P The total Lagrangian energy functions of the FGRP can be expressed as By using the Rayleigh-Ritz method, the total Lagrangian energy function can be minimized with respect to the unknown coefficients,
Thus, the vibration governing equation of the FGRP can be written as The stiffness matrix K, mass matrix M, and unknown coefficient matrix A are expressed by the following equations: If the right term F in Eq. (27) is zero, the natural frequencies of the FGRP can be obtained.
As seen in the above solving process, unlike the other methods, for the proposed method, all displacement components including boundary conditions and continuity conditions are included in a ARTICLE scitation.org/journal/adv polynomial. Therefore, the energy equation of the whole system can be expressed by one polynomial. In other words, the solution process can be seen as a unified process.

A. Convergency study
In this section, for the consideration of the free vibration characteristics of the FGRP with the proposed method, the convergency of numerical problems is investigated. For the convergency study, the geometrical dimensions and parameters of the FGRP are selected as follows: In the numerical calculation, the terms of polynomials must be truncated at an appropriate finite number by considering the computational efficiency and accuracy of the solution. Therefore, it is necessary to study the convergency for the determination of the number of polynomial terms, Mx and My. Figure 5 shows the convergence characteristics of the dimensionless frequencies of the FGRP with fully clamped boundary conditions when Mx and My increase. As can be clearly seen in Fig. 5, as Mx and My increase, the natural frequency of the plate converges to a constant value, and above a certain value, the natural frequency hardly changes. With the increase in Mx and My, the accuracy of the solution improves, but the computational efficiency decreases. Therefore, in consideration of the efficiency and accuracy of the calculation, in this study, the Mx and My for the vibration analysis of the FGRP will be uniformly selected as Mx = 6 and My = 6.
In order to study the convergency of the frequencies with respect to the number of divided segments, the variation of the frequency was investigated with the increase in the number of segments in the y direction. In this case, the number of segments in the x direction is Nx = 2. In the same way, convergence characteristics were investigated for the x-direction. The variation in frequencies following the increase in the segments is presented in Table I.
From Table I, it can be seen that the frequencies of the FGRP do not change at all when the number of segments reaches a certain value. Considering the efficiency and accuracy of calculation, in this paper, the number of segments is set as Nx = 4 and Ny = 4.
In this study, a virtual spring technique is introduced to model the boundary conditions arbitrarily. However, in the case of classical boundary conditions, a problem arises in which the stiffness value of the virtual spring should be determined. Therefore, the research is carried out to determine the stiffness value of the virtual boundary spring. As mentioned above, there are three linear spring groups  (ku, kv, kw) and two rotating spring groups (kφ, k θ ) on the four edges of the plate. In order to determine the stiffness value of the individual spring, the stiffness value of the spring to be considered is changed and the stiffness value of the other spring is set to infinity to investigate the change in frequency. Figure 6 shows the frequency change of the FGRP according to the change in the stiffness values of individual springs. As shown in Fig. 6, when the stiffness value of the boundary spring is less than 10 4 or greater than 10 12 , the frequencies of the plate are kept constant and the frequencies are increased rapidly in the interval [10 5 , 10 10 ]. Therefore, in order to set the classical boundary condition by the spring stiffness value, it is possible to set the free boundary condition when k l,e = 0 and the clamped boundary condition when k l,e is infinity. However, when k l,e is infinite, the calculation efficiency is low, so in this paper, the stiffness value for the clamped boundary condition is set to 10 14 . In addition, the k l,e values for determining the elastic boundary condition can be set at ranges between 10 5 and 10 10 .
The spring stiffness value corresponding to various boundary conditions of the FGRP is shown in Table II. For example, the boundary condition CSFE 1 means the overall boundary condition of the plate in which we have the clamped boundary condition at x = 0, the simply supported boundary condition at x = Lx, the free boundary condition at y = 0, and the elastic boundary condition at y = Ly.
Next, the convergence studies to determine the stiffness value of the connective spring have been performed in the same way as in the boundary spring, and the results are shown in Fig. 7.
As shown in this figure, if the stiffness value of the connecting spring is more than 10 12 , it can be considered that the individual segments are firmly connected. Therefore, in this study, the stiffness value of the spring connecting the individual segments is set to 10 12 .
Finally, the studies to determine the parameter λ of the ultraspherical polynomial are carried out.
The percentage error (Ω λ − Ω λ=0 )/Ω λ=0 of the dimensionless frequency of the FGRP for parameter λ of the ultraspherical polynomials is shown in Fig. 8. From Fig. 8, it can be seen that the change in the parameter λ of the ultraspherical polynomials does not affect the frequencies of the FGRP and the relative error does not exceed 3 × 10 −4 . Therefore, the parameter λ of the ultraspherical polynomials is set as λ = 0 in this study.

B. Numerical verifications
In this section, the accuracy of the proposed method is verified with the comparison of the results of the previous literature and the   results of the FEM. From Eq. (2), when the power-law exponent is zero or equal to infinity, the FGRP will be considered as an isotropic plate. 39 Table III comparatively shows the results of the prior literature (Fourier-Ritz method) and the dimensionless frequencies of the isotropic rectangular plate under several classical boundary conditions obtained using the proposed method.
In addition, the comparison results of the dimensionless frequency of the FGRP according to the several power-law exponents p are shown in Table IV. In this comparison study, the results of the presented method are compared with those of the local Kriging meshless method.
Through the comparative studies of Tables III and IV, the results of the natural frequencies of the FGRP by the proposed method are in good agreement with those of the previous literature, and these comparison results show that the proposed method is very efficient and accurate for free vibration of the FGRP.
To further prove the accuracy of the proposed method, the natural frequency results obtained by using the proposed method are ARTICLE scitation.org/journal/adv   Fig. 9. The geometrical dimensions and other parameters of the FGRP used in this comparison study are Lx = 1 m, Ly = 1 m, h = 0.1 m, for ceramic (Al 2 O 3 ), Ec = 380 Gpa, ρc = 3800 kg/m 3 , and νc = 0.3, and for metal (Al) Em = 70 Gpa, ρm = 2702 kg/m 3 , and νm = 0.3. It is assumed that the properties of the material are in accordance with FGM I (a = 1/b = −0.5/c = 2/p = 1). For the finite element method, the finite element analysis software ABAQUS (in Element, type-S4R, elements-15625) has been used.
In the case of the FEM, as the number of elements gets larger, the calculation time increases, and for the bigger number of elements, the solution does not converge but may diverge. However, the proposed method converges exactly when the number of segments increases. In addition, the FEM has a disadvantage that when the geometrical dimensions change, attribute definition and element division (called a mesh in the finite element analysis software) must be newly performed to construct a new model. Unlike the FEM, the aforementioned process is omitted in the proposed method, so it can be seen that our method is very efficient. In order to help the readers to understand the free vibration of the FGRP, some mode shapes of the FGRP under fully clamped boundary conditions and simply supported boundary conditions are shown in Fig. 10.

C. Parametric study
In this section, the effects of several geometrical dimensions and material parameters on the natural frequency of the FGRP are investigated through numerical examples. First, the effects of boundary conditions on the natural frequencies of the FGRP are investigated. The geometric dimensions of plate and material properties are the same as Fig. 9. Table V shows the dimensionless frequencies of the FGRP with several classical boundary conditions with the change in power-law exponent p. From Table V, it can be clearly seen that the boundary condition has a great influence on the natural frequencies of the FGRP. For the boundary condition considered in this study, natural frequencies of the FGRP with CCCC boundary conditions are the largest and are the smallest in the case of the SFSF boundary condition. In addition, for all boundary conditions, the natural frequencies of the FGRP are decreased as power-law exponent p increases. Table VI shows the dimensionless frequencies of the FGRP with several elastic boundary conditions with the change in powerlaw exponent p. In the case of the classical boundary condition, it can be seen that the variation of the boundary condition has a great influence on the natural frequencies of the FGRP under elastic boundary conditions. However, unlike the case of classical boundary conditions, when the power-law exponent p increases, the natural frequencies for some boundary conditions are increased.
Next, the effects of the geometric dimension on the natural frequencies of the FGRP are investigated. Figure 11 shows the effect of the thickness on the natural frequencies of the FGRP. The geometrical dimensions and material ARTICLE scitation.org/journal/adv  It can be seen from Fig. 11 that when the thickness of the plate increases, the natural frequencies increase distinctly for all boundary conditions except the E 1 E 1 E 1 E 1 boundary condition. In the case of the E 1 E 1 E 1 E 1 boundary condition, the increase in natural frequency with the increase in thickness is unclear.
Finally, the effects of distribution parameters of the material on natural frequencies of the FGRP are investigated. Figure 12 shows the change in the dimensionless frequencies of the FGRP according to the change in the distribution parameters under several boundary conditions. The geometric dimensions of the plate and material properties are the same as Fig. 11. As shown in Fig. 12, the natural frequencies of the FGRP are decreased as the power-law exponent p increases regardless of the boundary conditions and the three distribution parameters a, b, and c.
In particular, it can be clearly seen that the natural frequencies of the FGRP decrease rapidly in the first interval (0 < p < 10) and then slowly decrease outside this interval.

IV. CONCLUSION
In this paper, a unified solution method is presented to study the free vibrations of the FGRP with general elastic restraints in which the material properties of the FGRP are defined by four distribution parameters. The Ritz method in conjunction with a domain decomposition method is employed to derive the governing equations based on the first-order shear deformation plate theory. The displacement functions of the individual segment and boundary condition are expanded by using ultraspherical polynomials. The boundary conditions are generalized by introducing the virtual spring technique. The convergency, accuracy, and efficiency of the proposed method are verified through the numerical example. In addition, a study on parameters such as the boundary condition, geometric dimension, and material properties, affecting the natural frequency, was also performed. When the thickness of the plate increases, the natural frequencies increase distinctly for all boundary conditions. The natural frequencies of the FGRP are decreased as the power-law exponent p increases regardless of the boundary conditions and the three distribution parameters a, b, and c.

ACKNOWLEDGMENTS
The authors would like to express their gratitude to the unknown referees for carefully reading the paper and their helpful comments. The authors also gratefully acknowledge the support from the Pyongyang University of Mechanical Engineering of DPRK.
No potential conflict of interest was reported by the authors.

DATA AVAILABILITY
The data that support the findings of this study are available within the article.