Insight Into the Separation-of-Variable Methods for the Closed-Form Solutions of Free Vibration of Rectangular Thin Plates

1  Introduction

The free vibration problems of rectangular plates can be classified into three categories according to boundary conditions (BCs), Category 1: all edges are simply supported and/or guided; Category 2: only two opposite edges are simply supported and/or guided; Category 3: BCs not falling into any of the above two categories. This classification is also suitable for the eigenbuckling problems of rectangular plates.

The problems of Categories 1 and 2 have the well-known Navier and Levy types of exact solutions that rigorously satisfy characteristic partial differential equations (PDEs) and BCs [1–5]. For the problems of Category 3, it is hard or even impossible to find exact solutions. To obtain analytical solutions to the problems of Category 3, several analytical methods have been developed from different perspectives since the 1950s, such as the Kantorovich-Krylov (K-K) method [6–10], the separation-of-variable (SOV) methods [11–15], the superposition method [16–19], the series expansion-based methods [20–24], the symplectic eigenfunction expansion method [25–29], and the dynamic stiffness method [30–34] for the plate assemblies.

Since the present work focuses on the SOV methods obtaining closed form eigensolutions, the motivations of the SOV methods are presented below. In the Navier method (an inverse method) and Levy method (a semi-inverse method), the separable mode function has the form of w(x, y) = φ(x)ψ(y), in which both φ(x) and ψ(y), or any one of them is constructed based on the simply supported BCs, limiting the application scope of the inverse methods. To extend the application of the Levy method, the direct SOV method [11] also using w(x, y) = φ(x)ψ(y) is proposed from a new perspective, wherein both unknown φ(x) and ψ(y) are solved simultaneously, and the closed-form solutions for the free vibration of rectangular plates with clamped adjacent edges were obtained for the first time. Afterwards, a few more powerful SOV methods are proposed based on the Rayleigh quotient, including the imSOV method [13], the vSOV method [12], the iSOV method [12], and the eSOV method [14]. Among them, the imSOV, iSOV, and eSOV methods are general ones, implying that they are suitable for rectangular plates and cylindrical shells with arbitrary homogeneous BCs.

In the SOV methods, mode functions have separable forms, like w(x, y) = φ(x)ψ(y). The relationships of the eigenvalues corresponding to spatial and temporal coordinates are achieved according to characteristic PDEs or/and the Rayleigh quotient. Transcendental eigenvalue equations are attained with the homogeneous BCs. Although there are several publications for each SOV method, there are no works about the following things: the general theoretical frameworks, the practical calculation method for natural frequencies and modes, the accuracy of frequencies of different orders, and the problem of missing frequencies or solution integrity.

In this context, this work presents the general methodology of how to solve the SOV solutions of characteristic PDEs, and how to obtain closed-form analytical solutions through the Rayleigh quotient. In addition, bisection-based solution procedures are presented to solve eigenvalue equations, and nodal line patterns which are used to check the accuracy and integrity of the SOV solutions are investigated. Note that the free vibration analysis of an isotropic rectangular thin plate is employed to achieve the objective of the present work.

The rest of this manuscript is organized as follows. Section 2 presents the theoretical framework of the SOV methods. After presenting the designing purposes and nonlinear eigenvalue equations of each SOV method, Section 3 gives the bisection-based solution procedures. Section 4 investigates the relations of nodal line patterns and frequency orders, as well as the accuracy and integrity of solutions. Section 5 compares the results of the SOV methods with those of the finite element method (FEM) and other methods, and conclusions are drawn in Section 6.

2  Theoretical Frameworks of the SOV Methods

The SOV methods are used to solve characteristic PDEs, the stationary value problems of the Rayleigh quotient, or both simultaneously. This section first presents the theoretical framework of the SOV methods, including the introduction of the characteristic PDEs and the Rayleigh quotient, as well as the methodology of finding SOV solutions based on them, and then demonstrates that all SOV methods come from the framework.

Fig. 1 shows a rectangular plate with length 2a, width 2b, and thickness h. According to the Kirchhoff plate theory [35], the displacement functions of thin plates have the forms as

U(x,y,z,t)=−z∂W(x,y,t)∂xV(x,y,z,t)=−z∂W(x,y,t)∂y

W(x,y,z,t)=W(x,y,t)(1)

where x, y are the coordinates of the middle surface, z is the coordinate of the thickness direction, and the origin of the coordinates is at the center of the middle surface; U, V, and W represent the displacements in the x, y and z directions, respectively, and only the deflection W is the independent variable according to Eq. (1).

Figure 1: A rectangular plate and the Cartesian coordinates

2.1 Characteristic Partial Differential Equations and Boundary Conditions

For a harmonic motion, the deflection W(x,y,t)=w(x,y)eiωt, where i=−1, ω is the radial frequency, and w(x,y) is the mode function. The equation of the lateral translational motion for the free vibration of thin plates are [35]

∂Qx∂x+∂Qy∂y=−ρhω2w(2)

where ρ is the density, and h is the thickness; Qx and Qy are the transverse shear forces, which have the forms as

