Exact solutions for axisymmetric flexural free vibrations of inhomogeneous circular Mindlin plates with variable thickness

Circular plates with radially varying thickness, stiffness, and density are widely used for the structural optimization in engineering. The axisymmetric flexural free vibration of such plates, governed by coupled differential equations with variable coefficients by use of the Mindlin plate theory, is very difficult to be studied analytically. In this paper, a novel analytical method is proposed to reduce such governing equations for circular plates to a pair of uncoupled and easily solvable differential equations of the Sturm-Liouville type. There are two important parameters in the reduced equations. One describes the radial variations of the translational inertia and flexural rigidity with the consideration of the effect of Poisson’s ratio. The other reflects the comprehensive effect of the rotatory inertia and shear deformation. The Heun-type equations, recently well-known in physics, are introduced here to solve the flexural free vibration of circular plates analytically, and two basic differential formulae for the local Heun-type functions are discovered for the first time, which will be of great value in enriching the theory of Heun-type differential equations.


Introduction
Circular plates are one of the most important structural elements widely used in civil, mechanical, aeronautical, and electronic engineering [1][2][3][4][5] . In practical applications, the thickness, modulus, and mass density of circular plates usually vary continuously or stepwisely, in one spatial dimension or more, to meet the requirement for optimizing the self-weight distribution, improving the local strength, or tailoring the dynamical properties. In recent years, it becomes very popular to use functionally graded structures [6][7][8][9][10][11] to avoid the interface mismatch on mechanical or thermal properties efficiently.
The dynamical modal analysis for the inhomogeneous circular plates with non-uniform thickness has attracted much interest up to now, and many analytical solutions are available based on the classical Kirchhoff assumption. Conway [12] , Conway et al. [13] , and Lenox and Conway [14] obtained some analytical solutions for the Kirchhoff plates with variable thickness when Poisson's ratio took some particular values. Harris [15] derived the closed-form expressions of the natural frequencies and mode shapes for a circular Kirchhoff plate of lenticular thickness with absolutely sharp edge. Prasad et al. [16] studied the axisymmetric flexural vibrations of Kirchhoff plates with linearly varying thickness with the Frobenius method. Elishakoff and Storch [17] and Elishakoff [18] developed a novel method to calculate the required plate stiffness and the corresponding natural frequency of a given candidate mode shape where the density distribution of the circular Kirchhoff plate was simply supported at its boundary. Wang [19] derived the power series solutions for the axisymmetric vibration of circular Kirchhoff plates with power or shifted power variations in the thickness. Duan et al. [20] showed the generalized hypergeometric function solutions for the flexural free vibrations of the annular Kirchhoff plate with its thickness varying monotonically in arbitrary power. Caruntu [21] found a class of non-uniform circular Kirchhoff plates with the Jacobi polynomials as their mode shapes of transverse free vibrations.
However, there are very few analytical solutions available for the freely vibrating circular plates based on the Mindlin plate theory [22][23] . Xiang and Zhang [24] solved analytically the transverse free vibration of the circular Mindlin plates with multiple stepwise thickness variations. In this paper, we extend our recent work on the axially inhomogeneous beam structures with variable cross sections [25] , and propose an exact analytical method to solve the axisymmetric flexural free vibrations of the radially inhomogeneous circular Mindlin plates with variable thickness. The paper is outlined as follows. The basic equations governing the freely vibrating circular Mindlin plates are summarized in Section 2. The exact analytical method is proposed in Section 3 to solve the coupled governing equations analytically. Based on the above analytical method, several kinds of exact analytical solutions are given explicitly in Section 4 and Appendix A. In Section 4, the Heun-type equations governing the axisymmetric flexural free vibration of circular plates are derived, and two novel differential formulae on the local Heun-type functions are established. Concluding remarks are drawn in Section 5. Based on the exact analytical method in Section 3, an approximate analytical method is further developed in Appendix B to solve the low-frequency vibration problem.

Basic equations
We investigate the axisymmetric transverse free vibration of circular plates with radially varying thickness, modulus, and density in the polar coordinates (r, θ). The origin of the polar coordinate system coincides with the center of the circular plate. Based on the Mindlin plate theory [22][23] , the dynamic equilibrium equations governing the axisymmetric and harmonic free vibrations of the circular plate can be expressed in terms of the polar radius r as follows: where ω is the circular frequency of the harmonic vibration. Q is the shearing force per unit length in the circumferential direction defined by M r and M t are the radial and tangential bending moments per unit length in the circumferential direction, respectively, which are defined by In the above equations, W denotes the transverse deflection of the plate, ψ denotes the rotation angle due to bending, ν is Poisson's ratio, and D, M , I, and S are four combined parameters of the plates defined, respectively, by where h is the thickness of the plate, E, G, and ρ are Young's modulus, the shear modulus, and the mass density, respectively, and κ 2 is the shear correction factor. In general, all the above geometric/material parameters of the plate depend on the polar radius r.

