A Refined Higher-Order Hybrid Stress Quadrilateral Element for Free Vibration and Buckling Analyses of Reissner-Mindlin Plates

This paper is devoted to develop a new 8-node higher-order hybrid stress element (QH8) for free vibration and buckling analysis based on theMindlin/Reissner plate theory. In particular, a simple explicit expression of a refinemethodwith an adjustable constant is introduced to improve the accuracy of the analysis. A combined mass matrix for natural frequency analysis and a combined geometric stiffness matrix for buckling analysis are obtained using the refined method. It is noted that numerical examples are presented to show the validity and efficiency of the present element for free vibration and buckling analysis of plates. Furthermore, satisfactory accuracy for thin and moderately thick plates is obtained and it is free from shear locking for thin plate analysis and can pass the nonzero shear stress patch test.


Introduction
It is well-known that the finite element method for free vibration and buckling analysis of plates is highly significant in civil, mechanical, and aerospace engineering applications.The patch test [1], which can be used to examine the convergence of the element and construct a convergence element, has been seen as a criterion for assessing the convergence of the finite element for a long time.
The potential energy function of Mindlin plate element contains the displacements and the first-order derivatives of the displacements.According to the  0 continuity condition, it is quite easy to establish interpolation functions of deflection and rotation.Historically, the displacement-based approach was the first attempt in the formulation of effective Mindlin plate bending elements [2].However, as we all know, the original displacement element tends to cause the shearlocking phenomenon which induces overstiffness as the plate becomes progressively thinner for low-order interpolation polynomials in the Mindlin elements.During this time, the convergence characteristics of the Mindlin plate elements were performed by means of numerical computation of pure bending and pure torsion [3,4].In order to avoid shear locking, various numerical techniques and effective modifications have been proposed and tested, such as the reduced integration and selective reduced integration schemes proposed by Zienkiewicz et al. [5], Pugh et al. [6], and Hughes et al. [7]; assumed natural strain method introduced by Hughes and Tezduyar [8].It is acknowledged that the methods of reduced integration and selective integration are efficient approaches to prevent the appearance of the shear-locking phenomenon.However, it is found that such elements often exhibit extra zero-energy modes and also produce oscillatory results for some problems.Moreover, these solutions are not applicable to very thin plate; the thickness/span ratio of the plate is about restricted to 10 −6 .
Later, Belytschko et al. [9] proposed the stabilization procedure to remove the zero-energy modes by perturbing the stiffness.In addition, several efficient 9-DOF triangular elements based on the discrete Kirchhoff constraint and the equilibrium conditions were developed by Batoz et al. [10,11].These elements can eliminate locking phenomenon and converge towards the discrete Kirchhoff plate bending elements when the thickness of the plate is very thin.On the other 2 Mathematical Problems in Engineering hand, no element is free of shear locking in theory.Bathe et al. [12] introduced the MITC element and proposed strain energy patch test function to evaluate the convergence.Based on Timoshenko beam function, Soh et al. [13] proposed a triangular 9-DOF plate bending element which can be employed to analyze very thin plate (the thickness/span ratio of the plate is about 10 −11 ).Soon, Soh et al. [14] introduced a quadrilateral 12-DOF plate bending element.At this time, the progressively thinner plate which has the thickness/span ratio 10 −20 can be calculated.Wanji and Cheung [15] proposed the zero shear stress patch test functions.It is apparent that this patch test is more rigorous than the patch test using numerical computation of pure bending and pure torsion of a smallscale plate.Then, elements such as RDKQM [15], RDKTM [16], AC-MQ4 [17], and QC-P4 [18] that can pass the above patch test functions were proposed, indicating that the shearlocking problem is solved.All these elements can be used to solve the extremely thin plate problem (the thickness/span ratio of the plate can reach to 10 −30 ).In other words, these elements can accurately converge to thin plate finite element solution.
Chen proposed the enhanced patch test [19] and presented the zero shear deformation patch test and nonzero constant shear deformation test functions of Mindlin plate [20].Current patch test for Mindlin plate only satisfies the zero shear deformation condition.The patch test of nonzero constant shear for Mindlin plate problem cannot be performed.The convergence test should be performed during the process of developing finite element method.Only passing the rigorous nonzero constant shear stress patch test, the convergence can be completely guaranteed.The programs of this commercial software have no proof of convergence.The enhanced patch test is stronger than the original test; the original constant stress patch test is just a special case of it.This paper is devoted to establish Mindlin plate element which can pass the strict constant shear patch test.
Different from the classical Timoshenko beam function, Jelenić and Papa [21] proposed a new arbitrary-order Timoshenko beam function in 2011.So far, it is the only function which can be used to construct the functions of nonzero constant shear patch test for thick beam element.Since beam function can be regarded as a function on the boundary, the adopted hybrid stress method just requires the boundary function rather compared to the domain function.Because this beam function is arbitrary order, thus it has high enough order to perform the nonzero constant shear stress patch test.Since a complete cubic polynomial for the element function to pass the constant shear stress patch test is required, it was used to develop the higher-order hybrid stress triangular Mindlin plate bending element named TH6 [22] and quadrilateral Mindlin plate bending element named QH8 [23].The results of static analysis have proved that the TH6 element and QH8 element can pass the rigorous nonzero constant shear stress patch test and its accuracy is quite high.Only passing the rigorous nonzero constant shear stress patch test, the convergence can be completely guaranteed.
The purpose of this paper is to develop an 8-node Mindlin plate bending finite element for free vibration analysis and buckling analysis within an assumed stress formulation, whose main feature is that passing the rigorous nonzero constant shear stress enhanced patch test.To achieve this objective, the following steps have been taken.The first step concerns the choice of the variational framework with the adoption of complementary energy principle.Then boundary displacement interpolation function is established based on the new arbitrary-order Timoshenko beam function.Since the choice of the stress approximation is a crucial issue in developing reliable hybrid finite element, selecting a suitable stress approximation which satisfies the plate equilibrium equations is not trivial.In order to improve the performance of the constructed element, a refined mass matrix for calculation of the natural frequency and a refined geometric stiffness matrix for buckling analysis are developed by using refined element method [24][25][26].