Qx=∂Mx∂x+∂Mxy∂y,Qy=∂Mxy∂x+∂My∂y(3)

Mx=−D(∂2w∂x2+ν∂2w∂y2),My=−D(∂2w∂y2+ν∂2w∂x2),Mxy=−D(1−ν)∂2w∂x∂y(4)

Vx=Qx+∂Mxy∂y,Vy=Qy+∂Mxy∂x(5)

By substituting Eq. (4) into Eq. (3) and then into Eq. (2), the characteristic PDE is obtained as

D(∂4w∂x4+2∂4w∂x2∂y2+∂4w∂y4)−ρhω2w=0(6)

With the coordinate transformations ξ=x/a, η=y/b, Eq. (6) changes to

∂4w∂ξ4+2α2∂4w∂ξ2∂η2+α4∂4w∂η4=Ω4w(7)

where the aspect ratio α=a/b and the non-dimensional frequency Ω4=ω2(ρha4/D). All homogeneous BCs expressed in terms of w are listed in Table 1. Besides, symbolisms are used to denote the BCs of four edges. For example, “CSGF” indicates that the edges of ξ=−1, η=−1, ξ=1 and η=1 are clamped (C), simply supported (S), guided (G) and free (F) respectively, and other combinations of BCs can be interpreted in the same way.

2.2 Separation-of-Variable Solutions

This section introduces the general method of solving Eq. (7). Assuming that the closed-form mode function (deflection) w has the separable form as

w(ξ,η)=ϕ(ξ)ψ(η)(8)

then one can see from Table 1 that S, C, and G are separable, but F is not. Here ‘separable boundary condition’ implies that the boundary condition depends only on ϕ(ξ) or ψ(η). However, there is a special case: if a pair of opposite edges are S-S, G-G, or S-G, then even if the other two edges are free, their BCs are also separable. Section 3 will give an integration method to convert inseparable free BCs into separable ones.

Assuming the functions ϕ(ξ) or ψ(η) in Eq. (8) are as follows:

ϕ(ξ)=Aeμξ,ψ(η)=Beλη(9)

where μ and λ are the spatial eigenvalues concerning spatial coordinates ξ and η, respectively. With the substitution of Eq. (9) into Eq. (7), we have

(μ2+α2λ2)2=Ω4(10)

or

μ2+α2λ2=±Ω2(11)

which is the explicit relations of the two spatial eigenvalues μ and λ and the frequency Ω. One can see that there are four μ and four λ satisfying Eq. (10), and they can be written as

μ1,2≜±iα1,μ3,4≜±β1(12)

λ1,2≜±iα2,λ3,4≜±β2(13)

where α1, β1, α2, β2 are real numbers. Accordingly, the closed-form ϕ and ψ can be expressed as

ϕ(ξ)=A1cos⁡(α1ξ)+A2sin⁡(α1ξ)+A3cosh⁡(β1ξ)+A4sinh⁡(β1ξ)(14)

ψ(η)=B1cos⁡(α2η)+B2sin⁡(α2η)+B3cosh⁡(β2η)+B4sinh⁡(β2η)(15)

where A1−A4 and B1−B4 are the mode function coefficients. It should be emphasized that the mode functions in all SOV methods have the same expressions as given in Eqs. (8), (14) and (15), but each SOV method has a different technique to obtain α1, β1, α2, β2 and Ω. By substituting ϕ and ψ into the BCs of two pairs of opposite edges in Table 1, one can achieve A1−A4 and B1−B4 and two transcendental eigenvalue equations about α1, β1, α2, β2 and Ω. Solving for these five unknowns needs the other three equations, which are obtained below.

With the substitution of the four combinations (μ1,2,λ1,2), (μ1,2,λ3,4), (μ3,4,λ1,2), (μ3,4,λ3,4) into Eq. (10), one can have

α12+(αα2)2=Ω2(16)

−α12+(αβ2)2=±Ω2(17)

β12−(αα2)2=±Ω2(18)

β12+(αβ2)2=Ω2(19)

If the right-hand terms of Eqs. (17) and (18) are all equal to Ω2, or −Ω2, then using Eqs. (16)–(18) yields β12+α2β22=3Ω2, or β12+α2β22=−Ω2. These two equations contradict with Eq. (19), and this situation is denoted as Case 1. If the signs of the right-hand terms of Eqs. (17) and (18) are opposite, then Eq. (17) plus Eq. (18) and using (16) lead to Eq. (19), denoting this situation as Case 2.

Therefore, it can be seen that in Case 2, β12+α12=0 or β22+α22=0, having no non-trivial solutions, thus Case 2 is unrealistic and not considered further. In Case 1, the right-hand terms of Eqs. (17) and (18) should be equal to Ω2 since β12+α2β22=3Ω2 is more reasonable than β12+α2β22=−Ω2 in terms of Eq. (19), so Eqs. (17) and (18) should be

(αβ2)2=Ω2+α12(20)

β12=Ω2+(αα2)2(21)