Exact analytical method
In this section, we will solve the coupled governing differential equations (1a) and (1b) in an exact and analytical way. As the first step, we introduce an auxiliary function as follows: where λ rS gives the ratio of the rotation angle due to shearing to the rotation angle due to bending.
Substituting Eqs. (2), (3a), and (5) into Eqs. (1a) and (1b) yields We now simplify Eqs. (6a) and (6b) to the extent possible. The basic idea is that, if we choose λ(r) properly to make Eqs. (6a) and (6b) identical, i.e., 1 rD λ rD ν r then the coupled governing differential equations (1a) and (1b) can be simply replaced by a single equation (6a) or equivalently To achieve this, we must find suitable λ(r), D, M , I, and S satisfying Eqs. (7) and (8) simultaneously. Actually, Eqs. (7) and (8) can be further simplified as follows: where C ω is an integral constant in terms of the circular frequency ω. Because the four combined parameters D, M , I, and S defined in Eq. (4) are all independent of ω, Eq. (11) is equivalent to where α is a constant associated with the radially varying parameters D and M and Poisson's ratio ν, β is a constant reflecting the comprehensive effect of the rotatory inertia and shear deformation, and C + ω and C − ω are two distinct roots of the quadratic equation (11) with respect to C ω in view of Eqs. (12) and (13). Correspondingly, substituting Eq. (14) into Eq. (10) gives Besides, Eq. (13) is further equivalent to Substituting Eq. (15) and the first relation of Eq. (16) into Eq. (9) yields where In summary, if D and M satisfy Eq. (12) and I and S satisfy Eq. (13) or Eq. (16), then we can find two desired λ(x) from Eq. (15) so that the coupled governing differential equations (1a) and (1b) can be degenerated into Eq. (17) of the Sturm-Liouville type. Next, we derive the general solutions for ψ, Q, M r , M t , and W under the restrictive conditions in Eqs. (12) and (13). Let us denote the general solution of the first and second equations in Eq. (17) as ψ + and ψ − , respectively. The general solution for the rotation angle ψ due to bending can be constructed by the superposition of ψ + and ψ − , i.e., The shearing force Q can be constructed immediately from Eq. (5) by The radial bending moment M r and the tangential bending moment M t are obtained by substituting Eq. (19a) into Eq. (3b), i.e., The transverse deflection W is obtained by substituting Eq. (19b) into Eq. (1b), i.e., At the end of this section, we consider three degenerate models: (i) the Kirchhoff plate model, in which both the rotatory inertia and the shear deformation are neglected, i.e., I = 0 and 1/S = 0; (ii) the rotatory plate model, in which the rotatory inertia is considered while the shear deformation is neglected, i.e., I = 0 and 1/S = 0; (iii) the shear plate model, in which the shear deformation is considered while the rotatory inertia is neglected, i.e., 1/S = 0 and I = 0. Our exact analytical method is also applicable to the above degenerate models, where the two parameters β and ς are calculated by Eq. (13) and Eq. (16) in the limit sense as follows: for Mindlin plates, 0 for Kirchhoff plates, for rotatory plates,

Exact solutions for solid circular Mindlin plates
In this section, we will derive from Eq. (17) the exact analytical solutions for several kinds of solid circular Mindlin plates satisfying the restrictive conditions in Eqs. (12) and (13). For the sake of simplification, we assume that the four combined parameters D, M , I, and S of the Mindlin plate, as defined in Eq. (4), also satisfy an additional restrictive condition as follows: To satisfy Eq. (22), we further assume that both ν and κ 2 are constants. In this paper, we always set κ 2 = π 2 /12 [22][23] , which indicates that ς > 2.
Note that Γ ± , as defined in Eq. (18), no longer depend on r because of Eq. (22). Γ + is always positive for all nonzero ω, while Γ − becomes positive only when ω exceeds a critical value ω c , which is defined by Γ − = 0 as follows:

Constant {D, M , I, S}
To check the correctness of the exact analytical method in Section 3, we reconsider a homogeneous solid circular Mindlin plate with uniform thickness, the solution of which has already been obtained [3,22] . The four combined parameters in Eq. (4) are constants in this case, and thus the restrictive conditions in Eqs. (12) and (13) are obviously satisfied, leading to Substituting Eq. (24) into Eq. (17) yields which can be further transformed into the Bessel/modified Bessel equations [26] . Γ ± in Eq. (25) denotes Γ + or Γ − in Eq. (18). Let us define ω c = (ς − 1)/( √ ςβ) from Eq. (23). The general solution for the rotation angle ψ due to bending, without any singularity at r = 0, can be expressed as follows: for 0 < ω < ω c and for ω > ω c , where J 1 and I 1 are the first-order Bessel function and the first-order modified Bessel function of the first kind, respectively [26][27] . Equations (26a) and (26b) are fully consistent with those in Refs. [3] and [22], which verifies the correctness of our analysis.

{D, M , I, S} in exponential form
In this subsection, we investigate an interesting case in which the four combined parameters D, M , I, and S of the solid circular plate vary with r in the exponential form as follows: where p is an arbitrary dimensionless real constant, R is the radius of the circular plate, and D 0 , M 0 , I 0 , and S 0 are the particular values of D, M , I, and S at r = 0, respectively. The restrictive conditions in Eqs. (12) and (13) are exactly satisfied by Eq. (27), which leads to It is interesting to notice that, if the Mindlin plate model is considered, then Eq. (27) actually describes a class of circular plates with uniform thickness, whose radially varying Young's modulus E and mass density ρ are E 0 e p(r/R) 2 and ρ 0 e p(r/R) 2 , respectively. Here, E 0 and ρ 0 are the particular values of E and ρ at r = 0, respectively. Substituting Eqs. (27) and (28) into Eq. (17) gives the following confluent hypergeometric equation [27] : One solution of the confluent hypergeometric equation is called the Kummer function, which is defined by [27] F where k! is the factorial of the non-negative integer k, and (θ) k is the Pochhammer symbol defined by Keeping in mind that there exists no singularity at r = 0, from Eqs. (29) and (19a), we can obtain the general solution for the rotation angle ψ due to bending as follows: where c 1 and c 2 are unknown coefficients to be determined from the boundary condition on the edge of the circular plate. Once ψ is obtained, we can derive the shearing force Q, the radial bending moment M r , the tangential bending moment M t , and the transverse deflection W from Eqs. (19b)-(19d) as follows: In combination with the specific boundary conditions, the above general solutions can be further used to derive the so-called characteristic equation for the determination of the eigen-frequencies and associated eigen-modes of the axisymmetrically vibrating plate. The characteristic equation can be expressed as for a free edge with Q = M r = 0 at r = R, for a simply supported (S-S) edge with W = M r = 0 at r = R, and for a clamped edge with W = ψ = 0 at r = R. Let us set p = −1/2 so that the point at r = R is exactly the inflection point of e −r 2 /(2R 2 ) . Table 1 shows the first five dimensionless natural frequencies Ω = ωR 2 M 0 /D 0 with the accuracy to four decimal digits, in which a comparison is made between four kinds of plate models under three different boundary conditions. Here, ν = 0.3, and r g 0 = h 0 / √ 12 = 0.05R, where r g 0 and h 0 are the gyration radius and the thickness of the plate at r = 0, respectively.

{D, M , I, S} in power form
In this subsection, we turn to consider some other interesting cases, in which the four combined parameters D, M , I, and S vary with r in power form, and satisfy the restrictive conditions (12) and (13) exactly.

On hypergeometric equations
Suppose that where n is an arbitrary dimensionless real constant, p is a dimensionless nonzero constant ranging from −1 to +1, R is the radius of the solid circular plate, and D 0 , M 0 , I 0 , and S 0 are the particular values of D, M , I, and S at r = 0, respectively. Substituting Eq. (34) into Eqs. (12) and (13) yields It is interesting to notice that, if the Mindlin plate model is considered, then Eq. (34) actually describes a class of circular plates with the radially varying thickness h, Young's modulus E, and the mass density ρ expressed as follows: where h 0 , E 0 , and ρ 0 are the particular values of h, E, and ρ at r = 0, respectively. Substituting Eqs. (34) and (35) into Eq. (17) gives the following hypergeometric equation [27] : One solution of the hypergeometric equation is called the hypergeometric function, which is defined by [27] F (θ 1 , θ 2 , γ, z) Keeping in mind that there exists no singularity at r = 0, we can obtain from Eqs. (36) and (19a)-(19d) the general solutions for ψ, Q, M r , M t , and W as follows: where c 1 and c 2 are unknown coefficients to be determined from the specific boundary conditions on the edge of the circular plate, and As the first example, we consider a solid circular plate with the radius R and −1 < p < 1. After simplifications, we obtain the characteristic equation as follows: for a free edge, i.e., Q = M r = 0 at r = R, for an S-S edge, i.e., W = M r = 0 at r = R, and for a clamped edge, i.e., W = ψ = 0 at r = R. An important application of Eqs. (40a)-(40c) can be found in Appendix B [28] . As the second example, we consider a solid circular plate (p = −1 and n 0) with an absolutely sharp edge at r = R. In this extreme case, the plate edge cannot sustain any bending moment or shearing force, which means that the hypergeometric functions in Eqs. (38a)-(38d) must degenerate into the Jacobi polynomials in order to guarantee the convergence of ψ, Q, M r , M t , and W at r = R. This gives [27] The characteristic equation for determining the natural frequencies can thus be derived from Eq. (41) as follows: ςξ 2 Ω 4 − (2ξ(2(N + 1)(N + n + 2) + 2ςN (N + n + 3) + ς(1 + ν)(n + 2)) + 1)Ω 2 + 8(N + 1)(N + n + 2)(2N (N + n + 3) + (1 + ν)(n + 2)) = 0 for N = 0, 1, 2, · · · , (42) where For the Mindlin plate model, the characteristic equation (42) with respect to Ω always has two distinct positive roots since n 0, which indicates the existence of two independent naturalfrequency branches for the Mindlin plates. We denote them as Ω Mindlin branch-1 and Ω Mindlin branch-2 , respectively, where The associated eigen-modes with Ω Mindlin branch-1 or Ω Mindlin branch-2 can be easily obtained by setting c 1 = 0 or c 2 = 0 in Eqs. (38a)-(38d), respectively. However, for the Kirchhoff plate, rotatory plate, or shear plate model, substituting Eqs. (20) and (21) into Eqs. (42) and (43) yields only one positive Ω, i.e., a single natural-frequency branch, as follows: Ω Kirchhoff = 2 2(N + 1)(N + n + 2)(2N (N + n + 3) + (1 + ν)(n + 2)), Ω rotatory = 4 (N + 1)(N + n + 2)(2N (N + n + 3) + (1 + ν)(n + 2)) 2 + 8(N + 1)(N + n + 2)(M 0 R 2 ) −1 I 0 , where N = 0, 1, 2, · · · . We point out that the first expression in Eq. (44) for the Kirchhoff plates is exactly the same as the formula derived by Harris [15] . For very large N , Eqs. (42) and (44) can be written asymptotically as follows: That is to say, when N is very large, the first natural-frequency branch of the Mindlin plates approaches the single branch of the shear plates, while the second branch of the Mindlin plates approaches the single branch of the rotatory plates. As a quantitative study, we set in Eqs. (42) and (44) that n = 1, ν = 0.3, and r g 0 /R = 0.05 with r g 0 = h 0 / √ 12, where r g 0 and h 0 are the gyration radius and the thickness of the plate at r = 0, respectively. Figure 1 shows the natural-frequency branches for different plate models. It can be observed from Fig. 1 that (i) Ω Kirchhoff agrees well with Ω Mindlin branch-1 only for N = 0 and 1; (ii) Ω rotatory undergoes a transition from Ω Mindlin branch-1 to Ω Mindlin branch-2 when N increases; (iii) Ω shear is a good approximation to Ω Mindlin branch-1 for all N . The above observations are quite similar to those reported in Ref. [25] for beam structures.  (12) and (13) leads to Substituting Eqs. (46) and (47) into Eq. (17) gives The general solution of Eq. (48) for ψ without singularity at r = 0 can be expressed in terms of the hypergeometric functions [27] as follows: where The expressions for Q, M r , M t , and W are also derived but will not be listed here.

