A New Shear Locking-Free Quadrilateral Plate Element for Free Vibration and Eigen-Buckling Analysis

A new plate element has been formulated by simply adding one degree of freedom for each node on the basis of the classical ACM element, based on a modified first-order shear deformation plate theory, which assumes that the plate deflection has bending and shearing compo-nents, with two rotations defined by differentiations of bending deflections. This new element, with 16 degrees of freedom, exhibits no locking effect, and is appropriate for both thin and thick plates. To assess the performance of the newly developed plate element, it has been applied to the prediction of natural frequencies and critical bucking loads of square plates with various thickness ratios and boundary conditions. Comparisons with available analytical and numerical solutions in open publication show that the new element leads to highly accurate solutions in free vibration and eigen-buckling analysis with good convergence, and can be competitive with existing plate elements.


Introduction
In the last few decades, extensive efforts have been devoted to numerical schemes for flexural vibration and eigen-buckling of plates, which are widely used in engineering structures. For the eigenvalue problem of a rectangular geometry, exact or analytical solutions, which often serves as benchmarks for validation of numerical methods, do not always exist except for the special case where at least two parallel edges are simply supported. As for the other boundary conditions, many numerical solutions have been generated by approximate approaches, such as the Rayleigh-Ritz method [1,2], the differential quadrature method [3,4], the finite strip [5] and finite element methods [6,7], etc., among which the finite element method has proven to be very popular, and remains as a focused topic in the research area of computational mechanics [8][9][10][11][12] with growing publication.
It should be pointed out that in early days, most of the research concentrated on the classical plate theory (CPT) which requires displacement and its derivatives to be continuous across element boundaries, thus C 1 continuity is a necessity. However, the neglect of rotary inertia and transverse shear deformation may lead to overestimation of natural frequencies for moderately thick plates.
Efforts, therefore, have been devoted to developing elements based on the Mindlin plate theory (MPT) [13,14] which accounts for shear deformations, and requires only C 0 continuity for the displacement field. Comparative studies [15][16][17] have shown that MPT can predict natural frequencies highly consistent with three-dimensional elastic solutions for a variety of plate geometries. But unfortunately, for the finite element approach, formulating a plate bending element based on MPT faces the inevitable difficulty known as the shear locking effect, which arises due to the contradiction between the independence and non-independence among the transverse deflection and two rotations, leading to infinite shear rigidity and unreliable numerical solutions when thick plates are degenerated to thin ones.
In seeking a remedy to the shear locking effect, many researchers have proposed a large number of procedures, among which the reduced and selective reduced integration scheme [18][19][20][21] is the most widely used one to correct the locking effect. However, this approach cannot completely eliminate shear locking and it may cause singularity of the stiffness matrix which will often lead to spurious zero-energy modes. Utilizing a mixed formulation based on the multi-variable Hellinger-Reissner principle [22][23][24][25] has also proven to be effective in alleviating shear locking, but its practical use might be limited by the difficulty of matching between various field variables and higher requirements on computing resources due to increased complexity of element formulation.
In recent years, continuous efforts have been made to tackle the aforementioned challenges. Nguyen et al. [26] used a piecewise-linear shape function to derive a polygonal element, in which the assumed strains along the element edges are interpolated based on the Timoshenko's beam formulas so that the shear locking phenomenon can be naturally avoided. Liu et al. [27,28] proposed the smoothed finite element method (SFEM) by integrating the strain smoothing technique into the classical FEM, and generated various smoothed FEM models which are not only free of shear locking but also insensitive to mesh distortion. Senjanovic et al. [29] gave a shear locking free finite element by decomposing the Mindlin mathematical model of three DOFs into a single-DOF bending model and a double-DOF shear model.
Motivated by the idea of including the shear deformation effect in available thin plate elements [30,31], the present research proposes a new shear locking-free plate bending element by introducing an extra nodal degree of freedom taking the shear deformation into account on the basis of the classical ACM element [32], following a modified first-order shear deformation theory [33] which only involves two independent variables, namely, the bending deflection wb and the shearing deflection ws.
To examine the use of this new developed element in calculating natural frequencies and buckling loads of plates, numerical examples considering various thickness ratios and a range of boundary conditions are presented, and a comprehensive comparison is carried out with available analytical and numerical results in open literature.