The Eqs. (16), (20) and (21) are the three required equations. By using Eqs. (16), (20), (21) and two transcendental eigenvalue equations, one can determine α1, β1, α2, β2 and Ω. Besides, according to Eqs. (16), (20) and (21), we have

α12+β12=2Ω2(22)

(αβ2)2+(αα2)2=2Ω2(23)

Eq. (22) denotes the relation between Ω and two x-direction eigenvalues (α1, β1), and Eq. (23) stands for the relation between Ω and two y-direction eigenvalues (β2, α2).

For the problems of Category 1, Eq. (16) and two transcendental eigenvalue equations are solved together, and the solutions satisfy Eq. (7) and BCs exactly. For the problems of Category 2, Eqs. (16) and (20) or Eqs. (16) and (21) are solved together with two transcendental eigenvalue equations, and the results also satisfy Eq. (7) and BCs exactly. Therefore, the Navier method and the Levy method are the two special cases of the SOV methods.

For the problems of Category 3, three Eqs. (16), (20), (21) and two transcendental eigenvalue equations are solved concurrently, but Eq. (7) is not satisfied since Eq. (19) is not satisfied, implying that the obtained closed-form mode function w and the frequency Ω are not exact. There are two SOV methods solving Eqs. (16), (20), (21) and two transcendental eigenvalue equations. For the case without free edges, the constructed SOV method is the direct SOV method [11]. For the case with free edges, the constructed SOV method is the imSOV method [13], in which the expressions of separable free BCs are given in Section 2.3.

2.3 Rayleigh Quotient

The Rayleigh quotient is generally used to derive characteristic PDEs of continuous systems and generalized eigenvalue equations of discrete systems and also serves as the foundation of FEM for dynamics problems. However, the present authors and coworkers found that the closed-form SOV analytical solutions can also be achieved through the Rayleigh quotient.

For free vibration, the Rayleigh quotient [36] has the form as

ω2=stUmaxT0(24)

where ‘st’ denotes ω2 is the stationary value of the quotient (Umax/T0); Umax is the maximum potential energy of a harmonic vibration; ω2T0 is the maximum kinetic energy; T0 is called the reference kinetic energy or the kinetic energy coefficient. For thin plates, we have

Umax=D2∬[(∂2w∂x2)2+2ν∂2w∂x2∂2w∂y2+(∂2w∂y2)2+2(1−ν)(∂2w∂x∂y)2]dxdy(25)

T0=12∬ρhw2dxdy(26)

By using the non-dimensional coordinates ξ and η, Eqs. (25) and (26) can be transformed into

Umax=abD2∬[1a4(∂2w∂ξ2)2+2νa2b2∂2w∂ξ2∂2w∂η2+1b4(∂2w∂η2)2+2(1−ν)a2b2(∂2w∂ξ∂η)2]dξdη(27)

T0=ab2ρh∬w2dξdη(28)

With the substitution of Eqs. (27) and (28) into Eq. (24), one can obtain the same PDE as Eq. (7) and the same natural BCs as those in Table 1.

When seeking the SOV solutions based on the Rayleigh quotient, Eqs. (8), (9) and Eqs. (12)–(15) are also used, and we assume that ψ is known but ϕ is unknown, or ψ is unknown but ϕ is known. When ψ is known and ϕ is unknown, Ωx is used to denote the corresponding frequency. When ψ is unknown and ϕ is known, Ωy is used to represent the corresponding frequency.

For the case with known ψ and unknown ϕ, according to Eqs. (27), (28) and (24), we can obtain the governing equation for ϕ, as

d4ϕdξ4+2α2[νS2S1−(1−ν)S3S1]d2ϕdξ2+(α4S4S1−Ωx4)ϕ=0(29)

and the boundary bending moment and equivalent shear force for the edges ξ=−1 and ξ=1, as

Mξ=abDa4(S1d2ϕdξ2+να2S2ϕ)Vξ=abDa4[S1d3ϕdξ3+α2(νS2−2(1−ν)S3)dϕdξ](30)

S1=∫−11ψ2dη,S2=∫−11(ψd2ψdη2)dη,S3=∫−11(dψdη)2dη,S4=∫−11(d2ψdη2)2dη(31)

If an edge parallel to the y-axis is free, the separable free BCs can be obtained from Eq. (30), as

d2ϕdξ2+να2S2S1ϕ=0d3ϕdξ3+α2(νS2S1−2(1−ν)S3S1)dϕdξ=0(32)

By substituting the first expression in Eq. (9), or ϕ(ξ)=Aeμξ into Eq. (29), one can obtain the relation between μ and Ωx as

μ4+2α2[νS2S1−(1−ν)S3S1]μ2+(α4S4S1−Ωx4)=0(33)

the roots of which can be written in the same forms as those in Eq. (12), but here the spatial eigenvalues are

α1=α[νS2S1−(1−ν)S3S1]2−S4S1+Ωx4α4+[νS2S1−(1−ν)S3S1](34)

β1=α[νS2S1−(1−ν)S3S1]2−S4S1+Ωx4α4−[νS2S1−(1−ν)S3S1](35)

The expression of ϕ in terms of α1 and β1 has the same form as that in Eq. (14). Similarly, for the situation with unknown ψ and known ϕ, the governing equation for ψ can be achieved as

d4ψdη4+2α−2[νT2T1−(1−ν)T3T1]d2ψdη2+α−4(T4T1−Ωy4)ψ=0(36)

where

T1=∫−11ϕ2dξ,T2=∫−11(ϕd2ϕdξ2)dξ,T3=∫−11(dϕdξ)2dξ,T4=∫−11(d2ϕdξ2)2dξ(37)

The separable boundary bending moment and equivalent shear force for the edges η=−1 and η=1 are

Mη=abDb4(T1d2ψdη2+νT2α2ψ)Vη=abDb4[T1d3ψdη3+νT2−2(1−ν)T3α2dψdη](38)

The separable free BCs of an edge parallel to the x-axis can be obtained from Eq. (38), as

d2ψdη2+νT2α2T1ψ=0d3ψdη3+νT2−2(1−ν)T3α2T1dψdη=0(39)

With the substitution of the second expression in Eq. (9), or ψ(η)=Beλη into Eq. (36), the relation between λ and Ωy is obtained as

λ4+2α−2[νT2T1−(1−ν)T3T1]λ2+α−4(T4T1−Ωy4)=0(40)

whose roots have the forms of Eq. (13), and

α2=α−1[νT2T1−(1−ν)T3T1]2−T4T1+Ωy4+[νT2T1−(1−ν)T3T1](41)

β2=α−1[νT2T1−(1−ν)T3T1]2−T4T1+Ωy4−[νT2T1−(1−ν)T3T1](42)

The expression of ψ in terms of α2 and β2 is the same as that in Eq. (15).

The SOV methods based on the Rayleigh quotient include the iSOV method [12], the eSOV method [14], and the vSOV method [12]. In these three methods, the coefficients A1−A4 and B1−B4 and two transcendental eigenvalue equations can also be achieved by substituting ϕ and ψ into the simply support, clamp and guide BCs in Table 1 and the free BCs in Eqs. (32) and (39).

In the iSOV and eSOV methods, two transcendental eigenvalue equations and Eqs. (34), (35), (41), (42) are solved for α1, β1, α2, β2, Ωx and Ωy. In the vSOV method, Ωx=Ωy=Ω, and two transcendental eigenvalue equations and Eqs. (34), (35), (23), or two transcendental eigenvalue equations and Eqs. (41), (42), (22) are solved for α1, β1, α2, β2, Ω.

The iSOV and the eSOV methods are suitable for solving three Categories of problems. For the problems of Categories 1 and 2, the obtained solutions are exact; for the problems of Category 3, the solutions make the Rayleigh quotient take stationary value, implying that the solutions are the most accurate in the SOV function space. The vSOV method applies to the eigenproblem analysis of the rectangular plates without adjacent free edges.

3  SOV Methods and Bisection-Based Solution Procedures

In Section 2, the SOV methods have been formed based on the theoretical framework. This section first gives the purposes of proposing each SOV methods and then presents the corresponding eigenvalue equations. Then, the bisection-based solution procedures are provided for the imSOV method, the vSOV method, and the eSOV method respectively. Unlike the Newton iteration method, the bisection-based methods do not have the problem of choosing the initial values of solutions. The computational cost is discussed finally.

3.1 Improved SOV Method

The imSOV method [13] is proposed to solve characteristic partial differential Eq. (7), which is a general method since it is capable of dealing with arbitrary homogeneous BCs. If an edge is free, its separable BCs are given in Eqs. (32) and (39).

This method uses Eqs. (16), (20), (21) and two transcendental equations to solve α1, β1, α2, β2 and Ω. The two transcendental equations can be attained through the substitution of ϕ(ξ) and ψ(η) into two pairs of BCs, as shown in Section 2.2. Here, taking a CCSS rectangular plate as an example, the two transcendental equations can be derived as

α1tanh⁡2 β1=β1tan⁡2 α1(43)

α2tanh⁡2 β2=β2tan⁡2 α2(44)

The bisection-based solution procedure for the imSOV method is as follows:

Preparation

1) To express α2 and β2 in terms of α1 and β1.

Eliminating Ω from Eqs. (16) and (21) generates

α22=β12−α122α2(45)

Then, with Eqs. (16), (20) and (45), one can have

β22=β12+3α122α2(46)

2) Substituting Eqs. (45) and (46) into Eq. (44) yields

β12−α12tanh⁡2(β12+3α12)α2=β12+3α12tan⁡2(β12−α12)α2(47)

Therefore, the original five Eqs. (16), (20), (21), (43) and (44) reduce to two nonlinear equations Eqs. (43) and (47), and the unknowns are α1 and β1.

Solution procedure

The flowchart of solving Eqs. (43) and (47) together is presented in Figs. 2 and 3, in which the non-zero small number ε is the initial value of α1; Δ is a small increment, and f1(α1,β1)=0 denotes Eq. (47) and f2(α1,β1)=0 stands for Eq. (43).

Figure 2: The procedure finding solution domains in the imSOV method

Figure 3: The process of calculating accurate results in each solution domain

Step 1 Find solution domains

The flowchart in Fig. 2 is used to search solution domains of Eq. (47), and a domain is represented by α1∈(α1bot,α1up), β1∈(β1bot,β1up).

It can be seen that by increasing the value of α1 gradually, one can find a set of (α1,β1) satisfying Eq. (43); then, one can select the solution domains of Eq. (47) from the set of (α1,β1) satisfying Eq. (43).

Step 2 Calculate accurate results

For each domain, the procedure shown in Fig. 3 is employed to calculate accurate solutions with the bisection method, in which δ, the solution precision of α1 and β1, is a number much smaller than Δ. Note that δ is also the solution precision when using the bisection method, for example, to solve β1bot for a given α1bot in Fig. 2. Finally, α2, β2 and Ω can be calculated through Eqs. (45), (46), and one of Eqs. (16), (20), (21).

It is noteworthy that the above solution procedure including two steps, finding solution domains and calculating accurate solutions, is also used in other SOV methods. If the imSOV method is used to deal with the plates with free corners, the solution procedure is the same as that of the eSOV method, refer to Section 3.3.

3.2 Variational SOV Method

The vSOV method [12] is proposed to improve accuracy and extend the application of the direct SOV method [11]. The vSOV method can deal with rectangular plates without free corners (a free corner implies two adjacent edges are free). In addition to two transcendental eigenvalue equations, this method assumes Ωx=Ωy=Ω and solves either Eqs. (22), (41) and (42) or Eqs. (23), (34) and (35).

It can be seen that Eqs. (22) and (23) comes from Eq. (7) (a strong-form equation) and Eqs. (34) and (35), or (41) and (42) comes from Eq. (24) (a weak-form equation). So the vSOV method can be viewed as a mixed method, implying that the vSOV method is used to solve the strong-form Eq. (7) in the x or the y direction and obtain the stationary value of the Rayleigh quotient in the y or the x direction simultaneously.

It should be noteworthy that the vSOV method is the basis of the iSOV method [12] and the eSOV method [14]. Based on the bisection method, the solution procedures of the vSOV method are presented below for two cases.

Case1: Plates without free edges

Also take a CCSS plate as an example, and the transcendental equations are Eqs. (43) and (44). The coefficients A1−A4 of ϕ(ξ) are

A2=A1cot⁡α1,A3=−A1cos⁡α1cosh⁡β1,A4=−A1cos⁡α1sinh⁡β1(48)

Preparation

1) To express α2, β2 and Ω in terms of α1 and β1.

As shown in Eqs. (14) and (48), ϕ(ξ) is the explicit function of α1 and β1, then Eq. (37) shows that T1−T4 are the explicit functions of α1 and β1.

Besides, Eq. (22) shows that Ω is also the function of α1 and β1, so Eqs. (41) and (42) indicate that α2 and β2 are the explicit functions of α1 and β1.

2) Substituting Eqs. (41) and (42) into Eq. (44), one can have a transcendental equation about α1 and β1, which can be solved together with Eq. (43) for α1 and β1.

Solution procedure

The solution procedure is the same as that in Section 3.1, refer to Figs. 2 and 3.

Case2: Plates with one free edge or two opposite free edges

Here taking a CSSF plate as an example, the transcendental equation for the x direction is Eq. (43). According to the BCs of the edges η=−1 and η=1, the y-direction transcendental equation can be achieved as

β2(Ny2−α22)(Ny1+β22)tan⁡2 α2=α2(Ny2+β22)(Ny1−α22)tanh⁡2 β2(49)

where Ny1 and Ny2 are

{Ny1=νT2−2(1−ν)T3α2T1Ny2=νT2α2T1(50)

Preparation

1) To express α2, β2 and Ω in terms of α1 and β1.

Eq. (37) shows T1−T4 are the explicit functions of α1 and β1, and Ω can be calculated from Eq. (22) using α1 and β1, so α2 and β2 are the explicit functions of α1 and β1, refer to Eqs. (41) and (42).

2) By substituting Eqs. (41) and (42) into Eq. (49), one can obtain a transcendental equation regarding α1 and β1, which can be solved together with Eq. (43) for α1 and β1.

Solution procedure

The solution procedures of the vSOV method are the same as that of the imSOV method in Section 3.1.

3.3 Extended SOV Method

The eSOV method [14] is proposed to improve accuracy and extend the application of the vSOV method [12], and it applies to any homogeneous BCs, so it is a general solution method. The eSOV method obtains the closed form eigensolutions by finding the stationary values of the Rayleigh quotient.

The method has nothing to do with Eq. (7). In this method, four Eqs. (34), (35), (41) and (42) in conjunction with two transcendental eigenvalue equations are solved to achieve the six unknowns α1, β1, α2, β2, Ωx and Ωy, in which Ωx and Ωy are independent each other in mathematical sense, but they are the same in physical sense, which can be proved through the Rayleigh quotient itself.

To simplify solution procedures, the four unknowns of α1, β1, α2 and β2 are solved first. By eliminating Ωx from Eqs. (34) and (35), we have