On reduced confluent Heun equations
In this subsection, we suppose that Exact solutions for axisymmetric flexural free vibrations 517 where p, R, D 0 , M 0 , I 0 , and S 0 have the same definitions as in Eq. (34). Substituting Eq. (51) into Eqs. (12) and (13) gives Substituting Eqs. (51) and (52) into Eq. (17) yields which belongs to the reduced confluent Heun equation as follows [29][30][31] : The local solution of Eq. (54) at the regular singularity z = 0 with the characteristic exponent zero, i.e., the so-called local reduced confluent Heun function, is written as follows [29][30][31] : where the coefficients of the above series are determined by the three-term recurrence relation as follows: (56) Based on the above facts, the general solution for ψ without singularity at r = 0 can be expressed in terms of the local reduced confluent Heun functions, i.e., where c 1 and c 2 are undetermined coefficients. The shearing force Q can be further obtained from Eq. (19b) by It will be very attractive to express M r , M t , and W in terms of the local reduced confluent Heun functions. To achieve them, we establish an entirely new differential formula for an arbitrary n as follows: d dz (z 2 (1 − z) n+1 H c (n + 3; θ, 3, n + 2; z)) = 2z(1 − z) n H c (0; θ, 1, n + 1; z), and As the first example, we consider a solid circular plate with the radius R and a free edge. The characteristic equation in this case can be derived as follows: Let us fix ν = 1/3 (so that n = 1) and r g 0 /R = 0.05. The first six dimensionless natural frequencies Ω = ωR 2 M 0 /D 0 for several discrete values of p are shown in Table 2 as the benchmarks. As the second example, we consider a solid circular plate with the radius R and an absolutely sharp edge (p = −1). The edge condition in this extreme case is equivalent to the finiteness of ψ at r = R. For a given n, θ has infinitely many discrete values, which can make the series H c (n + 3; θ, 3, n + 2; 1) convergent. The first five discrete values of the above parameter θ for n = 1 (i.e., ν = 1/3) are calculated and listed in the order of 7.367 9, 19.563 8, 36.724 1, 58.840 5, and 85.906 6. Once θ is known, the characteristic equation for determining Ω = ωR 2 M 0 /D 0 can thus be expressed as follows: As a special but very interesting case, given ν = 1/3, the natural frequencies of a linearly tapered and homogeneous solid circular Kirchhoff plate with an absolutely sharp edge are easily obtained, i.e., Ω = θ, by setting β = 0 in Eq. (62).

On Heun equations
In this subsection, we deal with a more complex case with where R is the radius of the circular plate, and D 0 , M 0 , I 0 , and S 0 are the particular values of D, M , I, and S at r = 0, respectively. p and q are two dimensionless nonzero constants satisfying q 2 = 4p. Let us denote the two distinct roots of p(r/R) 2 + q(r/R) + 1 = 0 with respect to r as r 1 and r 2 , respectively. We further assume that R |r 1 | < |r 2 |. Substituting Eq. (63) into Eqs. (12) and (13) gives Substituting Eqs. (63) and (64) into Eq. (17) yields which belongs to the Heun equation as follows [29][30]32] : The local solution of Eq. (66) at the regular singularity z = 0 with the characteristic exponent zero, i.e., the so-called local Heun function, is written as follows [29][30] : and As an interesting example, we investigate a solid circular plate with the positive radius r 1 and an absolutely sharp edge (p < 0, q < 0, and p + q = −1). The edge condition is equivalent to the finiteness of ψ on the edge r = r 1 . With given n, p, and q, Γ has infinitely many discrete values, which can make the series H(σ, q; θ 1 , θ 2 , 3, n + 2; x) convergent, where θ 1 and θ 2 are related to Γ. According to Eq. (70), we have By setting β = 0 in Eq. (75), we can derive the dimensionless natural frequency of the Kirchhoff plate as follows: Ω = (τ − 2(n + 1))τ p 2 .
As the closure of Section 4, we point out that the four combined parameters D, M , I, and S of the circular Mindlin plates, as described in Eqs. (27), (34), (46), (51), and (63), indicate the power-law relationship between Young's modulus E and the mass density ρ. In fact, many porous/cellular materials reported in Refs. [33]- [36] really obey the power law E = kρ n , where k and n are constants. For the circular plate structures made of the above porous/cellular materials, the radial inhomogeneity for the modulus and density can be achieved by changing the porosity or microstructure parameters continuously along the polar radius r.

