1 Introduction

Circular plates are one of the most important structural elements widely used in civil, mechanical, aeronautical, and electronic engineering^[1-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-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.

2 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:

(1a)

(1b)

where ω is the circular frequency of the harmonic vibration. Q is the shearing force per unit length in the circumferential direction defined by

(2)

M_r and M_t are the radial and tangential bending moments per unit length in the circumferential direction, respectively, which are defined by

(3a)

or

(3b)

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

(4)

where h is the thickness of the plate, E, G, and ρ are Young's modulus, the shear modulus, and the mass density, respectively, and κ² is the shear correction factor. In general, all the above geometric/material parameters of the plate depend on the polar radius r.

3 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:

(5)

where 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

(6a)

(6b)

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.,

(7)

(8)

then the coupled governing differential equations (1a) and (1b) can be simply replaced by a single equation (6a) or equivalently

(9)

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:

(10)

(11)

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

(12)

(13)

(14)

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

(15)

Besides, Eq. (13) is further equivalent to

(16)

Substituting Eq. (15) and the first relation of Eq. (16) into Eq. (9) yields

(17)

where

(18)

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.,

(19a)

The shearing force Q can be constructed immediately from Eq. (5) by

(19b)

The radial bending moment M_r and the tangential bending moment M_t are obtained by substituting Eq. (19a) into Eq. (3b), i.e.,

(19c)

The transverse deflection W is obtained by substituting Eq. (19b) into Eq. (1b), i.e.,

(19d)

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:

(20)

(21)

4 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:

(22)

To satisfy Eq. (22), we further assume that both υ and κ² are constants. In this paper, we always set κ²=π ²/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:

(23)

4.1 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

(24)

Substituting Eq. (24) into Eq. (17) yields

(25)

which can be further transformed into the Bessel/modified Bessel equations^[26]. Γ^± in Eq. (25) denotes Γ⁺ or Γ^- in Eq. (18). Let us define from Eq. (23). The general solution for the rotation angle ψ due to bending, without any singularity at r=0, can be expressed as follows:

(26a)

for 0＜ω＜ω_c and

(26b)

for ω >ω_c, where J₁ and I₁ 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.

4.2 {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:

(27)

where p is an arbitrary dimensionless real constant, R is the radius of the circular plate, and D₀, M₀, I₀, and S₀ 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

(28)

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₀e^{p (r/R)2} and ρ₀e^{p (r/R)2}, respectively. Here, E₀ and ρ₀ 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]:

(29)

One solution of the confluent hypergeometric equation is called the Kummer function, which is defined by^[27]

(30)

where k! is the factorial of the non-negative integer k, and (θ)_k is the Pochhammer symbol defined by

(31)

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:

(32a)

where c₁ and c₂ 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:

(32b)

(32c)

(32d)

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

(33a)

for a free edge with Q=M_r=0 at r=R,

(33b)

for a simply supported (S-S) edge with W=M_r=0 at r=R, and

(33c)

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^-r2/(2R2). Table 1 shows the first five dimensionless natural frequencies 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 , where r₀^g and h₀ are the gyration radius and the thickness of the plate at r=0, respectively.

4.3 {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.

4.3.1 On hypergeometric equations

Suppose that

(34)

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₀, M₀, I₀, and S₀ are the particular values of D, M, I, and S at r=0, respectively. Substituting Eq. (34) into Eqs. (12) and (13) yields

(35)

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₀, E₀, and ρ₀ 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]:

(36)

One solution of the hypergeometric equation is called the hypergeometric function, which is defined by^[27]

(37)

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:

(38a)

(38b)

(38c)

(38d)

where c₁ and c₂ are unknown coefficients to be determined from the specific boundary conditions on the edge of the circular plate, and

(39)

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:

(40a)

for a free edge, i.e., Q=M_r=0 at r=R,

(40b)

for an S-S edge, i.e., W=M_r=0 at r=R, and

(40c)

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]

(41)

The characteristic equation for determining the natural frequencies can thus be derived from Eq. (41) as follows:

(42)

where

(43)

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 natural-frequency branches for the Mindlin plates. We denote them as Ω_branch-1^Mindlin and Ω_branch-2^Mindlin, respectively, where

The associated eigen-modes with Ω_branch-1^Mindlin or Ω_branch-2^Mindlin can be easily obtained by setting c₁=0 or c₂=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:

(44)

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:

(45)

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/R=0.05 with , where r₀^g and h₀ 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 Ω_branch-1^Mindlin only for N=0 and 1; (ii) Ω^rotatory undergoes a transition from Ω_branch-1^Mindlin to Ω_branch-2^Mindlin when N increases; (iii) Ω^shear is a good approximation to Ω_branch-1^Mindlin for all N. The above observations are quite similar to those reported in Ref. ^[25] for beam structures.

4.3.2 On transformed hypergeometric equations

Suppose that

(46)

where p, R, D₀, M₀, I₀, and S₀ have the same definitions as in Eq. (34). Substituting Eq. (46) into Eqs. (12) and (13) leads to

(47)

Substituting Eqs. (46) and (47) into Eq. (17) gives

(48)

The general solution of Eq. (48) for ψ without singularity at r=0 can be expressed in terms of the hypergeometric functions^[27] as follows:

(49)

where

(50)

The expressions for Q, M_r, M_t, and W are also derived but will not be listed here.

4.3.3 On reduced confluent Heun equations

In this subsection, we suppose that

(51)

where p, R, D₀, M₀, I₀, and S₀ have the same definitions as in Eq. (34). Substituting Eq. (51) into Eqs. (12) and (13) gives

(52)

Substituting Eqs. (51) and (52) into Eq. (17) yields

(53)

which belongs to the reduced confluent Heun equation as follows^[29-31]:

(54)

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-31]:

(55)

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.,

(57)

where c₁ and c₂ are undetermined coefficients. The shearing force Q can be further obtained from Eq. (19b) by

(58)

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:

(59)

which can be verified by direct substitution. With the aid of Eq. (59), we can obtain from Eqs. (19c) and (19d) the elegant expressions for M_r, M_t, and W as follows:

(60a)

(60b)

and

(60c)

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:

(61)

Let us fix υ=1/3 (so that n=1) and r₀^g/R=0.05. The first six dimensionless natural frequencies 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 can thus be expressed as follows:

(62)

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).

4.3.4 On Heun equations

In this subsection, we deal with a more complex case with

(63)

where R is the radius of the circular plate, and D₀, M₀, I₀, and S₀ are the particular values of D, M, I, and S at r=0, respectively. p and q are two dimensionless nonzero constants satisfying q² ≠ 4p. Let us denote the two distinct roots of p (r/R)²+q (r/R)+1=0 with respect to r as r₁ and r₂, respectively. We further assume that R≤|r₁|＜|r₂|. Substituting Eq. (63) into Eqs. (12) and (13) gives

(64)

Substituting Eqs. (63) and (64) into Eq. (17) yields

(65)

which belongs to the Heun equation as follows^{[29-30, 32]}:

(66)

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]:

(67)

where the coefficients of the above series are determined by

(68)

Based on the above facts, the general solution for ψ without singularity at r=0 can be expressed as follows:

(69)

where c₁ and c₂ are undetermined coefficients, and

(70)

The shearing force Q is thus obtained from Eq. (19b) as follows:

(71)

In order to express M_r, M_t, and W in terms of the local Heun functions, we establish another entirely new differential formula for an arbitrary n as follows:

(72)

which can be proved by direct substitution. Using Eq. (72), we obtain from Eqs. (19c) and (19d) the expressions for M_r, M_t, and W as follows:

(73a)

(73b)

and

(73c)

As an interesting example, we investigate a solid circular plate with the positive radius r₁ 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₁. With given n, p, and q, Γ has infinitely many discrete values, which can make the series H(σ, q; θ₁, θ₂, 3, n+2; x) convergent, where θ₁ and θ₂ are related to Γ. According to Eq. (70), we have

(74)

Let us set υ=1/3 so that n=1. The first five discrete values of τ for p=q=-0.5 are calculated and listed in the order of 7.670 0, 23.866 5, 46.597 2, 75.854 7, and 111.634 8. Once the parameter τ is obtained, the characteristic equation for the determination of can be derived as follows:

(75)

By setting β=0 in Eq. (75), we can derive the dimensionless natural frequency of the Kirchhoff plate as follows:

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ρ ⁿ, 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.

5 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.

Acknowledgements This paper is dedicated to the late Professor Y. H. PAO (1930-2013), who had been the supervisor, collaborator, and mentor of the two authors.