α12−β12=2α2[νS2S1−(1−ν)S3S1](51)

and eliminating Ωy from Eqs. (41) and (42) generates

α22−β22=2α−2[νT2T1−(1−ν)T3T1](52)

Then we only need to solve Eqs. (51), (52) and two transcendental equations for α1, β1, α2 and β2. After obtaining α1, β1, α2 and β2, one can use any one of Eqs. (34), (35), (41) and (42) to obtain the frequency Ω because of Ω=Ωx=Ωy.

Case1: Plates without free edges

Also taking the CCSS plate as an example, and the two transcendental equations are Eqs. (43) and (44). The mode function coefficients A1−A4 are given in Eq. (48), and B1−B4 are

B2=B1cot⁡α2,B3=−B1cos⁡α2cosh⁡β2,B4=−B1cos⁡α2sinh⁡β2(53)

Preparation

1) To express β2 in terms of α1, β1, and α2.

Eq. (52) can be reformed as

β22=α22−2α−2[νT2T1−(1−ν)T3T1](54)

which shows that β2 is an explicit function of α1,β1, and α2. Then ψ(η) in Eq. (15) and S1−S4 in Eq. (31) are also the explicit functions of α1,β1,α2.

2) By substituting Eq. (54) into Eq. (44), we have

α2tanh⁡2α22−2α−2[νT2T1−(1−ν)T3T1]=α22−2α−2[νT2T1−(1−ν)T3T1]tan⁡2α2(55)

Since T1−T4 are the explicit functions of α1 and β1, thus the y-direction transcendental Eq. (55) is an equation for α1,β1 and α2,

Solution procedure

After preparation, we can solve Eqs. (55), (43) and (51) for α1,β1,α2. Eqs. (43) and (55) are the transcendental equations, and Eq. (51) denotes the relation between the two-direction eigenvalues.

The procedure to find solution domains of Eqs. (55), (43) and (51) is presented in Fig. 4, wherein Eq. (43) is denoted by f2(α1,β1)=0, Eq. (55) is denoted by f1(α1,β1,α2)=0 and Eq. (51) is denoted by f3(α1,β1,α2)=0. The process of calculating accurate results is similar to that in Fig. 3. More explanations about Fig. 4 are given below:

1) Since Eq. (43) is only about α1 and β1, so to facilitate the solution procedure, for a given α1bot Eq. (43) is solved firstly to obtain β1bot with the bisection method.

2) With the known (α1bot,β1bot), Eq. (55) is a nonlinear equation of α2, which can be solved by the bisection method to obtain α2bot.

3) After having a group of (α1bot,β1bot,α2bot) satisfying Eqs. (43) and (55). By gradually increasing the value of α1 with the increment Δ, one can determine a group of (α1up,β1up,α2up), then the solution domain of Eq. (51) is obtained, which is α1∈(α1bot,α1up), β1∈(β1bot,β1up) and α2∈(α2bot,α2up).

Figure 4: The procedure finding solution domains in the eSOV method

Case2: Plates with one free edge or two opposite edges free

The situation with one free edge is considered first. Similar to the vSOV method, the CSSF plate is taken as an example too. Eqs. (43) and (49) are the two transcendental equations. The coefficients A1−A4 are given in Eq. (48), and the coefficients B1−B4 are

B2=B1cot⁡α2,B3=−B1cos⁡α2cosh⁡β2Ny2−α22Ny2+β22,B4=−B1cos⁡α2sinh⁡β2Ny2−α22Ny2+β22(56)

Preparation

1) As shown in Section 3.1, β1 can be achieved from Eq. (43) for a given α1 with the bisection method, and T1−T4 are the explicit functions of α1 and β1.

2) Ny1 and Ny2 are also the explicit functions of α1 and β1 according to Eq. (50), from which we have

Ny1+Ny2=2α−2[νT2T1−(1−ν)T3T1](57)

From Eqs. (52) and (57), one can obtain

β22=α22−Ny1−Ny2(58)

Substituting Eq. (58) into Eq. (49) to eliminate β2 generates

α22−Ny1−Ny2(Ny2−α22)2tan⁡2α2=α2(Ny1−α22)2tanh⁡(2α22−Ny1−Ny2)(59)

which is a nonlinear equation for α1,β1,α2. Then from Eqs. (58), (56), (15) and (31), one can see that β2, B1−B4, ψ(η) and S1−S4 are also the explicit functions of α1,β1,α2.

Solution procedure

Solving Eqs. (43), (59) and (51) can give α1,β1,α2. Among these three equations, Eqs. (43) and (59) are the transcendental equations corresponding to two coordinate directions respectively, and Eq. (51) establishes the eigenvalue relation of two directions. The process for finding the solution domains of these three equations is the same as that shown in Fig. 4. Then a similar process as that in Fig. 3 can be used to calculate accurate results.

Case3: Plates with free corners

Taking a SSFF plate as an example, the y-direction transcendental equation is Eq. (49), and the x-direction transcendental equation is similar to Eq. (49), as