Modified shear deformation theory
It is widely acknowledged that, CPT can be viewed as an extension of the classical Bernoulli-Euler beam theory (CBT). Likewise, MPT is a two-dimensional equivalent of Timoshenko beam theory (TBT). Both CPT and CBT have one independent variable and one governing equation, however, a Timoshenko beam has two, while MPT has three. Based on this, a two-variable first-order shear deformation theory [33] was developed which decomposes the total deflection into the bending deflection wb and the shearing deflection ws, and defines two rotations by the first derivative of wb. Consider a rectangular plate (length a, width b, uniform thickness h) with its middle surface located in the x-y plane of a Cartesian coordinate system shown in Fig. 1.
The displacement field of this modified first-order plate theory is: where: u, v and w denote the displacements in the x-, y-and z-axes, respectively; wb and ws are the deflections caused by the bending and shear deformation, and rotations of a transverse normal with respect to the y-and x-axes are denoted by ∂wb/∂x and ∂wb/∂y, respectively.
These stains can be divided into two groups, one in terms of bending and the other of shearing: Similarly, the stress field can be written as: , , , , , where the bending and shear rigidity matrices Db and Ds are defined, respectively, as: (5) in which υ is the Poisson's ratio; E and G = E/2(1+υ) are Young's modulus and the shear modulus, while κ represents the shear correction factor. The governing differential equations of motion can be derived from the following Hamilton's principle to guarantee variationally consistency: where: Π, U, V and T denote the total potential energy, the strain energy, the potential energy due to external forces and the kinetic energy, respectively, and δ is the symbol of variation.
For free vibration problem, the kinetic and potential energies can be expressed in terms of wb and ws as: where: J = h 3 /12 is the axial moment of inertia, and ρ is the density of the plate. Taking wb and ws as two independent variables and utilizing the Hamilton's principle, the governing differential equations of motion for free vibrations of the plate can be obtained as: where: D = Eh 3 /12(1−υ 2 ) is the bending rigidity; S = κGh is the shear rigidity, and 2  is the Laplacian operator. It can be easily deduced from Eq. (8) that for thin plates where both ws and J can be taken as zero, the governing equations become: Eq. (9) is the same as the governing equation of CPT. Thus the present theory is consistent with CPT, that is, it would reduce to CPT without considering the shearing deflection and rotary inertia.
It needs to be pointed out that in this modified theory, two rotations of the transverse normal are related to the bending deflection, and the shear deformation energy will vanish when a moderately thick plate degenerates to a thin one, thus avoiding the shear locking effect.
The following conclusions can be drawn: 1) This modified two-variable shear deformation plate theory is a two-dimensional equivalent of the Timoshenko beam theory. The extension from the one-dimensional beam theory to its two-dimensional counterpart conforms to unified regularity, i.e., both of them have the same number of independent variables and governing equations.
2) In the present modified theory, the total deflection is composed of the bending and shearing defections, with rotations only related to the bending deflection, and the effect of rotary inertia is involved. In comparison, CPT only considers the bending deflection, and ignores the shearing deflection as well as rotary inertia.
3) In MPT, the rotations ψx and ψy cannot be determined solely by the deflection w, thus three independent variables exist at the median surface. The present model defines the rotations of the transverse normal by bending slopes ∂wb/∂x and ∂wb/∂y, which conforms to that of the CPT but violates the hypothesis of MPT. This is a trade-off for a simpler two-variable shear deformation approach versus the more comprehensive hypotheses of MPT involving three variables.

Finite element formulation
From Eq. (7), it can be seen that the derivatives with respect to wb and ws are of the second and the first order, respectively. Therefore, wb must be twice differentiable and C 1 -continuous, whereas ws once differentiable and C 0continuous. And the degrees of freedom required at each node of the element are wb, ∂wb/∂x, ∂wb/∂y and ws.
A new four-noded quadrilateral element with 16 degrees of freedom is thus developed to discretize the plate geometry. Dimensionless coordinates ξ and η are introduced as shown in Fig. 2. 00 , , where: x0 and y0 are the coordinates of the center of the plate in the x-y plane; ex and ey are half of the dimensions of the element in the x-and y-axes, respectively.

Fig. 2 Dimensionless coordinates of elements
The generalized displacements wb and ws are then approximated by using the following interpolation: where: Ni, Nxi, Nyi are the shape functions of a non-conforming rectangular thin plate element (ACM element) [32], and Mi the Lagrange interpolation functions, which are also the shape functions of a bilinear quadrilateral four-node (Q4) element [12]. The shape functions are: Eq. (11) can be converted into a matrix form: , , (13) in which:    (14) And the strains can be expressed as: , , (15) in which: ; ; 2 , Using the chain rule of differentiation, the derivatives with respect to x and y coordinates can be expressed as: xy xy x y x y x y x y (18) Note that the following relation is satisfied: where: |J1| is the determinant of the Jacobian matrix J1.

