Forced free vibrations of a square plate

The vibrations of a square plate with free edges that is forced at its center has been investigated experimentally and numerically. A two-step finite volume scheme is proposed and is tested against an exact solution to a related problem. The governing equation and numerical solution procedure were constructed to closely mimic the experiment. Good agreement was found between the numerical simulations and the experiments.


Introduction
Owing to the numerous engineering applications there have been a plethora of investigations devoted to vibrating rectangular plates since the original mathematical approach formulated by Euler in 1766. The documented studies in the literature are thoroughly summarized in the works of Leissa [1,2] for the various distinct combinations of boundary conditions including clamped, simply-supported or free. Interest and research in vibrating plates gained a lot of attention after the German physicist and musician Ernst Chladni, who lived from 1756 to 1827, conducted his famous experiments [3]. Chladni is often referred to as the father of acoustics for his ground-breaking work on vibrating plates. He visualized the patterns formed when sand was sprinkled on a vibrating horizontal plate. Chladni observed that when the plate was allowed to resonate the sand particles formed stationary nodal lines which are now known as Chladni figures. He would demonstrate his discovery by strumming his bow along the edge of a plate and watching the sand dance on the surface. The sand would eventually settle along the nodal lines which are places on the plate which remained stationary. The complex patterns changed as the frequency of vibration was varied [4]. Apparently, very fine particles were later observed to settle along the anti-nodes, forming what are known as inverse patterns, due to air currents that are induced by the vibrating plate which drag the particles to the anti-nodes. In a vacuum it was believed that all particles would migrate to the nodes regardless of their size. However, Gerner et al. [5] showed that air currents are not the only mechanism leading to inverse patterns. They discovered that when the acceleration of the resonating plate did not exceed the acceleration due to gravity, the particles would always roll to the anti-nodes, irrespective of their size, and even in the absence of air.
Chladni also noticed that for a circular plate f ∼ (m + 2n) 2 where f is the frequency, m is the number of nodal diameters and n is the number of nodal circles. This result is known as Chladni's Law [6] and Kirchhoff [7] showed that this is an asymptotic relation valid for large f. The quest to uncover the theory behind the Chladni figures was ignited by Emperor Napolean when he observed a live demonstration delivered by Chladni during a presentation in Paris in the year 1809 [8]. In the same year he commissioned the French Academy of Science to offer a prize to whomever could formulate the mathematical theory that could explain the observed patterns. Famous mathematicians including Lagrange, Laplace, Legendre and Poisson served on the judging panel. In the first competition only one submission was received; it was from Sophie Germain, a French mathematician who lived from 1776 to 1831 and made important contributions in number theory as well as the theory of elasticity [9,10]. Unfortunately, errors were found in her derivation and no prize was awarded. A second competition was announced. Again, only one submission was received and it was from Sophie Germain. This time she received an honourable mention but still no prize because there were mistakes in her analysis. A third competition then took place. Finally, in 1816 Sophie Germain received the prize and was the first woman to receive an award of this importance from the French Academy. Her contributions to this problem have been published in the following articles [11,12]. To recognize her achievement the French Academy established the Sophie Germain Prize in her honour and it has been awarded annually since 2003. Final corrections and improvements to the mathematical formulation were made by Lagrange, Poisson and Kirchhoff [7].
As a first approximation one can model the plate as a tightly stretched thin elastic membrane. If we further assume that the deflection from equilibrium, denoted by u, is small, the tension ( ) and membrane density ( ) are uniform, and there is no resistance to motion, then u(x, y, t) satisfies the two-dimensional wave equation given by Here, c refers to the wave speed which is related to the tension and mass density through c 2 = ∕ . For free edges we impose the following boundary conditions By separation of variables the general solution is easily found to be where with U 0 denoting the initial membrane displacement while U 0 refers to the initial velocity of the membrane. The angular frequencies mn are given by mn = c a √ m 2 + n 2 . Note that for a specified angular frequency mn there may be one or more distinct modes of vibration having the same angular frequency. For example, 12 = 21 = √ 5 c∕a and the modes are Hence, the general spatial form of the membrane deflection for this frequency will be given by where C, D are constants. Because of symmetry we expect that C = D and hence the nodal lines satisfy [13] The nodal lines are plotted in Fig. 1 (left) and are compared to the observed experimental pattern (right) corresponding to f = 190Hz . As explained in a later section, the experiments were carried out using a square plate which cos m x a cos n y a A mn cos( mn t) + B mn sin( mn t) ,  is forced at its center to vibrate at a set frequency. Apart from the diamond shape in the left-hand plot the agreement is remarkable. Unfortunately, the wave equation model is highly simplified and cannot exactly reproduce the observed patterns.
The main goal of this investigation is to bring the experimental observations in harmony with theory. The previous studies by Gander and Kwok [14] and Gander and Wanner [15] tackled this problem by solving an associated eigenvalue problem using the Ritz method [16] to compute the corresponding eigenpairs. As noted by Kirchhoff [7], the eigenvalues are related to the resonant frequencies while the eigenfunctions represent the nodal patterns. The approach adopted here is significantly different and new; to our knowledge the proposed solution procedure has not been carried out before. As we will see in Sect. 2 the equation governing the motion of the plate and the accompanying boundary conditions are far more complicated than the above wave equation model. Since an exact solution to the equation is still out of reach a numerical solution procedure that closely mimics the experiment is outlined in Sect. 3. So rather than solving an eigenvalue problem the governing equation will be solved with special attention paid to enforcing the non-standard boundary conditions. The experimental and numerical results are presented and discussed in Sect.4. Finally, in Sect. 5 the key findings of the investigation are summarized.