β1(Nx2−α12)(Nx1+β12)tan⁡2α1=α1(Nx2+β12)(Nx1−α12)tanh⁡2 β1(60)

where Nx1 and Nx2 are the variables from the free BCs of the edges ξ=−1 and ξ=1, as

Nx1=α2νS2−2(1−ν)S3S1(61)

Nx2=α2νS2S1(62)

The mode function coefficients B1−B4 are listed in Eq. (56), and A1−A4 are

A2=A1cot⁡α1,A3=−A1cos⁡α1cosh⁡β1Nx2−α12Nx2+β12,A4=−A1cos⁡α1sinh⁡β1Nx2−α12Nx2+β12(63)

Preparation

1) Deal with the eigenvalue equations in the x direction.

From Eqs. (61) and (62), we have

Nx1+Nx2=2α2[νS2S1−(1−ν)S3S1](64)

Then from Eqs. (51) and (64), one can obtain

α12−β12=Nx1+Nx2(65)

By substituting Eq. (65) into Eq. (60), we have

β1(Nx2−α12)2tan⁡2α1=α1(Nx2+β12)2tanh⁡2 β1(66)

which is a quadratic equation of Nx2. By using the root formula, Nx2 can be obtained, which is the explicit function of α1,β1, then the coefficient A1−A4 and the integrals T1−T4 are the explicit functions of α1,β1.

2) Deal with the eigenvalue equations in the y direction, and the work is the same as that in Section 3.2.

3) From Eqs. (65) and (66), we know that both Nx1 and Nx2 are the explicit functions of α1,β1. It can be seen from 2) that S1−S4 are the explicit functions of α1,β1andα2. By substituting S1−S4,Nx1,Nx2 into Eqs. (61) and (62), we have two nonlinear equations for α1,β1,α2.

Solution procedure

In the above preparation, one has obtained three nonlinear equations, including Eqs. (59), (61) and (62). Among them, Eq. (59) is the y-direction transcendental equation, and Eqs. (61) and (62) relate the eigenvalues of the x and y directions.

The solution domains of these three equations can be found with the procedure in Fig. 5, in which f1=0, f2=0 and f3=0 stand for Eqs. (59), (62) and (61), respectively. By increasing the value of α1 gradually, one can find a set of (α1,β1,α2) satisfying Eqs. (59) and (62); then by using the bisection method again, one can find (α1,β1,α2) satisfying Eq. (61) from the set of (α1,β1,α2). Finally, the closed-form solutions are achieved.

Figure 5: The procedure finding solution domains for SSFF plates in the eSOV method

In addition to the above solution procedures, an iteration method can also be used here. If (α1, β1, Ωx) are obtained with the assumed (α2, β2, Ωy), then (α2, β2, Ωy) can be updated with the obtained (α1, β1, Ωx), then updating (α2, β2, Ωy) again, the iteration ends if |Ωy−Ωx|/Ωx is small enough. This iteration solution method is the iSOV method [12].

3.4 Computational Cost

The SOV methods are analytical methods, which can be used to solve characteristic PDEs or the stationary value of the Rayleigh quotient, not requiring to discretize spatial domain like FEM. To obtain frequencies and modes, only a few closed-form transcendental and algebraic eigenvalue equations need to be solved, as shown in Sections 3.1–3.3. So, the SOV methods are efficient and economical, and solving programs can be implemented on regular PCs.

In previous publications (for example, [11–15]), the Newton iteration method was used to solve eigenvalue equations. If initial values are appropriate, convergence can be reached after 3 or 4 iterations. But the initial values of solutions should be close to solutions, otherwise the solutions cannot be found, or the solutions are not expected, implying that selecting initial values is a troublesome issue especially when solving solutions of many orders.

To avoid the issue of selecting initial values, the bisection-based solution procedures are presented in this work as alternative methods to the Newton method, as shown in Sections 3.1–3.3. The procedures include two steps, searching solution intervals and calculating accurate solutions, and the bisection approach is used in the two steps.

4  Nodal Line Patterns and the Solution Integrity

This section gives the nodal line equations, from which the nodal line geometrical configurations or patterns can be obtained. Besides, with the help of the nodal line patterns, the solution integrity of the SOV methods is qualitatively discussed.

4.1 Nodal Line Patterns

The SOV methods have the advantage that they have separable and explicit equations of nodal lines. The nodal line equations can be achieved by assuming the mode function w(ξ,η) in Eq. (8) is equal to zero, that is

w(ξ,η)=ϕ(ξ)ψ(η)=0(67)

then the nodal line equations are

ϕ(ξ)=0(68)

ψ(η)=0(69)

Solving the above two equations yields the positions of nodal lines within a plate. Since ϕ(ξ) is just a function of ξ and ψ(η) is just a function of η, then the nodal lines determined by Eqs. (68) and (69) are straight lines parallel to plate edges. That is, the nodal lines of ϕ(ξ)=0 are perpendicular to the x-axis, and the ones of ψ(η)=0 perpendicular to the y-axis. Besides, at the intersections of the nodal lines, ϕ(ξ)=ψ(η)=0. In the SOV methods, Eq. (68) and Eq. (69) are applicable to arbitrary homogeneous BCs.

According to Eqs. (14), (15), (68) and (69) can be rewritten as

ϕ(ξ)=A1cos⁡(α1ξ)+A2sin⁡(α1ξ)+A3cosh⁡(β1ξ)+A4sinh⁡(β1ξ)=0(70)

ψ(η)=B1cos⁡(α2η)+B2sin⁡(α2η)+B3cosh⁡(β2η)+B4sinh⁡(β2η)=0(71)

where A1−A4, B1−B4, α1,α2, β1,β2 can be obtained through the procedures provided in Section 3, so the forms of these two equations are explicit. A mode of a plate has no or many nodal lines, and the feature of a mode shape depends on the positions of nodal lines. The nodal line pattern of a mode is a geometry graph including all nodal lines inside a plate.

In the following, a simply supported plate is considered to show how to determine the positions of nodal lines. According to Eqs. (14) and (15) and the simply supported BCs, see Table 1, we can obtain the eigenvalue equations and mode function coefficients as follows:

sin⁡2 α1=0sin⁡2 α2=0}⇒α1=mπ2α2=mπ2}(72)

A1=sin⁡α1,A2=cos⁡α1,A3=A4=0

B1=sin⁡α2,B2=cos⁡α2,B3=B4=0(73)

where m,n=1,2,3,⋯ are the numbers of half waves in the x and y directions, respectively. Then using Eqs. (70) and (71) generates the equations of nodal lines as

ϕ(ξ)=sin⁡mπ2cos⁡(mπ2ξ)+cos⁡mπ2sin⁡(mπ2ξ)=sin⁡mπ2(1+ξ)=0(74)

ψ(η)=sin⁡nπ2cos⁡(nπ2η)+cos⁡nπ2sin⁡(nπ2η)=sin⁡nπ2(1+η)=0(75)

where ξ∈(−1,1), η∈(−1,1). Given m and n, one can calculate the nodal line positions and the number of nodal lines from Eqs. (74) and (75). For example, if m = 2, then ξ=0 and there is only one nodal line; if m = 3, then ξ=±1/3 and there are two nodal lines. The number of nodal lines is equal to the number of half waves minus one.

The nodal line patterns of the SSSS and CCSS isotropic square plates (α=1) are shown in Figs. 6 and 7 respectively, where the values (i, j) = (m – 1, n – 1), wherein i and j represent the numbers of nodal lines perpendicular to the x-axis and the y-axis, respectively. By observing the nodal line patterns and the frequency orders, one can find that for these two square plates, the correspondences between the numbers of nodal lines and the frequency orders are the same, but the positions of nodal lines are different. For the SSSS plate, the nodal lines distribute evenly within the plate, while for the CCSS plate, the nodal lines no longer distribute evenly, further away from the clamped edges parallel to them, refer to the two patterns of (2, 2) in Figs. 6, 7. It can be seen that the sizes of square ① and square ② in Fig. 6 are the same, while in Fig. 7 square ① is larger than square ②.

Figure 6: Nodal line patterns of the SSSS square plate for the first 15 modes

Figure 7: Nodal line patterns of the CCSS square plate for the first 15 modes

For comparison, Tables 2, 3 present the relationships between the numbers of nodal lines and the frequency orders for the SSSS and CCSS rectangular plates (α=1.5) respectively. For the two rectangular plates, the correspondences between the nodal line patterns and frequency orders are different. For example, the frequency of order 15 for the SSSS plate corresponds to the pattern (0, 3), while the frequency of the same order for the CCSS plate corresponds to the pattern (5, 0).

Besides, Figs. 6, 7 show that the nodal lines in the SOV methods are straight, even though for repeated frequencies, such as the two frequencies with patterns of (1, 0) and (0, 1). But the combinations of two modes having straight nodal lines that are parallel to plate edges can produce the modes having diagonal or curved nodal lines, also refer to the modes plotted in Section 5. Here the SSSS square plate is taken as an example to show the combining results. For clarity’s sake, the first pair of repeated frequencies corresponding to the patterns of (1, 0) and (0, 1), or (m, n) = (2, 1) and (1, 2), is taken into account below. The two-mode functions are

w(ξ,η)|(m,n)=(2,1)=ϕ(ξ)|m=2ψ(η)|n=1=−sin⁡(πξ)cos⁡(π2η)w(ξ,η)|(m,n)=(1,2)=ϕ(ξ)|m=1ψ(η)|n=2=−sin⁡(πη)cos⁡(π2ξ)(76)

Then by letting the two combinations of the two modes equal to zero, we have

w(ξ,η)|(m,n)=(2,1)+w(ξ,η)|(m,n)=(1,2)=−4cos⁡πη2cos⁡πξ2cos⁡π(ξ−η)4sin⁡π(ξ+η)4=0w(ξ,η)|(m,n)=(2,1)−w(ξ,η)|(m,n)=(1,2)=−4cos⁡πη2cos⁡πξ2cos⁡π(ξ+η)4sin⁡π(ξ−η)4=0(77)