Free vibration analysis
For free vibration problem, the potential energy due to applied loads equals zero, thus the kinetic energy and the strain energy can be transformed into a simpler form as: 11 , , 22 eT e eT e UT == d Kd d Md (20) in which the stiffness matrix K and the mass matrix M are defined as: Applying Eq. (20) into the Hamilton's principle as Eq. (6) leads to the governing equations of free vibrations:

Eigen-buckling analysis
Consider a plate subjected to constant initial membrane stress per unit length 0  x in the x-direction, 0  y in the y-direction and a shear stress component 0  xy in x-y plane, as shown in Fig. 3.
where: 1  and 2  are the scaling parameters and σcr is the dimensionless critical buckling load.
Neglecting terms with third and higher-order derivatives, the strain field of the membrane stresses for the assumed displacement field is obtained as: For the eigen-buckling problem, the potential energy V caused by the membrane stresses 0  x , 0  y and 0  xy , can be expressed as the following integration form: where: Define the geometric stiffness matrix KG as: in which: Utilizing the Hamilton's principle as Eq. (6), the governing equations of the eigen-buckling problem can be expressed as: where: λ is the buckling factors defined as: 22 .

Boundary conditions
For the present displacement-based finite element formulation, only the displacement boundary conditions need to be strictly satisfied, which are given as: where: subscript s and n denote the tangent and the normal of the plate edge, respectively. It is noteworthy that this new quadrilateral element is established by simply adding one degree of freedom (i.e. the shear deformation ws) for each node on the basis of the ACM element. This newly developed element takes the effect of shear deformation into account, and does not exhibit shear locking at all, thus is appropriate for both thin and thick plates. Moreover, programming implementation of the new element is simple and efficient, with only minor modification to existing thin plate elements.

Numerical results and discussions
The effectiveness of this newly developed element for vibration and eigen-buckling analysis are demonstrated by examples of square plates with edge length a subjected to various boundary constraints. It is noteworthy that though the results presented in this section are for plates of square geometry only, the present element formulation can be readily applied to calculating rectangular plates. Plates are described by a symbolism denoting the corresponding boundary conditions of the four sides. For instance, SCSF implies the following boundary combinations of the four edges: S (x = 0) C (y = 0) S (x = a) F (y = a).

Free vibration analysis
In the following numerical examples, three values of shear correction factor κ (5/6, π 2 /12 and 0.8601) are used to facilitate comparison with corresponding results from other analysis published in open literature. For SSSS plates, κ is taken as 5/6. For CCCC and CCCF plates, κ = 0.8601. And for SCSC plates, κ = π 2 /12. It is found that a limited variation in the value of κ has little effect on the calculated natural frequencies.
A non-dimensional natural frequency is defined as: where: indices m and n denote the number of half-waves in the x-and y-axes, respectively. Tables 1 to 4 tabulate the first six natural frequencies of square plates obtained using the new element developed in this study, covering four types of boundary conditions, i.e., SSSS, SCSC, CCCC and CCCF. Analytical or published numerical results are also listed for comparison. All finite element solutions in these tables are generated using mesh of 20 × 20 elements.
In Tables 1 and 2, comparisons are made between the present solutions and those using Q4 elements [34], as well as the analytical solutions based on MPT [35] for SSSS and SCSC plates. An overall good agreement can be observed, with the present solutions closer to analytical solutions than those by Q4 elements. Good results of the thin plate verify that the proposed new element does not exhibit the locking effect.
Natural frequencies of CCCC plates are tabulated in Table 3, together with solutions by Q4 elements, the Rayleigh-Ritz method [36] and a meshfree method by Liew [37]. Table 4 compares present results for CCCF plates with the solutions using Q4 elements and a nine-node, quadrilateral Mindlin plate element [38], as well as results by the finite strip method (FSM) [36]. Good agreements are also observed, particularly for thin plates. For the thick geometry, the present method generates solutions slightly higher than Rayleigh-Ritz or FSM solutions for CCCC and CCCF plates, of which analytical solutions do not exist. It can be seen that the newly developed element in the present study imposes more geometrical constraints on the clamped boundary of thick plates, yielding higher natural frequencies for thick plates with clamped boundaries.
Figs. 4 to 6 provide the first four modes of vibration for square plates of different boundary conditions, with thickness ratio h/a = 0.1, illustrating the physical patterns of the four modes.
Overall, the results of the present analysis exhibit good accuracy for square plates with different thickness ratios and boundary conditions, confirming the effectiveness of present formulation. Additionally, the present element is very stable and yields no spurious zero-energy modes.  Figs. 7 to 9 show the convergence of the fundamental frequencies for SSSS, CCCC and CCCF plates of both thin (h/a = 0.01) and thick (h/a = 0.1) geometries in terms of element numbers, and good agreements with analytical solutions [35] and Rayleigh-Ritz solutions [39] are obtained. To measure the error in the finite element solution, we introduce the percentage error with respect to reference values, in which ωref is the analytical solution for SSSS and SCSC plates, and the Rayleigh-Ritz solution for CCCC and CCCF plates, ω denotes the finite element solution.
The percentage errors are also depicted in Figs. 7 to 9. It is shown that the maximum percentage error of present results is less than 1%. Well-behaved convergence characteristics of the present element are demonstrated, and the rate of convergence to the true value is fast and monotonic.