A Brief Introduction of the Higher-Order
Hybrid Stress Quadrilateral Mindlin Plate Bending Element where , V,  are displacements along the -, -, and -axes, respectively,   ,   are the rotations of the transverse normal about the and -axes, and  is the transverse displacement field.
The geometric equations can be written as follows: where , , and  are, respectively, the rotations, the curvatures, and the shear strains: Mathematical Problems in Engineering 3 and operators B  , B  and Î are given by the following: The equilibrium equations can be obtained from the strain energy in the following form: where vectors M and S are, respectively, the moment and shear resultants: The boundary forces can be written as follows: where  is the angle between the normal of edge and the local -axis of element.For a linearly elastic material, the constitutive equations can be written as follows: where D  and D  are the elasticity matrices of bending and transverse shear moduli.In the isotropic case, the elasticity matrices specialize as where  is Young's modulus,  the shear modulus,  Poisson's ration, and  = 5/6 a correction factor to account for nonuniform distribution of shear stresses through the thickness.

Hybrid Stress Formulation.
Based on a modified complementary energy principle, the assumed stress hybrid formulation pioneered by Pian [28,29] can be used to avoid the difficulty of forming the displacement field interpolation functions, in particular, after the work of Malkus and Hughes [30] on the equivalence between reduced integration displacements and mixed/hybrid stress models.This kind of approach based on hybrid stress element method became very useful in recent years [31][32][33][34][35][36].A good number of effective elements which are free from shear locking have been developed by authors such as Tong [37], Bathe and Dvorkin [38], Ayad et al. [39], Brasile [40], and Li et al. [22,23].The higher-order hybrid stress quadrilateral Mindlin plate bending element QH8 is based on complementary energy principle.The complementary energy principle can be written as where  is the stress vector, D is the elasticity matrices, Τ is the vector of boundary force, u = [w θ θ ]  is the boundary displacement vector, w is the transverse displacement, and θ , θ are the rotations of the transverse normal about the and -axes.
The approximation for stress and boundary displacements can now be incorporated in the functional.The stress field is described in the interior of the element as follows: where P is matrix of stress interpolation functions and  is the unknown stress parameters.
The boundary force T can be represented as follows: where R is the combination of direction cosine for the boundary normal.
The boundary displacement field is described by where L are interpolation functions and q is nodal displacement parameters.Substituting the stress equation (11), boundary force equation (12), and displacement approximations equation (13) into the functional (10), where The form of ( 15) is directly amenable to numerical integration (i.e., Gauss quadrature).
Then the internal strain energy can be expressed as follows: By means of Π  /, we obtained Consequently, Substitution of  in ( 16), the internal strain energy reduces to Compared with  = (1/2)q  Kq, the element stiffness matrix can be taken as The solution of the system yields the unknown nodal displacement q.After q is determined, element stress or internal forces can be recovered by use of ( 18) and (11).Thus  locking; it cannot solve the problem of passing the nonzero constant shear patch test.This problem has not resolved for many years.In 2011, Jelenić and Papa [21] presented a new arbitrary-order Timoshenko beam function as follows: where  is the beam length,   and   are the values of the displacements and the rotations at the  nodes equidistantly spaced between the beam ends,   are the standard Lagrange polynomials of order  − 1, and   = / for  = 1 and   = 1 − ( − 1)/( − 1)(/) otherwise, in which  is the coordinate along the beam.An 8-node quadrilateral element was designed as given in Figure 1.If any quadrilateral side is taken as a beam element, take the 1-2 boundary as example and the deflection w and rotations θ , θ can be derived as follows: where is the length of 1-2 boundary and  is the coordinate along the 1-2 edge.

Assumed Stresses.
In this section, a strategy to select the stress approximation in a rational way is presented.In practice, initial polynomials are usually assumed for the stress after which the equilibrium equations are applied to these polynomials yielding relations between the 's and ultimately the final form of P. The number of stress parameters, which is the number of columns in P, must be at least equal to the number of degrees of freedom of the element less the number of degrees of freedom necessary to prevent rigid body motion.
In a series of studies of hybrid element, there is a general consensus on the selection of  parameter that is () ≥   −   , while the optimal selection is () =   −   ,   ,   being the degrees of freedom of the element and the number of allowed rigid body motions, respectively.The element has 24 DOF; therefore, a stress field with at least 21 parameters is needed to describe the stress field and without spurious zero-energy modes.Numerical results show that the stiffness matrix has spurious zero-energy modes and converges slowly when () = 21.The finds that gradually increasing the number of {  } until () = 28, there is a proper rank for the stiffness matrix and the absence of spurious zero-energy modes.Moreover, the calculation results are accurate and converged faster.And they are better if the number of {  } reaches   = 39 (complete quartic polynomial).Numerical experimentations [23] indicate that this 39-parameter selection of stress field is somewhat more accurate and has no spurious zero-energy modes.For this reason, 39 are chosen as the assumed stress field.
The stress field is described in the interior of the element as follows: where P is matrix of stress interpolation functions as shown in Table 1, and  is the stress parameters.From (12), the boundary force T can be expressed as Then, along the boundary 1-2, G can be obtained: Similarly, along the other boundaries, G 23 , G 34 , G 41 can also be obtained by cyclic permutation.

The Refined Mass Matrix and Refined Geometric Stiffness Matrix for Vibration and Buckling Analysis of Plates
A refined method has been developed by the authors [24,25], and the method has been applied to develop 12-DOF [15], 9-DOF [16], and 5-DOF [49] Mindlin elements and so forth.This method is applied here to formulate QH8-RCV with a mass matrix called refined mass matrix for calculation of the natural frequency and in the same way QH8-R with a refined geometric stiffness matrix for buckling analysis.

Free Vibration Analysis.
It is well-known that the equation of free vibration can be expressed as follows: where K is the global stiffness matrix, M is the global mass matrix, Φ is the mode shapes, and  is the natural frequencies.The major complication of the standard eigenvalue problem introduced by the finite element method is the mass matrix M. Two types of the mass matrix are of interest: lumped mass matrix and consistent mass matrix.

The Lumped Mass Matrix.
To use the diagonal elements of the consistent mass matrix and form a diagonal matrix by scaling these entries so that the total mass of the element is conserved, the method listed by Hinton et al. [50] produces excellent results.For the cases of eight-node elements, the method leads to the lumping schemes shown in Figure 2. The figures at each node show the proportion of the total mass at that node.For convenience, the total mass has been chosen to be 36.Therefore the following lumped mass matrix will be used to calculate the natural frequencies of plates: where   is the element mass.

The Refined Consistent Mass Matrix.
In order to improve the accuracy of the natural frequencies, the refined mass matrix has been used in the eigenproblem of the finite element method.We can give herein a new approach for improving the consistent mass matrix which can be obtained from a combination of the element displacement interpolations.
As for free vibration analysis, 8-node shape function is used to interpolate the displacement function components ,   ,   as where Another 4-node shape functions are used to interpolate the deflection part of displacement function as where So the new deflection function part  * can be described as where  is an adjustable constant which is used to improve the accuracy of the vibration analysis.
Based on ( 31), (33), and (35), the shape function of refined consistent mass matrix is of the following form: and the element consistent mass matrix M   can be expressed as where m is the matrix containing the mass density of the material  and thickness  as follows: ) . (38)

Buckling Analysis Problem and Refined Geometric Stiffness
Matrix.Consider a plate under an initial stress which may induce instability.Similar to the free vibration problem of plate, the equation of stability analysis can be expressed as where K is the global stiffness matrix, K  is the geometric stiffness matrix of the whole structure,  is the critical load, and Φ is the mode shape of buckling.The element geometric stiffness matrix K   can be written as where the terms of ( 40) are and vector ] stands for the initial stress of plate,  0  ,  0  are unit tractions in the direction of  and , and  0  is the unit shear stress.Similarly, we can obtain the refined geometric stiffness matrix from the following relation: where  0 is an adjustable constant used to improve the accuracy of the stability analysis, the displacement interpolation functions ,   ,   are the same as 8-node shape function (31), and   is in terms of (33).
Based on (42), we only modify the interpolation function  and geometric stiffness matrix part K  1 as The refined geometric stiffness matrix can be written as

Numerical Examples
To investigate the accuracy and reliability of the proposed QH8 element for free vibration analysis and buckling analysis, several numerical examples with different geometry and boundary condition are calculated in this section.For convenience, the boundaries of plates are denoted as follows: clamped supported (C), hard-type simply supported (S), softtype simply supported (S * ), and free (F), so a boundary condition can be written as CCCC, SSSS, CFFF, SCSC, SFSF, and so on for a rectangular plate.For example, the symbol, SFCF, represents simply supported edge along  = 0, a free edge along  = 0, a clamped edge along  = , and a free edge along  = .
4.1.Patch Test: Consistency Assessment.Consistency of the developed elements is tested for the constant strain and stress states on the patch example with five elements, covering a rectangular domain of a plate as shown in Figure 3.The size of the domain is 0.24 × 0.12.The mechanical properties of the plate are chosen as  = 10 3 ,  = 0.25,  = 5/6.For zero shear deformation patch test, the test function [20,23] can be written as follows: For nonzero shear deformation patch test, the test function can be expressed as follows: where   , ( = 0, . . ., 9) are arbitrary constants.In the present paper, we assume that   =  + 1.
The mesh of patch test is presented in Figure 3.The displacements at the nodes on the boundaries are imposed according to (45) and (46), given the displacements and rotations at the boundary nodes (8 displacements and 16 rotations), while all the internal nodal displacements and rotations are to be calculated by the finite element solution procedure.
Test results show that the element QH8-39 passes both the constant bending with zero shear stresses   ,   (Table 2) and the strict patch test with nonzero constant shear stresses   ,   (Table 3).No matter the inner element edges are straight or curved, and no matter the shapes of the elements are convex or concave, the exact results of the displacements and stresses at inner node obtained using the present QH8-39 element agree well with the test functions.This demonstrates that the new element passes the strict constant stress patch test, thus ensuring solution convergence.

Free Vibration Analysis.
In this section, we examine the efficiency and the applicability of the 8-node quadrilateral hybrid stress element in analyzing natural frequencies of plates.The proposed 8-node quadrilateral hybrid stress elements with lumped mass matrix (30) and refined consistent mass matrix (37) in vibration analysis are denoted by QH8-LV and QH8-RCV, respectively.
The effects of various slenderness ratios, skew angles, and boundary conditions on the frequencies are also discussed.Unless other stated, the following material parameters are used throughout the paper: Young's modulus  = 2.0 × 10 11 , Poisson's ratio  = 0.3, the mass density of the material  = 8000, and the shear correction factor  = 5/6 is taken if not specified otherwise.In the numerical calculation, the dimensionless solutions of natural frequencies with a multiplier  2 / 2 √/ are obtained. is the width of the edge,  =  3 /12(1 −  2 ).Unless otherwise specified, the actual computations are done using dimensionless units.

Square Plate.
For the hybrid stress element with the lumped mass matrix (QH8-LV) and refined mass matrix (QH8-RCV), the fundamental natural frequencies or free vibration problems of a square plate with different meshes and different boundary conditions are calculated.
The parameter  in (36) for the present element QH8-RCV (8-node hybrid element with refined mass matrix) is determined by comparison using numerical examples.A square plate (the length of the plate  and the width of the plate : / = 1) is divided into mesh of 16 × 16 elements.The thickness-to-width aspect ratios / = 0.001, 0.1, boundary conditions SSSS, CCCC, and  = −1, −0.5, −0.25, −0.08, 0 are chosen.The results of six lowest modes are given in Tables 4-7 and Figure 4, in which the exact solutions can be found in [27].We advise  = −0.08 of the present QH8-RCV() for following numerical examples.
The convergence of frequencies results of thin and thick square plates with SSSS and CCCC boundary conditions for QH8-RCV(−0.08) and QH8-LV are shown in Table 8 and Figure 5.The square plate is discretized by  ×  uniform meshes with  = 4, 8, 16.The results agree well with [27] when  = 8, 16.
Then, natural frequencies analyses are further carried on the thin (/ = 0.001) and thick (/ = 0.1) square plate with the meshes of 16 × 16 elements and various boundary conditions such as SSFF, SFSF, SCSC, CFFF, CCFF, and CFCF.The results are summarized in Table 9 and agree well with [27].The six lowest modes of clamped and simply supported (CCCC, SSSS) thin square plate obtained by the QH8-RCV are also investigated and plotted in Figure 6.

Free Vibration of Circular Plate.
Free vibration of a circular plate (Figure 7) of  = 5 with the clamped boundary is then studied in this subsection.The material parameters of the circular plate are given as follows: the density mass  = 8000, Young's modulus  = 2.0 × 10 11 , and Poisson's ratio  = 0.3.60 proposed elements and 205 nodes as shown in Figure 7, Mesh A, are used to discretize the whole plate for S8R of Abaqus, QH8-RCV(−0.08) and QH8-LV, while the compared solutions given by other elements employ even finer meshes.The frequencies given by MIN3 and NS+ES-FEM [51] employ a mesh of 394 elements and 222 nodes; the solutions of NS-DSG and ES-FEM+DSG3 employ a much finer mesh of 848 elements and 460 nodes and ANS4 [52] using 432 quadrilateral plate elements (or 864 triangular elements).
The first ten modes of thin (/2 = 0.01) and thick (/2 = 0.1) circular plates are shown in Tables 10 and 11, respectively.Figure 8 plots the results obtained from different methods.It is observed that the QH8-RCV and QH8-LV are good competitors to ANS4 and exact solutions [42,43] and is more accurate than MIN3, ES-FEM+DSG3, and NS+ES-FEM.
First six lowest modes of this clamped circular obtained from QH8-RCV using the mesh of Figure 7 Mesh B are shown in Figure 9, which exactly describes the real physical modes of the circular plate.

Free Vibration of Triangular
Plates.QH8-RCV ( = −0.08)and QH8-LV are then applied in the free vibration analyses of a cantilever triangular plate (CFF) with various shape geometries just as shown in Figure 10 nondimensional frequency parameter  =  2 √// 2 for thin (/ = 0.001) and thick (/ = 0.2) triangular plates against different skew angles  = 0 ∘ , 15 ∘ , 30 ∘ , 45 ∘ , and 60 ∘ is undertaken and shown in Tables 12 and 13, respectively.The first four nondimensional frequency parameters for the thin and thick triangular plates are plotted in Figures 11 and 12.
The NS+ES-FEM [51] with the mesh of 170 triangular elements and 108 nodes and DSG3, MIN3, and ES-FEM+DSG3 are also used for comparison.The reference solution was obtained by the ANS method [52] using a mesh of 796 triangular elements.For thick triangular plate, QH8-RCV's performance is the best, but QH8-LV provides worse result when the angle becomes larger for the third and fourth modes in Figure 12, because of the mesh with small number of elements (45 ∘ , 26 elements and 60 ∘ , 28 elements) in some degree.For the thin triangular plate, QH8-LV is better than DSG3 and MIN3 and competes well with NS+ES-FEM and Rayleigh-Ritz [53].QH8-RCV's performance is still the best and a strong competitor of Rayleigh-Ritz.Six lowest modes of cantilever triangular plates (CFF) obtained from QH8-RCV are illustrated in Figure 13.

Critical Load of Buckling.
In the following tests, we perform a series of plate buckling analyses to assess the convergence characteristics and accuracy of the developed QH8 element in predicting a critical buckling load.The results are presented in terms of a nondimensional buckling load intensity factor.
with  being the edge width of the plate that is used throughout the paper and  =  3 /12(1 −  2 ).

Rectangular Plates Subjected to an In-Plane Compressive
Load.The parameter  0 in (42) for the present element      A shear-locking phenomenon is firmly involved when using thick plate theories to analyze thin plates.The present QH8 element is free of shear-locking in static calculation.A square plate is considered for this purpose and all other related parameters are taken exactly the same as above.Various boundaries such as SSSS, SFSC, SFSS, and SFSF  Table 14: The factor of uniaxial buckling loads along the -axis of rectangular plates with length-to-width aspect ratios / = 1 and thicknessto-width aspect ratios / = 0.01.are investigated.Several thickness-span aspect ratios / = 0.1, 0.2 (moderately thick plate) and / = 0.001, 0.05 (thin plate) are considered, respectively, and the results are shown in Table 16 and Figure 15, which are compared with the analytical solutions [45] and other existing numerical approaches such as the SSM [54] and the shear-locking-free and mesh-free method (Shear-MK, MK) [55].The element S8R and QH8-R use a regular pattern of 11 × 11 element (408 Table 15: The factor of uniaxial buckling loads along the -axis of rectangular plates with length-to-width aspect ratios / = 1 and thicknessto-width aspect ratios / = 0.1.Through the achieved results shown in Table 16 and Figure 15, the shear-locking phenomenon, on the one side, is found very clearly because of losing the accuracy significantly when the standard MK is applied to deal with thin plates, that is, / = 0.001, 0.05, but conversely matches well with thick plates, that is, / = 0.1, 0.2.Higher buckling loads by the standard MK are obtained because the plates may exhibit highly stiffer behavior.On the other hand, a good agreement and free of shear-locking are found when the Shear-MK, S8R, and the present QH8-R are used.

Compression Buckling Behavior of Skew Plates.
The buckling problem of skew plate has also been considered by Kitipornchai et al. [46] and Wang et al. [56] using the pb-2 Ritz method.Skew plates, as shown in Figure 20, with skew angle, , thickness-to-width ratio, /, and different combinations of edge support conditions, are examined.The plate is modelled with 6 × 6, 9 × 9 and 16 × 16 distributed particles for QH8-R ( 0 = −0.17)element.
Numerical results for thick plates with thickness-to-width ratios, / = 0.05, 0.1, skew angles,  = 0 ∘ , 15 ∘ , 30 ∘ , 45 ∘ , and various boundary conditions of S * S * S * S * , CCCC, CFCF, and S * FS * F, are compared with those of mesh-free method [48] (with 17×17 distributed particles) and Kitipornchai et al. [46] in Table 20.S * , C, and F denote soft-type simply supported, clamped, and free boundary conditions, respectively.Based on the favorable comparisons observed in Table 8, it can be concluded that the plate buckling load intensity factors are well approximated by the QH8-R ( 0 = −0.17)element.

A Square Plate with a Hole Subjected to Different
In-Plane Loads.Lastly we consider a hard-type simply supported (SSSS) and clamped (CCCC) square plate with a hole.For the purpose of the comparison, the material parameters are taken the same as in the above example but the length  =  = 10 and the thickness  = 1 are employed instead.The irregularly scattered pattern of 84 elements and 308 nodes shown in Figure 21 is used.The plate is subjected to three different in-plane load cases, namely, an uniaxial compressive load in -direction,    21.The results agree very well with those computed by the EFG [57] involving both the thirdorder shear deformation plate theory (TSDT) and the FSDT.And the Shear-MK [55] with 400 nodes is also used for comparison.Additionally, it is also interesting to view the buckling modes of the problem and, thus, Figure 22 presents the first buckling modes under the three types of buckling loads for four different boundary conditions.

Conclusions
The higher-order hybrid stress quadrilateral Mindlin plate element QH8 for solving free vibration and buckling analysis problems is formulated in this paper.Due to numerical results given above, some concluding remarks can be drawn as follows.
The QH8 element is an efficiency method for solving the free vibration problems and buckling analysis of plate due to its high accuracy.It is free of shear locking and can pass both the constant bending with zero shear stresses and the strict patch test with nonzero constant shear stresses, no matter the element edges are straight or curved and no matter the shapes of the elements are convex or concave.
The refined nonconforming element method is a very efficient method, which can be used to improve the accuracy of nonconforming elements.It has been used to construct the QH8 element successfully.With the refined mass matrix and refined geometric stiffness matrix, the accuracy of vibration and stability analysis has been improved just by changing the value of  and  0 , and it really works.
From the numerical results, the present plate element has performed well in most situations.It is shown that the present plate element is applicable to either thin or thick situations with enough accuracy.Note.QH8-R ( 0 = −0.17Note.QH8-R ( 0 = −0.17).a Shear correction factor  = 0.823045.b Shear correction factor  = 5/6.
Table 19: The factor of uniaxial buckling loads along the -axis of rectangular plates (SSSS) with various length-to-width ratios and various thickness-to-width ratios.

Figure 2 :
Figure 2: Mass lumping at the nodes.

Figure 4 :
Figure 4: Error of six lowest nondimensional frequency parameters.

Figure 8 :
Figure 8: Ten lowest parameterized natural frequencies of a clamped circular plate.

Figure 11 :
Figure 11: Convergence study of the lowest four parameterized natural frequencies of triangular plates with / = 0.001.

Figure 12 :Figure 13 :Figure 14 :
Figure 12: Convergence study of the lowest four parameterized natural frequencies of triangular plates with / = 0.2.

Figure 19 :
Figure 19: Axial buckling modes of simply supported rectangular plates with thickness-to-width ratios / = 0.01 and various length-to-width ratios /.

Figure 21 :
Figure 21: Square plate with a circular hole and the distributions (308 nodes, 84 elements).

Figure 22 :
Figure 22: The first uniaxial, biaxial, and pear shear buckling modes of the plate with a circular hole associated with different boundaries as SSSS, CCCC (thickness  = 1.0).

Table 4 :
The six lowest nondimensional frequency parameters of a CCCC thin square plate (/ = 0.001).

Table 5 :
The six lowest nondimensional frequency parameters of a CCCC thick square plate (/ = 0.1).

Table 6 :
The six lowest nondimensional frequency parameters of a SSSS thin square plate (/ = 0.001).

Table 7 :
The six lowest nondimensional frequency parameters of a SSSS thick square plate (/ = 0.1).

Table 9 :
Frequency parameters for square Mindlin plates with different boundary conditions.

Table 12 :
The lowest six parameterized natural frequencies of triangular plates with / = 0.001.

Table 13 :
The lowest six parameterized natural frequencies of triangular plates with / = 0.2.
0  , pure shear load, and biaxial compressive loads in  and  directions, where

Table 17 :
The factor of uniaxial buckling loads (  0 ) along the y-axis of rectangular plates (/ = 1) with various boundary conditions and various thickness-to-width ratios (/).

Table 21 :
The uniaxial, biaxial, and pure shear buckling load factor for a square plate with a circular hole and different boundary conditions.