A Dynamic Coefficient Matrix Method for the Free Vibration of Thin Rectangular Isotropic Plates

The free flexural vibration of thin rectangular plates is revisited. A new, quasi-exact solution to the governing differential equation is formed by following a unique method of decomposing the governing equation into two beam-like expressions. Using the proposed quasi-exact solution, a Dynamic Coefficient Matrix (DCM) method is formed and used to investigate the free lateral vibration of a rectangular thin plate, subjected to various boundary conditions. Exploiting a special code written on MATLAB , the flexural natural frequencies of the plate are found by sweeping the frequency domain in search of specific frequencies that yield a zero determinant. Results are validated extensively both by the limited exact results available in the open literature and by numerical studies using ANSYS and in-house conventional FEM programs using both 12and 16-DOF plate elements. The accuracy of all methods for lateral free vibration analysis is assessed and critically examined through benchmark solutions. It is envisioned that the proposed quasi-exact solution and the DCMmethod will allow engineers to more conveniently investigate the vibration behaviour of two-dimensional structural components during the preliminary design stages, before a detailed design begins.


Introduction
Many vibrating airframe structural components could be modelled as thin plates.Not only that do these structural elements transmit various internal and external loads that may affect their stiffness but they are also frequently in close proximity to vibrating components such as engines.Therefore, it is of utmost importance to device and develop solution techniques to study the vibrational characteristics of these structures during preliminary design stages.Such vibrational analyses would allow the designers to investigate the effects of various boundary conditions the structural elements would be subjected to during its operation and the vibrational characteristics of the component before progressing to advanced stages of design.Using these results designers could alter the geometry or the materials used to avoid resonance and gain a favourable outcome.
Among the many methods available for vibration analysis, the analytical and semianalytical methods yield the highest accuracy but one major hurdle in using these methods is that they require the closed form solution to the governing partial differential equation.This can be a very tedious process if at all a tractable one.To circumvent this problem, many simplifying assumptions have been incorporated into the existing exact methods and as a result they exhibit many limitations.Having lost their generality, these exact methods are then only applicable to specific plate shapes, geometries, and those subjected to certain boundary conditions.
The orthogonality, completeness, and stability of Fourier series expansions have resulted in their frequent application to plate vibration problems [1].The Navier [2] and Levy methods [3,4] are two of the most common analytical procedures available for plate vibration analysis that incorporate such Fourier series expansions, where the former exploits a double Fourier series to solve the governing differential equation, the latter is based on a single Fourier series.However, both methods have a common drawback in that they are only applicable to plates having at least two simply supported boundaries.In addition, the Levy method is also limited to rectangular shaped plate configurations and is incapable of taking into account the effects of bending-twisting coupling.In addition to the above weaknesses, all methods that are based on conventional Fourier series expansions consist of a convergence problem along the boundaries arising as a result of discontinuities in displacement and its derivatives [1].Therefore, both of these methods are unsuitable for most aerospace applications as they could only tackle simple and special cases.In order to overcome the discontinuity in displacement and its derivatives along the boundaries, the Improved Fourier Series Method (IFSM) [5] was later proposed.Although IFSM possesses a higher rate of convergence and is more readily applicable to a host of plate configurations and boundary types, it is still inadequate to study problems comprising material and geometric nonlinearity.
The Rayleigh-Ritz method is another very popular exact method that has been exploited by many researchers in the past.It was first introduced by Rayleigh [6] and later improved by Ritz [7] by assuming a set of admissible trial functions, each of which had independent amplitude coefficients; thus, it is termed the Rayleigh-Ritz method or Ritz method.Young [8] and Warburton [9] used the Ritz method to study the vibration behaviour of rectangular plates.Later, Vijayakumar and Ramaiah [10] studied the vibration of clamped square plates using the Rayleigh-Ritz method.The flexural vibration of simply supported rectangular plates was investigated by Dickinson [11,12] using Rayleigh's method.One of the most comprehensive studies on thin isotropic rectangular plate vibration was carried out by Leissa [13,14] using the Rayleigh-Ritz method.Warburton [15] later extended the Rayleigh-Ritz method for the response calculation of a damped rectangular plate.The vibration of rectangular plates with elastically restrained edges was studied by Warburton and Edney [16].The Rayleigh-Ritz method was again used to study the vibration of rectangular plates using plate characteristic equations as shape function by Rajalingam et al. [17].However, the Ritz method in general is based on the weak form of the governing equations and is only applicable to self-adjoint problems.Furthermore, the choice of test functions in formulating the weak form is restricted to the approximation functions and it is required that the test and approximation functions are defined across the full domain of the problem, which is a major disadvantage.
The Galerkin method is also an analytical method which falls under the category of indirect classical variational methods.The Galerkin method has also been extensively used by researchers around the world.Although being somewhat similar in nature to the Rayleigh-Ritz method and belonging to the wide group of weighted residual methods, there are some distinct differences between the two techniques.Unlike the Rayleigh-Ritz method the Galerkin method commences with the weighted integral equations that are not comprised of boundary conditions.Thus, comparatively, the Galerkin method demands higher order approximation functions.Secondly, the Galerkin method does not require the system to be self-adjoint.But both methods take the test and approximation functions to be equivalent.Among many who exploited the Galerkin method for plate vibration analysis purposes, the transverse vibration of a rectangular plate was studies by Galin [18].Munakata [19] used the Galerkin method to investigate the vibration and elastic stability of a rectangular plate clamped at its four edges.Aynola [20] and Stanisic [21] also studied the vibration behaviour of rectangular plates using the Galerkin method.Laura and Saffell [22] investigated the small-amplitude vibration of clamped rectangular plates.Later Laura and Duran [23] applied the Galerkin method to determine the vibration characteristics of a clamped rectangular plate subjected to forced vibration.Nevertheless, one of the biggest drawbacks associated with classical variation methods in general such as Rayleigh-Ritz and Galerkin methods is the difficulty involved in accurately developing the approximating functions for arbitrary domains.This difficulty associated with constructing the arbitrary test and approximate functions that should satisfy essential edge conditions, smoothness levels, linear independence, and completeness and continuity conditions is a massive limiting factor and the complicatedness of the problem becomes even more severe in magnitude for difficult geometries commonly found in most aerospace structures.Therefore, the lack of a credible method to formulate proper approximation functions for a specific geometry drastically reduces the convergence quality and applicability of classical variation methods.
The method of superposition is also a very powerful approximate analytical method that has been used extensively by many researchers in the past to obtain highly accurate results for problems involving plate vibrations.It was developed by Gorman [24] who utilised it to analyse the vibrational behaviour of thin isotropic rectangular plates.In this method, the plate is divided into a number of subsystems, termed building blocks, under different boundary conditions compared to the global system, and subjected to a distributed force, moment, rotation, and translation [24].The steadystate response of each subsystem is then superimposed.Unlike most other exact methods, this method is applicable to a variety of plate types, which include orthotropic, hybrid, and laminated plates.The superposition technique also allows for the application of various classical and nonclassical boundary conditions as well as loading configurations and is readily applicable to thin plates, thick Mindlin plates, transverse shear deformable laminated plates, and open cylindrical shells.Furthermore, throughout the entire domain of the plate, the governing differential equations are satisfied exactly by all the solutions [24].Gorman and Sharma [25] used the superposition method to conduct a free vibration analysis of rectangular plates.A free vibration analysis of cantilevered plates was also carried out by Gorman [26] using the superposition method.Later, Gorman [27] also conducted a study on the free vibration analysis of completely free rectangular plates using the superposition-Galerkin method.However, the main problem with the method of superposition is that, for mixed boundary types, it has not been verified yet if the results yielded are an upper bound or a lower bound.Thus, this uncertainty may well be a problem when trying to estimate the error of the results.
Among the exact methods commonly used for the vibration analysis of plates is the dynamic stiffness method (DSM), which was first presented by Kolousek [28] in the forties.Later Boscolo and Banerjee [29] applied DSM to study the vibration behaviour of plates using both the Classical Plate Theory (CPT) and the First-Order Shear Deformation Theory (FSDT).Banerjee and Papkov [30] also presented a DSM solution of a rectangular plate for the most general case.Subsequently, the free vibration of plates subjected to arbitrary boundary conditions was investigated by Liu and Banerjee [31] using a novel spectral dynamic stiffness method.However, the DSM method is cumbersome to use when applied to complex, real-life plate configurations consisting of material and geometric nonlinearity.
Thus, the objective of this work is twofold.Firstly, the authors wish to develop a new quasi-exact solution to the plate governing equation by treating the governing equation as a sum of two beam-like expressions, an approach that does not incorporate any simplifying assumptions, thus, preserving the generality of the solution and which, to the best of the authors' knowledge, has not been explored before.
The second objective will be to develop a new Dynamic Coefficient Matrix (DCM) method for the modal analyses of thin rectangular plates, having any aspect ratio, based on the new quasi-exact solution.To the best of the authors' knowledge, the new DCM method built upon a quasi-exact Dynamic Coefficient Matrix has also not been developed and presented in the open literature.What distinguishes the DCM method from other classical exact methods is the frequencydependent nature of the resulting system's matrix and most importantly the fact that its generality is not compensated by any simplifying assumptions.Together, the new DCM method and the quasi-exact solution would, upon further development in the future, provide researchers with the flexibility to study the vibration of thin rectangular plates of any dimension or thin isotropic plate assemblies modelled using rectangular elements, subjected to any boundary condition.

Theoretical Background
Consider a linearly elastic, homogeneous, isotropic, thin rectangular plate, as shown in Figure 1, having length , width , and thickness ℎ.The thickness ℎ is assumed to be much smaller compared to the other characteristic dimensions as well as the wavelength.Thus, Classical Plate Theory is used for the purpose of this study.As a result, during vibration only small deflections are assumed and the rotary inertia and shear effects are neglected.
The governing partial differential equation for the plate [32] will take the following form: where  is the flexural displacement,  is the density,  is the time, and  represents the plate modulus and it is defined as follows: In order to obtain a quasi-exact solution to (1), the roots have to be determined.To this end, a new approach is taken here, which is to decompose the plate equation into two separate beam-like expressions representing each spatial coordinate direction of the plate.The main steps of this procedure are outlined below.If the solution is assumed to take the following form, then the characteristic equation will be as follows: This could be rewritten as follows: where  1 and  2 are the mass distribution constants along the and -directions, respectively.These constants were introduced to decompose the plate governing equation in to the two beam-like expressions above.Through careful observation it could be seen that simply plugging numerical values in place of these constants will allow one to reconstruct and rewrite the plate governing equation into its original form.The numerical values of  1 and  2 can be anything between 0 and 1 (0 <  1 and  2 < 1); however, the sum of the two mass distribution constants should be unity ( 1 +  2 = 1).For example,  1 and  2 will both be equal to 0.5 for a square plate.They will assume other values for other rectangular plate shapes.The term ( * ) represents the -direction and the term ( * * ) is for the -direction of the plate.In both expressions,  is the coordinate in the x-direction and  is the coordinate in the -direction.The terms ( * ) and ( * * ) are treated as two different beam equations for determining roots.Furthermore, in expression ( * ),  can vary and  is held constant and, for the term ( * * ),  is held constant and  is allowed to vary.The quadratic formula was then applied on the expressions ( * ) and ( * * ) separately.Simplification resulted in the following roots for expression ( * ) of the plate governing equation.
An identical procedure resulted in the roots for the term ( * * ) in ( 5) above.
The detailed mathematical manipulations for both simplification processes are not included here for brevity.Thus, from (6) it could be seen that there will be four roots,   , ( = 1, 2, 3, 4) for the expression ( * ), as defined in ( 8) and ( 9), of which two are real and two are imaginary.
Similarly, from (7) it could be found that there are four roots,   , ( = 1, 2, 3, 4) for the expression ( * * ), of which two are real and two are imaginary, as defined in the following: It is important to note here that the not only do roots shown in ( 8) through (11) satisfy their individual expressions separately, but together any real-real or imaginary-imaginary combination (  ,   ) of these roots also satisfy (5) as a whole.Thus, each real-real and imaginary-imaginary pair of roots (  ,   ) is an exact solution to the plate governing equation.
As the solution to the plate equation was assumed to take the form shown in (3), the following expressions were constructed using the roots shown in ( 8) through (11).
=  1 sin (  ) +  2 cos (  ) +  3 sinh (  ) +  4 cosh (  ) (13) where  1 to  4 and  1 to  4 are unknown coefficients.Since the solution is assumed to take the form defined by (3) the final 16-term quasi-exact solution for a thin plate could be derived by multiplying (12)  where,   , in (14), are the new unknown coefficients defined as follows: Thus, the non-nodal flexural displacement (, ) anywhere in the plate could be written in the matrix form as follows: where the row vector ⟨⟩ is the solution vector which contains the roots to the plate governing differential equation and the column vector {} is the vector of unknown coefficients.
The slope along the -direction could then be written as follows.
In ( 17) above the row vector ⟨  ⟩ is determined by differentiating the solution vector ⟨⟩ with respect to .Similarly, the slope along the -direction could be expressed as follows: where the row vector ⟨  ⟩ is obtained by taking the derivatives of the roots  1 to  16 contained within the with solution vector, with respect to .The curvature of the plate   could also be represented as follows: where the row vector ⟨  ⟩ is determined by obtaining the derivatives of the solution vector with respect to both  and .
The boundary conditions for the displacements are as follows: By applying the boundary conditions for displacements, i.e., substituting (20) into ( 16) through (19), the following matrix relationship is obtained.
The expression in (21) above could be written in the shorthand form as fp;;pws: where [  ()] is the symmetric 16 x 16 Dynamic Coefficient Matrix (DCM) of the system.The stiffness matrix in (22) consists of the essential requirements to compute the natural frequencies for a thin rectangular plate subjected to any boundary condition.To obtain the natural frequencies using the Dynamic Coefficient Matrix (DCM) method, boundary conditions are applied on the dynamic stiffness matrix and a determinant sweep is conducted by sweeping the frequency domain in search of frequencies at which the determinant of the DCM will be zero; i.e., |  ()| = 0.

Results
Numerical checks were performed to confirm the predictability, accuracy, and practical applicability of the proposed Dynamic Coefficient Method (DCM) method, programmed in a MATLAB5 code.In what follows, an illustrative example of homogeneous, rectangular, thin plate is examined.At first, the natural frequencies for the plate with one edge clamped and other three edges free (C-F-F-F) were investigated, where the exact results from reference [13], together with the frequency data obtained using ANSYS5 and in-house conventional FEM programs, based on both 12and 16-DOF plate elements, were used as the benchmarks for comparison and to validate the DCM solution method.For further studies, 10 more different sets of boundary conditions were considered, where the DCM results were validated against only exact results from reference [13].
The aspect ratio, in this case, is (/) = 1.5.However, as explained before, the DCM formulation can be applied to any thin rectangular plate configuration with any aspect ratio.In what follows, the natural frequencies of such a plate, determined using the new DCM method outlined in Section 2, are presented for various sets of boundary conditions.The notation S-F-S-F, for example, will identify a rectangular plate whose edges  = 0,  = 0,  = , and  =  are subjected to pinned, free, pinned, and free boundary conditions, respectively.The results of these modal analyses are included below.
In Table 1, the DCM plate natural frequencies for one edge clamped and other three edges free (C-F-F-F) boundary conditions are presented alongside and are compared with the exact data [13] and those obtained using various conventional FEM formulations.As can be seen, the first five natural frequencies obtained from DCM are in perfect match with the exact values reported in [13]; i.e., zero relative error.The ANSYS5 results, obtained from a 196-element mesh model, show slight differences with DCM/exact data, with the maximum error of less than 1% for the highest mode.These slight discrepancies can be attributed to the fact that while using ANSYS5, the 3D, 4-noded, SHELL-181 element was used to model the system.This element is a shell element, which has 6 DOFs per node and these are the three translations along, and three rotations about the x-, y-, and z-axes.
The comparison is also made between the DCM results and those obtained using the in-house FEM code, where again 196-element mesh models of 12-and 16-DOF conventional FEM plate elements are used.In general, when compared with the DCM/exact frequencies, the results obtained from the 16-DOF plate elements show the lowest differences (a max of 0.77%), followed by ANSYS (a max 0.99%) and those evaluated using 12-DOF plate element (a max of 1.13%), respectively.
As can be seen from Tables 2-11, the presented DCM method produces exact results for the first five natural frequencies of a thin rectangular plate, subjected to any type of boundary conditions.Thus, the accuracy of the unique solution procedure adopted in determining the roots of the plate governing equation and subsequently forming the quasi-exact solution is validated.Unlike most exact methods available, which are limited to certain configurations and special boundary conditions, the DCM method presented here is a powerful tool that can be used to study the vibration behaviour of square or rectangular thin plates of any dimension and subjected to any set of boundary conditions.