Mathematical formulation
We consider the transverse deflections of a square plate. The plate is taken to be composed of an isotropic, homogeneous material and is assumed to be thin; that is, the side length, 2a, is much larger than the uniform thickness, l. The material properties of the plate include the density (mass per unit volume), , Young's modulus, q, and the Poisson ratio, . A Cartesian coordinate system (x, y, z) is adopted whereby the centre of the plate coincides with the origin, the transverse deflections occur in the z direction, and the edges of the plate are located at x = ±a and y = ±a . For small deflections from the equilibrium position z = 0 linear elasticity theory predicts that the displacement of the plate, denoted by u(x, y, t), will satisfy the following partial differential equation given by [17] where D 2 is the flexural rigidity defined by and ∇ 4 is the biharmonic operator given by For free vibrations equation (1) is to be solved subject to the following boundary conditions [17] We next cast the problem in dimensionless form. In order to achieve this we introduce H as a representative scale for the deflections and apply the following transformation in terms of dimensionless quantities identified by an asterisk With these scalings in place, and dropping the asterisks for notational convenience, the dimensionless governing equation and boundary conditions become Lastly, we also impose initial conditions which in dimensionless form are Here, U 0 refers to the initial plate displacement while U 0 is the initial velocity of the plate. Although the governing equation (8) is linear, what makes this problem particularly difficult are the free edge boundary conditions (9)(10)(11). As noted by Rayleigh [13], this is why an exact solution to this problem has not yet been found. This motivated Ritz [16] to formulate a (3) 2 u y x = 0 at the corners (x, y) = (±a, ±a) .
technique, known as the Ritz method, that allowed one to construct an approximate solution. He was able to produce the first numerical solution of Chladni figures on a square plate in 1909. It is worth noting that an exact solution to eq. (8) satisfying the boundary conditions can be found. Here, and it is easy to show that w also satisfies equation (8). The above boundary conditions are inspired by those used for solving the wave equation and the exact solution to eq. (8) subject to the boundary conditions (13)(14) and initial conditions (12) with U 0 = 0 is given by where The solution bears a close resemblance to that of the wave equation; the only difference lies in the frequency mn which is the square of the frequency found for the wave equation when cast in dimensional form. Further, the above solution also satisfies the corner conditions (11). This exact solution will be used to test our numerical scheme outlined in the following section.
Since the experiment was configured so that the plate is forced at its center, in order to mimic this we propose solving the equation where Q(x, y, t) = (x) (y) cos( 0 t) represents a forcing term concentrated at the origin which coincides with the center of the plate. Here, (x) and (y) denote the Dirac delta functions and 0 is the angular forcing frequency. Boundary conditions (9-11) and initial conditions (12) still apply. Approximations to the Dirac delta function with varying degrees of usefulness include: with similar expressions for n (y) . The actual expression used in our simulations will be discussed later in Section 4.2.

