A Framework for Extension of Dynamic Finite Element Formulation to Flexural Vibration Analysis of Thin Plates

Dynamic Finite Element formulation is a powerful technique that combines the accuracy of the exact analysis withwide applicability of the finite element method. The infinite dimensionality of the exact solution space of plate equation has been a major challenge for development of such elements for the dynamic analysis of flexible two-dimensional structures. In this research, a framework for such extension based on subset solutions is proposed. An example element is then developed and implemented in MATLAB software for numerical testing, verification, and validation purposes. Although the presented formulation is not exact, the element exhibits good convergence characteristics and can be further enriched using the proposed framework.


Introduction
Dynamic stiffness modeling is a well-established technique in vibrational analysis of structural elements.These methods seek to propose formulations that have the accuracy of the exact solutions and wider applicability of the finite element methods by incorporating some form of the closed form solutions of governing equations instead of polynomials used by classic finite element method (FEM).
Two of the most famous dynamic stiffness formulations, mainly applied to various beam-structures, are dynamic stiffness matrix (DSM) and Dynamic Finite Element (DFE) methods, which produce accurate results with much coarser mesh compared to FEM formulations.
Extensively developed in 1970s, and pioneered again by Banerjee and his coworkers since 1990s [1][2][3][4][5][6][7][8], the beam DSM formulation dates back to much earlier times [9].The 1941 work by Kolousek [9] is probably the first to derive the dynamic stiffness matrix for the Euler-Bernoulli beam.The DSM formulation is performed by developing a dynamic relation between force and displacement matrix for a problem domain, where the displacement functions are combination of trigonometric and hyperbolic frequencydependent functions, which satisfy the Euler-Bernoulli or Timoshenko beam equations.In this method, the boundary and loading conditions of the problem in hand are taken into consideration in developing the force-displacement relationship and thus this process yields a one-step formulation to the problem but is very case specific [1][2][3][4][5][6].Using truncated Taylor series expansion and assigning penalty for different boundary conditions, the DSM has been more recently generalized by Pagani and his coworkers [7,8] to vibrational analysis of composite beams and solid and thin-walled structures subjected to various boundary conditions.This approach was also further extended to use higher-order 1D dynamic stiffness elements for the vibration analysis of composite plates [10].
On the other hand, the DFE technique follows the general procedure of FEM element development by using the general displacement solution of beam equation as basis functions.Next, the element dynamic shape functions are developed and supplied to a residual minimization scheme, such as Galerkin's technique, to develop the frequency-dependent element matrix.The matrix thus developed is not case specific and can be assembled.The boundary and loading conditions are applied in a manner similar to FEM [11][12][13][14].
Both DSM and DFE methods generate one matrix per problem which has the effects of mass and stiffness combined.In fact, the mass terms are distributed over the element displacement as a dynamic force using D' Alembert's principle.The assembly of the element dynamic stiffness matrices and the application of the system boundary conditions then lead to a nonlinear eigenvalue problem, which is then solved using a root finding technique, such as Wittrick-Williams algorithm [15].
Both the DSM and DFE formulations have been very successful for one-dimensional elements such as rods, beams, and beam-structures under various loading and boundary conditions.The expansion of such techniques to twodimensional (2D) elements such as plates, however, has led to cumbersome mathematics with limited applicability which served as a motive for current research.The particular problems with such extension arise from the fact that the plate governing equations are two-dimensional partial differential equations with infinite number of solutions.This is unlike beam and rod situation where the frequency-dependent solution of the ordinary differential equation can be used directly as the basis functions of approximation space to develop set of dynamic shape functions [16].
The other problem in the case of DFE formulation comes from evaluating the resulting integral equations, written in terms of the dynamic shape functions, to develop element matrices.Since these functions are transcendentally dependent on frequency, the integrations are difficult to handle given the current computational power of computers especially for infinite dimensional domains.This difficulty is particularly important when extension of such elements to irregular shapes is required.In such situations, the area integrals used in mapping of the element to its natural coordinate system will involve the absolute value of Jacobian mapping matrix that is difficult to handle for such complicated displacement functions.
Because the Dynamic Finite Element matrices are a result of these integrations, use of numerical integration schemes such as Gauss quadrature formula will require introduction of numerous evaluation points for an acceptable result.This is unlike FEM extension to plates of arbitrary shapes where the polynomial shape functions provide descent result with small number of evaluation points.Besides, FEM shape functions are much easier to evaluate numerically as they do not have any dependency on vibrational frequency and produce constant matrices of numerical values that can be manipulated with ease.Therefore, much of the available literature for exact solution of plates was focused on solving these equations for simple geometries and loading cases, as presented by Leissa [17][18][19][20][21][22], or involved variation error minimization techniques, such as finite element method [23][24][25][26][27][28][29][30][31][32][33][34].
Recently, there have been more researches to extend the dynamic stiffness matrix scheme to plate vibrations.Since this method involves force-displacement relationship development to arrive at the element DSM, it requires developing differential relationships, avoiding the integration problem.Casimir et al. [16] developed the DSM for rectangular thin plates and demonstrated that accurate results can be achieved with limited number of elements.However, their method involved solving infinite dimensional matrices, limited to simple cases, and provided problem specific formulation that must be reformulated for each new configuration.This reformulation for plate elements would therefore involve solving large matrices and is less favourable to be used as a general analysis tool.More recently, Liu and Banerjee [35][36][37] presented a novel Spectral-DSM (S-DSM) method to generalize the DSM using truncation of Fourier series and the application of spectral analysis.The S-DSM, when applied to the free vibration analysis of complex problems involving isotropic laminated plates under various boundary conditions, demonstrated significantly improved efficiency over the conventional FEM formulations.
In this paper, the focus is to present an alternate dynamic stiffness formulation in the form of DFE for the vibration analysis of Kirchhoff plates.Particularly, the focus is to generate a model that can be easily extended to plates of different shapes without loss of accuracy of the original formulation.The mathematical procedure of beam DFE formulation served as a guideline for development of the current plate DFE design, as presented in the following section.An example 4-node 12-DOF element is then developed and validated against FEM, where its accuracy and efficiency are demonstrated through the free flexural vibration analysis.

Mathematical Modeling
To develop Dynamic Finite Element of thin plates, we begin with the partial differential equation governing the free flexural vibration [38]: where  is the flexural displacement,  is the mass density,  is the thickness, and represents bending coefficient.The terms  and ] denote Young's modulus and Poisson's ratio of the plate material, respectively.Applying Galerkin's weighted residual minimization scheme [39], the integral form of the governing equation ( 1) is written as where (, ) is the residual weighting function and  is the differential of the area within closed element boundary.
Green's theorem is a conversion mechanism between area and line integrals performed over a closed boundary.For functions  and  Green's theorem states here,   and   are components of the unit normal vector  to the element boundary, and  is differential segment along the boundary.Since Galerkin's method uses the same approximation for both field variable  and weighting function , by applying Green's theorem (4) to the above area integrals twice and taking the conformity of virtual displacements and loads across element boundaries into consideration, one can change Galerkin's integral equation ( 3) to its equivalent weak form: where (, ) = [(, )][  ] can also be interpreted as the virtual displacement.Equation ( 5) is the formulation generally used for finite element development of plates.The DFE derivation will deviate from FEM at this point.
Let (, , ) = [  (, )][  ] be the dynamic flexural displacement, where [  (, )] is a row vector containing the dynamic shape function, derived from the frequencydependent (dynamic) basis functions of approximation space, and [  ] is the vector of nodal displacements.The dynamic basis functions, in this case, are defined to be the solutions of the Kirchhoff 's plate governing equation.
Applying Green's theorem (4) to the above area integrals twice, ( 5) can be written as Note that the last 4 area integrals can be combined to generate Kirchhoff 's plate equation in terms of virtual displacement  and thus can be written as Assuming, as mentioned earlier, the dynamic shape functions are derived from the dynamic basis function, which are in turn solutions of Kirchhoff 's equation, the terms in the bracket [ 4 [  ]/ 4 + 2( 4 [  ]/ 2  2 ) +  4 [  ]/ 4 − ( 2 /)[  ]] will yield zero.Therefore, the DFE stiffness matrix [  ] can be obtained using only evaluation of line integrals performed over element boundaries, as follows: Although not explicitly visible, the above matrix is symmetric as it was obtained using mathematical manipulation of the symmetric Galerkin's weak integral form.

Extension to Arbitrary Shapes.
The absence of area integrals in DFE stiffness matrix evaluation enables exact integration of the boundary line integrals when transformed between different coordinate systems using variable substitution.Consider, for example, the quadrilateral case given in Figure 1(b).
The variables  and  from element coordinate system can be related to the variables  and  from natural coordinate system using bilinear relations: where the coefficients   and   can be obtained by evaluating these relations for the nodes of the element.For   , The DFE stiffness matrix (8) defined in last section was developed for any arbitrary closed area in element domain.Using relationships between element and natural coordinate systems (9) developed above, differentials of DFE matrix can be replaced with appropriate natural coordinate corresponding terms.The first integral, for example, can be evaluated in the natural coordinate system as In the above integrations, the shape function vector, [  ], and its element coordinate system (Figure 1(b)) derivatives must be converted to the natural coordinate system (Figure 1(a)) as well before the integration is performed.This can be done by evaluating the derivatives in the element coordinate system and then replacing the variables  and  with the respective definitions in natural coordinate system.Alternatively, one can transfer the shape function to natural coordinate system and then define the relation between derivatives in natural coordinate system and element coordinate system using application of chain rule.It is important to note that the relationship between the natural and element coordinate system does not need to be bilinear and the geometry is not restricted to quadrilaterals.The same procedure is applicable to any geometry and any relationship.

Development of Example Element
The governing differential equation of Kirchhoff plates is a two-dimensional homogenous partial differential equation.The solution to such equations forms an infinite dimensional Hilbert space, which is usually prescribed in terms of Fourier series for dynamic formulations of thin plate vibrations.Although the DFE mathematical formulation presented does not have any restriction on the dimension of the shape function vector used, to generate a widely applicable scheme, a subset of the infinite solution space is used to develop example element.In this approach, instead of using the general Fourier trigonometric functions, a solution space is generated considering the physics of the element equation rather than a pure mathematical approach.Since in the development of the DFE the entire solution space is not used, the solution will not be exact for one element but would rather converge to an analytical answer.Special focus will be given to solutions that have both the following characteristics: (1) Being symmetric with respect to  and  (2) Being  ∞ continuous.
Referring to the differential equation of motion (1), governing a harmonically oscillating Kirchhoff (thin) plate, while there are multiple solutions possible to this equation, one simple subset is obtained by considering the following reduction: Note that this simplification has made the plate equation similar to that of a thin beam in  direction.For this simplification, it is necessary that the solution considered must make the summation  4 / 4 + 2( 4 / 2  2 ) equal to zero.One such solution set can be obtained when each of the terms of the summation is zero individually.These solutions can be found using separation of variables and solving the equation above as an ordinary differential equation.Assigning The first set is identical to beam differential equation solutions in  direction.Similar 8 other solutions can be developed for a reduction in  direction, written as cos (  ) , sin (  ) , cosh (  ) , sinh (  )  cos (  ) ,  sin (  ) ,  cosh (  ) ,  sinh (  ) .

(15)
Now, the FEM development can be followed to produce the shape functions using these solutions as basis function.It is worth noting that although the above 16 basis functions developed here are linearly independent in a continuous Hilbert space, when evaluated over some finite nodal points, they are found to produce linearly dependent values.This problem is resolved by mixing these solutions to create 12 new basis shape functions which in low frequency values (i.e., when the frequency tends to zero; static behavior) will converge to 12 Hermitian polynomial basis functions commonly used in the conventional FEM formulation of thin plates [38], that is, 1, , ,  2 , ,  2 ,  3 ,  2 ,  2 ,  3 ,  3 ,  3 .
The new dynamic basis functions are written as The nonnodal displacement function   can be defined as a linear combination of these shape functions, written as By evaluating these displacements and corresponding derivatives at nodes of an element (Figure 1(a)), the following relations are obtained: The constant matrix [] can be obtained as where for a 2 by 2 square element, [ 12×12 ] matrix is obtained as Shock and Vibration When premultiplied by the basis function matrix [  ], the inverse of matrix [ 12×12 ] generates the shape function vector [  (, )], containing the 12 new dynamic (frequency-dependent) shape/interpolation functions: Note that although the shape functions obtained through this method are variables of the rotational vibration frequency  on the element domain, they still maintain the characteristics of Hermitian shape functions at nodes in that they produce a value of 1 at their corresponding nodes and zero at other nodes.Also, these shape functions converge to their FEM counterparts for small values of vibrational frequency, when  tends to zero.The expressions for dynamic shape functions are very lengthy and are omitted here for brevity.For demonstration purposes, the first resulting dynamic interpolation function is presented graphically in Figure 2, for two values of   = 10 − 6 (approaching zero, Figure 2 As can be seen from Figures 2(a) and 2(c), as the frequency, , approaches zero, the dynamic shape function,   (, ), becomes the same as the Hermitian shape function, 1 [40], that is, Using the displacement and weighting functions ( 16) expressed in terms of these dynamic shape functions and introduced in the DFE boundary integral equation ( 8), the DFE frequency-dependent element dynamic stiffness matrix, [()], is obtained.
Having this matrix, next an illustrative example of a square plate with three sides simply supported and one side free boundary condition is analyzed.The natural frequency is normalized over plate material and mechanical properties using the relation √/.The element length and width are equal to 2 units and Poisson's ratio is taken as 0.3.A conventional FEM analysis using Hermitian shape functions [40] is also performed to highlight the improvements obtained from DFE formulation.
The element matrices are assembled, where applicable, to find the system's global dynamic stiffness matrix [()].The system boundary conditions are then applied by modifying the DFE matrix, where the constrained degrees of freedom (DOF) are enforced by eliminating the corresponding rows and columns, in a similar way to the conventional FEM.This produces a nonlinear eigenvalue problem, written as where {} is the vector of system's DOFs.The system natural frequencies are then evaluated by setting the determinant of DFE matrix equal to zero, |()| = 0.The determinant of [()] is a transcendental function of natural frequency of vibration that must then be swept across the frequency domain to find the values for which determinant function is zero.For a plate with 3 sides simply supported and one side free (SS-SS-SS-F) boundary condition, modeled using a single DFE, the determinant matrix is plotted in Figure 3.The determinant function also produces undefined regions, known as poles [41], where the determinant value approaches infinity.These regions represent the denominator of the determinant function [12].
The plate's first five nondimensionalized natural frequencies (√/) were evaluated using FEM and DFE techniques and compared with the analytical data from [17] (see Table 1).
As demonstrated in Table 1, for a single element mesh, the performance of the presented DFE plate surpasses that of its FEM counterpart, especially for the fundamental frequency where DFE produces very small error.It is also worth noting that Hermitian FEM vibrational analysis produces a linear eigenvalue problem in form of a polynomial in terms of frequency.This polynomial, in turn, has finite number of roots (in the case of above boundary condition, two roots only) and therefore, detection of higher modes requires mesh refinement.The DFE element, however, leads to a nonlinear eigenproblem through which one can extract multiple frequencies, even when a single-element mesh is used.The 2 × 2 square plate is further analyzed using both FEM and DFE formulations and under three different boundary conditions, for which analytical results were reported by Leissa [17].These boundary conditions include Clamped-Clamped-Free-Free (C-C-F-F) boundary condition, Clamped-Free-Simply Supported-Free (C-F-SS-F) boundary condition, and Simply Supported-Free-Simply Supported-Free (SS-F-SS-F) boundary condition.Since FEM can capture different number of frequencies for different boundary conditions, only the fundamental natural frequencies are reported here for consistency in comparison.The result of this analysis is reported in Table 2.
Besides the performance of a single DFE element, it is important to analyze the convergence characteristics of the developed DFE plate.The convergence of the DFE element was compared against the 12-DOF plate finite element, which utilizes Hermitian polynomials.For this study, the mesh density was increased from a single element to 4 elements (2 by 2), and to 25 elements (5 by 5).For illustration purposes only, the result for the second and fifth mode of vibration of a plate with (SS-SS-SS-F) boundary condition is compared in Figure 4 for DFE and FEM formulations as examples of monotonic and nonmonotonic convergence to the analytical solution reported in [17].Similar convergence was observed for other modes as well.
As can be seen from Figure 4, the DFE plate has comparable convergence rate with FEM.However, as mentioned earlier, since DFE produces a transcendental determinant equation, it is capable of providing potentially infinite number of natural frequencies, while FEM plates create constant mass and stiffness matrices which can only generate a finite number of solutions through linear eigenvalue problem.Also, the plate DFE presented here was produced based on one subset of solution of the plate equation, which reduced the plate behavior to two beam-like equations.Many other solution sets can be found, for example, based on distribution of dynamic load over the element strain, and combined with these solutions to enrich element physics further.The technique can also be extended to generate plate elements with higher number of nodes and DOFs.

Conclusion
The well-established framework for development of beam Dynamic Finite Element (DFE) formulation, based on exact solutions of governing differential equation, is extended to two-dimensional (2D) plate elements.Unlike other dynamic plate formulations proposed in past, which exploit an infinite dimensional solution, an element with limited number of solutions was developed.The solutions obtained were founded on physics-based simplification of Kirchhoff 's equation to beam-like equations and demonstrated comparable performance with classic FEM formulations.The research is underway to further enrich the presented 2D element by using the proposed framework to potentially achieve a quasiexact plate Dynamic Finite Element (DFE) formulation and further extend DFE plates to arbitrary geometries of various distortions.

Figure 1 :
Figure 1: Natural coordinate system (a) is used to map arbitrary shapes from element coordinate system (b) on a 2 by 2 square.

Figure 2 :Figure 3 :
Figure 2: Variation of dynamic shape functions with vibrational frequency and comparison to Hermitian shape functions.The plot of first DFE shape function,  1 (, ), for   = 0.000001 (a) and   = 2 (b).The first Hermitian shape function, 1, is plotted in (c).

Figure 4 :
Figure 4: Comparison of convergence of √/ values to analytical solution provided by [20] between FEM and DFE models for second mode (a) and fifth mode (b) of vibration of a plate under SS-SS-SS-F boundary condition.

Table 1 :
Comparison of DFE and FEM formulations under SS-SS-SS-F boundary condition against analytical solution for values of √/.

Table 2 :
Comparison between DFE and FEM formulations against analytical solution values of √/ for fundamental natural frequency under different boundary conditions.