Conclusion
A new, quasi-exact, frequency-dependent solution was developed for the free flexural vibration of thin (Kirchhoff) rectangular plates using a distinctive procedure of splitting the thin plate governing equation in to two beam-like expressions.Using these quasi-exact solutions to the governing equation the Dynamic Coefficient Matrix (DCM) of the thin plate was developed.The boundary conditions of the system were applied using a special code written in MATLAB5 and the natural frequencies of a thin plate subjected to various sets of boundary conditions were investigated to validate the accuracy of the new quasi-exact solution and the DCM method.When investigating the system's first five natural frequencies, the results were found to match perfectly with exact results from the open literature.Further research is being carried out to extend the DCM method, based on the new quasi-exact solutions, to thick and multilayered plates.

Table 1 :
Natural frequencies for a plate with one edge clamped and three edges free (C-F-F-F).

Table 2 :
Natural frequencies for fully pinned plate (S-S-S-S).

Table 3 :
Natural frequencies for two opposite edges pinned and two edges clamped plate (S-C-S-C).

Table 4 :
Natural frequencies for three edges pinned and one edge free plate (S-S-S-F).

Table 5 :
Natural frequencies for two opposite edges pinned and two edges free plate (S-F-S-F).

Table 6 :
Natural frequencies for all edges clamped plate (C-C-C-C).

Table 7 :
Natural frequencies for three edges clamped and one edge free plate (C-C-C-F).

Table 8 :
Natural frequencies for two adjacent edges clamped and two edges pinned plate (C-C-S-S).

Table 9 :
Natural frequencies for two adjacent edges clamped and two edges free plate (C-C-F-F).

Table 10 :
Natural frequencies for two opposite edges clamped and two edges free plate (C-F-C-F).

Table 11 :
Natural frequencies for three edges clamped and one edge pinned plate (C-C-C-S).