Eigen-buckling analysis
In this section, a number of numerical examples are carried out to access the new element in buckling analysis of isotropic plates. Four values of shear correction factor κ (0.822, 5/6, π 2 /12 and 0.823045) are used for a proper comparison. Both thin and moderately thick plates with various boundary combinations such as SSSS, SCSC, SCSS, SCSF, SSSF and SFSF are investigated. Different selection of χ1 and χ2 are considered, including monoaxial compressive loads in the x (χ1=-1, χ2=0) directions, and biaxial compressive loads (χ1=-1, χ2=-1).  Table 5 shows the convergence of the buckling factors for uniaxially loaded (χ1=-1, χ2=0) thin square plates (h/a=0.01) with various boundary conditions, and good agreements with analytical solutions [40] are obtained. Well-behaved convergence characteristics of present finite element formulation is clearly demonstrated, and a very good convergence rate for all the considered boundary conditions is shown. It is also shown that for the thin square plate, the ultimate convergence will be to levels slightly lower than the analytical values.
Tables 6 and 7 present the comparisons of the buckling factors with analytical solutions [41], those by Ferreira [34] using Q4 elements, and closed-form solutions by Xing and Xiang [42], for uniaxially and biaxially loaded square plates with thickness ratios h/b =0.001, 0.05, 0.1, and 0.2.
It is observed that a satisfactory agreement between them is achieved, which validates the accuracy of the present element. As expected, the decrease of critical buckling factor is observed with the increase of thickness ratio. This indicates that the shear deformation plays more and more significant roles as a plate gets thicker. Besides, the present results are generally lower than analytical solutions when only simply supported and clamped edges are involved. It is also found that the present results are generally higher than analytical solutions when free edges are involved, and decrease faster than the latter, which indicates that the present formulation weakens the stability of the plate with free edges more intensely as it gets thicker.
Buckling modes for different boundary conditions under uniaxial compression are plotted, see Fig. 10.

Conclusions
The present newly developed finite element formulation, which is based on a modified first-order shear deformation plate theory, avoids the shear-locking problem seen in other plate elements commonly used. The new formulation reduces to the thin plate theory when the rotary inertia and shearing deflection ws can be treated as negligible.
Moreover, the present formulation can be easily coded, by including the effect of shear deformations in existing thin plate elements. The model is effective and efficient in programing implementation, requiring only slight modification to existing programing of thin plate elements.
This simple four-noded element for free vibration and eigen-buckling analysis of plates manifests the same or even higher level of accuracy compared with available plate bending elements in literature or some other numerical methods widely acknowledged.
The present research may also be extended to develop more complex elements with more sophisticated geometric shapes and increasing number of nodes, so as to achieve higher level of precision for broader applications.

A NEW SHEAR LOCKING-FREE QUADRILATERAL PLATE ELEMENT FOR FREE VIBRATION AND EIGEN-BUCKLING ANALYSIS
S u m m a r y A new plate element has been formulated by simply adding one degree of freedom for each node on the basis of the classical ACM element, based on a modified first-order shear deformation plate theory, which assumes that the plate deflection has bending and shearing components, with two rotations defined by differentiations of bending deflections. This new element, with 16 degrees of freedom, exhibits no locking effect, and is appropriate for both thin and thick plates.
To assess the performance of the newly developed plate element, it has been applied to the prediction of natural frequencies and critical bucking loads of square plates with various thickness ratios and boundary conditions. Comparisons with available analytical and numerical solutions in open publication show that the new element leads to highly accurate solutions in free vibration and eigen-buckling analysis with good convergence, and can be competitive with existing plate elements.
Keywords: finite element analysis, plates and shells, shear locking, natural frequencies, buckling.