Concluding remarks
In this paper, we concentrate on the axisymmetric flexural free vibrations of the radially inhomogeneous circular Mindlin plates with non-uniform thickness, and propose an exact analytical method to reduce the two coupled governing differential equations with variable coefficients into a pair of uncoupled Sturm-Liouville equations under two restrictive conditions on the geometrical/material parameters of the Mindlin plates. The above two restrictive conditions yield two important constants. One describes the radial variations of the translational inertia and flexural rigidity with the consideration of the effect of Poisson's ratio. The other reflects the total effect of the rotatory inertia and shear deformation.
New exact analytical solutions are derived for the first time and expressed in terms of the well-known special functions when the four combined parameters of the circular Mindlin plates vary with the polar radius r in the complicated exponential or power form. The Heun-type equations are introduced to deal with the freely vibrating circular plates exactly and analytically, and meanwhile two basic differential formulae for the Heun-type functions are discovered for the first time, which will be of great value in enriching the theory of Heun-type differential equations.
For several interesting cases, the lower-order dimensionless natural frequencies with the accuracy to four decimal digits are tabulated as the benchmarks. It can be found that the shear plate model is an excellent approximation to the Mindlin plate model below the cutoff frequency.
Based on the exact analytical method in the main text, an approximate analytical method is developed in the appendix under the assumption of low-frequency free vibrations, which allows us to take into account the total effect of the rotatory inertia and shear deformation in the average sense. The validity and accuracy of the approximate analytical method are fully verified by comparing it with available numerical approaches.

Appendix A
In this appendix, we turn to investigate the annular circular plates satisfying the restrictive conditions in Eqs. (12) and (13). Three cases are considered.
Case 1 On Bessel equations Suppose that " n , where R2 is the outer radius of the annular circular plate, and D2, M2, I2, and S2 are particular values of D, M , I, and S at r = R2. Substituting Eq. (A1) into Eqs. (12) and (13) gives Let us define ωc = (ς − 1)/( √ ςβ) according to Eq. (23). As an example, when n is an integer, the general solution for ψ with 0 < ω < ωc can be expressed as follows [26] : where J and Y are the Bessel functions of the first and second kind, respectively, I and K are the modified Bessel functions of the first kind and the second kind, respectively. The explicit expressions for Q, Mr, Mt, and W are also derived, but will not be listed here. For the special case of the Kirchhoff plates with ν = 1/3, the above solution degenerates into that reported by Conway [12] . Case 2 On Eulerian equations Suppose that " n , where n is an arbitrary dimensionless real constant, and R2, D2, M2, I2, and S2 are defined the same as those in Eq. (A1). Substituting Eq. (A5) into Eqs. (12) and (13) gives Substituting Eqs. (A5) and (A6) into Eq. (17) yields Equation (A7) is of the Eulerian type. Thus, we can get the general solution for ψ as follows:  The corresponding expressions for Q, Mr, Mt, and W are also derived, but will not be listed here. For the special case of the Kirchhoff plates, the above solution degenerates into that reported by Conway [12] . Case 3 On hypergeometric equations Suppose that 8 > < > : D = DR(r(2 −r)) n+2 , M = MR(r(2 −r)) n , I = IR(r(2 −r)) n+1 , S = SR(r(2 −r)) n+1 , n = 1 ν − 2,r = r R , where R is a length scale, and DR, MR, IR, and SR are the particular values of D, M , I, and S at r = R, respectively. Substituting Eq. (A9) into Eqs. (12) and (13)  (A11) Equation (A11) belongs to the well-known hypergeometric equation [27] . The explicit expressions for ψ, Q, Mr, Mt, and W are derived, but will not be listed here. For the special case of the Kirchhoff plates, the above analysis degenerates into that reported by Caruntu [21] .
freely vibrating circular plates with α = constant and β = β(r). We point out that the above approach for circular plates is similar to the analysis for the beam structures in Ref. [25].
To check the validity and accuracy of the above approximate analytical method, we now investigate the homogeneous solid circular Mindlin plate with the thickness h = h0(1 + p(r/R) 2 ) (−1 < p < 1), and compare the obtained results with the available numerical results [28] , where h0 is the plate thickness at r = 0, and R is the plate radius. The above thickness variation of the plate indicates that 8 > > < > > : Substituting the second relation of Eq. (B4) into Eq. (B2) gives Using the expression ofβ in Eq. (B5), we calculate the first three dimensionless natural frequencies ω p ρ0R 2 (1 − ν 2 )/E0 by Eqs. (40a)-(40c) for three different kinds of edge conditions of the Mindlin plate. Our final results for p = 0.5 are listed in Table B1 with the numerical results in Ref. [28]. It can be observed from Table B1 that the maximum relative error between them is less than 1.5%, which fully verifies the validity and accuracy of our approximation.