Numerical solution procedure
The finite volume method was used to solve equation (16). The square domain enclosed by the lines x = ±1 and y = ±1 was discretized and we have adopted the notation u k i,j to represent the numerical approximation to u( is the uniform grid spacing in the x, y directions and Δt is the time step. A two-step scheme was used to solve eq. (16). We first solve for w where and then solve the forced wave-type equation for u given by When discretizing eqs. (17 and 18) we must consider three cases: interior points, edge points and the corners. We will follow the approach implemented by Gander and Kwok [14] in their discretization of the associated eigenvalue problem and begin with an interior grid point. In this case all the grid points will lie inside the domain and the finite volume discretization will be identical to the standard finite differencing discretization. Using second-order central differencing this leads to the following explicit scheme [18] An advantage of using a finite volume method is that it is well suited to handle the singular forcing term. At all interior points excluding the central grid point the integral of Q over a control volume vanishes, while when Q is integrated over the central control volume we obtain cos( 0 t k ). We next consider an edge point shown in Fig. 2. We will present the details for the edge x = −1 with the understanding that the other edges are dealt with in a similar manner. Integrating equation (18) over the control volume yields where the first term is approximated as To evaluate the second term we apply the two-dimensional divergence theorem to obtain and approximate the fluxes over each segment of V ij . For example, As noted in [14], along Γ 4 we must incorporate the boundary condition as follows and from (9) and (17) we have that where The ghost points appearing in the above (i.e. having a subscript−1) can be determined by making use of the first condition in (9) which when discretized gives Making use of this expression the discretized form of eq. (19) along the edge becomes Lastly, we consider a corner shown in Fig. 3. Again, we will present the details for the corner located at (x, y) = (−1, −1) knowing that the remaining corners are handled in a similar fashion. This time the first term in (21) yields while the second term after repeating the calculations outlined above simplifies to Fig. 2 The control volume for an edge grid point Fig. 3 The control volume at a corner In arriving at this result we have made use of the corner condition (11) as well as the fact that w(±1, ±1, t) = 0 . The mixed derivative terms above can be approximated using formulae similar to those presented earlier.
We will next show that this scheme is stable provided Δt ≤ h 2 ∕2 by performing a stability analysis [19]. We consider the one-dimensional version of eq. (8) given by which in discretized form using second-order central differencing becomes Thus, in order for u k i to satisfy (22) the quantity in the brackets must vanish; solving for e − Δt yields Δt so for the solution u k i to remain bounded as t k → ∞ we require that which means that −1 ≤ ≤ 1 . In order to guarantee that this holds for all and Δx we must have that ≤ 1∕4 . Thus, which is also the convergence criterion for the one-dimensional heat equation

Results and discussion
We next present the experimental and numerical results obtained and discuss the comparisons.

Experimental results
The experiments were conducted using the PASCO Scientific Mechanical Wave Driver (model no. SF-9324) and the Chladni Plates Kit (model no. WA-9607). The wave driver was a standard electromagnetic device resembling a speaker with a steel rod attached to it, and the center of the plate was fastened to the rod. The aluminum square plate had a side length of 2a = 24cm and a thickness of l = 0.09cm ; the material properties of aluminum are listed in Table 1. Using these values the flexural rigidity, D 2 , given by eq. (2) has a numerical value of approximately 1.936 m 4 ∕s 2 and thus a 2 ∕D ≈ 0.010s . The experimental setup is illustrated in Fig. 4. The wave driver was connected to a Digital Function Generator/ Amplifier (model no. PI-9598) which forced the center of the plate to vibrate over frequencies ranging from 1 Hz up to 1 kHz. Following Chladni's technique sand was sprinkled on the plate in order to make the nodal lines visible. As the frequency of vibration was varied the plate was excited and various modes were observed revealing a variety of standing wave patterns. The observed patterns, which are readily available on the internet, are shown in Fig. 5. As expected, the patterns

Numerical results
In dimensionless form the problem is completely characterized by the Poisson ratio and the angular forcing frequency 0 . In our simulations we set = 0.33 for aluminum. The uniform grid spacing in the x, y directions was taken to be h = 0.01 and a constant time step of Δt = 0.000025 was used. In order to test the numerical scheme we made use of the exact solution given by eq. (15). First we solved (8) subject to the boundary conditions (13 and 14) using the initial conditions U 0 (x, y) = cos( x) cos(2 y) + cos(2 x) cos( y) ,U 0 = 0 at t = 0 .
As expected, the solution did not change with time and the nodal pattern was identical to that shown in Fig. 1 (left). Then we solved the forced eq. (16) using the boundary conditions (13 and 14) with a forcing frequency of 0 = 12 = 21 = 5 2 and initial conditions U 0 =U 0 = 0 . This time we observed the nodal pattern approach that shown in Fig. 1 as time advanced. Figure 6 displays the computed nodal patterns at t = 1, 5, 25 . Again, this is the expected outcome because we are forcing the plate at a resonant frequency, and although all the modes are excited, the modes having the same frequency will eventually dominate.
We also plotted the displacement of a corner of the plate and ran a long simulation. The output illustrated in Fig. 7 shows a beat pattern emerging; the left panel shows a blown up portion of the right panel and displays the rapid oscillations happening within the much longer oscillations taking place in the right panel. The occurrence of the beat pattern can be explained as follows. The discretized eqs. (19 and 20) actually approximates the coupled set of partial differential equations given by Since 1 ≫ 2 we can ignore the term 2 4 0 in the above, and so Hence, the resonant angular frequency of the discretized equation, 0 , is slightly different from the theoretical resonant angular frequency, mn = 2 (m 2 + n 2 ) , of the governing equation. Because of this, we should expect a beat pattern to form. Associated with the beat pattern are short and long periods, T short and T long , respectively, where 2 4 0 − 2 0 + 4 (m 2 + n 2 ) 2 − 2 1 6 (m 2 + n 2 )(m 4 + n 4 ) = 0 .
) . Substituting the values m = 1, n = 2 we obtain T short ≈ 0.13 and T long ≈ 911 which are in close agreement with the periods shown in Fig. 7. Displayed in Fig. 8 are time variations of the displacement of the corner of the plate using different grid spacings and time steps. As suggested by the expressions for 0 and T long the long period of the beat pattern decreases by a factor of four when the grid spacing is doubled. For example, with h = 0.02, T long ≈ 228 while with h = 0.04, T long ≈ 57 . On the other hand, the short period, T short , remains relatively the same as h varies. We also observe that as the grid gets coarser the time variation of the corner displacement begins to deviate slightly from the purely oscillatory beat profile shown in Fig. 7.
Having validated the numerical scheme we are now ready to present comparisons between the experiments and the simulations. In the experiment the forcing was applied by an oscillating cylindrical rod having a finite radius. This suggests using the following polar coordinate approximation for the Dirac delta function In the experiment the rod diameter was 12.7mm and the side length of the plate was 24cm. Thus, the ratio of the rod diameter to the side length is 0.0529, and the crosssectional area of the rod is negligible compared to the plate area. Because of this we have decided to model the forcing as a point source. Contrasted in Figs. 9,10,11,12,13, are nodal patterns for frequencies f = 190, 340, 490, 800 and 955Hz, respectively. The simulated nodal patterns shown are at a later time ( t = 40 ) when a steady pattern emerged. The agreement between the experiments and the simulations is quite good. For f = 190Hz the angular forcing frequency, 0 , needed to generate the simulated nodal pattern shown in Fig. 9 is close to that predicted by equation (23). Here, 0 = 13 was used in the simulation while equation (23) yields 0 = 12.4 . As the frequency increased the difference in 0 also increased. For example, for f = 955Hz equation (23) gives 0 = 62.1 while in the simulation 0 = 72 was needed to reproduce the nodal pattern in Fig. 13. Although a full explanation for this discrepancy is unknown, it is partly due to the experimental error in determining the resonance frequency f. To further validate the angular forcing frequencies, 0 , used in the simulations the following analysis was conducted. We first set u(x, y, t) = f (x, y) cos( 0 t) and substituting this into the governing eq. (8) yields Next, eq. (27) was discretized using second-order centraldifference approximations. Along and near the boundaries discretized versions of the boundary conditions (9)(10)(11) were used to relate the values of f at grid points lying outside the domain (i.e. ghost points) to those lying inside the domain. Doing this leads to an eigenvalue problem given by MF = F where = h 4 2 0 with h = Δx = Δy = 2∕N . Here, F is the column vector consisting of the (N + 1) 2 values f i,j = f (x i , y j ) in the following order and M is the (N + 1) 2 × (N + 1) 2 block pentadiagonal matrix having the form where B 1 , B 2 , B 3 , D 1 , D 2 , D 3 , I 1 , I 2 are (N + 1) × (N + 1) sparse matrices given in the Appendix. Thus, the eigenvalues of matrix M are related to the angular forcing frequencies. The MATLAB routine eig(M) was implemented to compute the eigenvalues. It was found that using N = 10 was sufficient in size to yield eigenvalues that are consistent with the angular forcing frequencies used in the simulations. For example, the angular forcing frequencies used in the simulations were 0 = 13(190Hz), 26(340Hz), 40(490Hz), 61(800Hz), 72(955Hz) . With N = 10 the following angular forcing frequencies were computed: 14.0, 26.5, 39.7, 62.6, 72.2. Of course, many more frequencies were obtained. Here, we list only those that are close to the values used in the simulations and the agreement is good. As N was increased the values of these frequencies changed but not by much, and more physically unrelevant frequencies emerged.

Summary
Investiagted in this paper was the free vibrations of a square plate that is forced at its center. The experiments were found to be in harmony with the numerical simulations that were conducted. A numerical solution procedure which is well suited to handle the complicated boundary conditions and forcing was presented; it consisted of a two-step finite volume scheme. The numerical scheme was checked by making comparisons with a known exact solution to a related problem. For the related problem the plate displacement at the corners displayed the formation of beats which were shown to be the result of a slight difference between the forcing frequency and the natural frequency of the discretized equation.

Conflict of Interest
The author wishes to report that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright