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

Authors: Tan Li; Xu Ma; Wanji Chen; Xili Jing
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

This paper is devoted to develop a new 8-node higher-order hybrid stress element (QH8) for free vibration and buckling analysis based on the Mindlin/Reissner plate theory. In particular, a simple explicit expression of a refine method with 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.


---

## Body

## 1. 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 theC0 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 shear-locking 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 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 small-scale 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 shear-locking 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–26].


## 2. A Brief Introduction of the Higher-Order Hybrid Stress Quadrilateral Mindlin Plate Bending Element

### 2.1. Fundamental Equations of Mindlin Plate

Consider a plate referred to as a Cartesian coordinate frameo,x,y,z, with the origin o on the mid-surface Ω and the z-axis in the thickness direction, -h/2≤z≤h/2, where h is the plate thickness. Let ∂Ω be the boundary of Ω. The Reissner-Mindlin theory, that is, the first-order shear deformable theory, is employed. Thus it is assumed that (1)ux,y,z=zθy,vx,y,z=-zθx,wx,y,z=wx,y,where u, v, w are displacements along the x-, y-, and z-axes, respectively, θx, θy are the rotations of the transverse normal about the x- and y-axes, and w is the transverse displacement field.The geometric equations can be written as follows:(2)χ=Bbθ,γ=Bsw+I^θ,where θ, χ, and γ are, respectively, the rotations, the curvatures, and the shear strains: (3)θ=θxθy,χ=χxχyχxy,γ=γxzγyz,and operators Bb, Bs and I^ are given by the following:(4)Bb=0∂∂x-∂∂y0-∂∂x∂∂y,Bs=∂∂x∂∂y,Ι^=01-10.The equilibrium equations can be obtained from the strain energy in the following form:(5)-BbTM+I^TS=0,BsTS=0,where vectors M and S are, respectively, the moment and shear resultants:(6)M=MxMyMxy,S=SxSy.The boundary forces can be written as follows:(7)T=RMS,R=0-sin⁡α-cos⁡α00cos⁡α0sin⁡α00000cos⁡αsin⁡α,where α is the angle between the normal of edge and the local x-axis of element.For a linearly elastic material, the constitutive equations can be written as follows:(8)M=Dbχ,S=Dsγ,where Db and Ds are the elasticity matrices of bending and transverse shear moduli. In the isotropic case, the elasticity matrices specialize as(9)Db=Eh3121-μ21μ0μ10001-μ2,Ds=κGh1001,where E is Young’s modulus, G the shear modulus, μ Poisson’s ration, and κ=5/6 a correction factor to account for nonuniform distribution of shear stresses through the thickness.


### 2.2. 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–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(10)Πe=∫Ωe12σTD-1σdΩ-∫∂ΩeΤTu-dS,where σ is the stress vector, D is the elasticity matrices, Τ is the vector of boundary force, u-=w~θ~xθ~yT is the boundary displacement vector, w~ is the transverse displacement, and θ~x, θ~y are the rotations of the transverse normal about the x- and y-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:(11)MxMyMxySxSyT=Pβ,where P is matrix of stress interpolation functions and β is the unknown stress parameters.The boundary forceT can be represented as follows:(12)T=RPβ,where R is the combination of direction cosine for the boundary normal.The boundary displacement field is described by(13)u-=w~θ~xθ~yT=Lq,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),(14)Πe=12βTHβ-βTGq,where(15)H=∫ΩePTD-1PdΩG=∫∂ΩeRPTLdS.The form of (15) is directly amenable to numerical integration (i.e., Gauss quadrature).Then the internal strain energy can be expressed as follows:(16)U=∫Ωe12σTD-1σdΩ=12βTHβ.By means of∂Πe/∂β, we obtained(17)Hβ=Gq.Consequently,(18)β=H-1Gq.Substitution ofβ in (16), the internal strain energy reduces to(19)U=12qTGTH-1Gq.Compared with U=1/2qTKq, the element stiffness matrix can be taken as(20)K=GTH-1G.The solution of the system yields the unknown nodal displacementq. After q is determined, element stress or internal forces can be recovered by use of (18) and (11). Thus(21)σ=PH-1Gq.


### 2.3. The Displacement Interpolation Function of QH8 Element

Euler-Bernoulli beam function has been successfully employed in the construction of refined thin plate elements. It is well-known that when constructing Mindlin plate element, both thick and thin plates should be taken into account, and it is necessary to eliminate shear-locking phenomenon. To seek out such element displacement function is definitely very difficult. Note that a closed form solution for both thick and thin beams exists in the form of the Timoshenko beam function, and it is possible to use it to derive more efficient Mindlin plate elements [41]. However, the use of Timoshenko beam function is capable of solving the problem of shear 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:(22)w=∑i=1nIiwi-Ln∏j=1nNj∑i=1n-1i-1n-1i-1θiθ=∑i=1nIiθi,where L is the beam length, wi and θi are the values of the displacements and the rotations at the n nodes equidistantly spaced between the beam ends, Ii are the standard Lagrange polynomials of order n-1, and Nj=x/L for j=1 and Nj=1-n-1/j-1x/L otherwise, in which x is the coordinate along the beam.An 8-node quadrilateral element was designed as given in Figure1. If any quadrilateral side is taken as a beam element, take the 1-2 boundary as example and the deflection w~ and rotations θ~x, θ~y can be derived as follows:(23)w~=I1w1+I2w2+I5w5-I0L-θn,1+2θn,5-θn,2=I1w1+I2w2+I5w5-I0-θx,1+2θx,5-θx,2b+-θy,1+2θy,5-θy,2aθ~x=I1θx,1+I2θx,2+I5θx,5θ~y=I1θy,1+I2θy,2+I5θy,5,where I1=L12L1-1, I2=L22L2-1, I5=4L1L2, I0=1/3L1L2L2-L1, L1=1-s/L, L2=s/L, L is the length of 1-2 boundary and s is the coordinate along the 1-2 edge.Figure 1
Eight-node quadrilateral plate element.The displacement componentsu- along the 1-2 boundary was given as follows:(24)u-1-2=w~θ~xθ~yT=L1-2q1-2,where(25)L1-2=I1I0bI0a0I1000I1I2I0bI0a0I2000I2I5-2I0b-2I0a0I5000I5and q1-2=q1q2q5T, qi=wiθxiθyiT(i=1,2,5).


### 2.4. 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 Nβ≥nd-nr, while the optimal selection is Nβ=nd-nr, nd, nr 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 Nβ=21. The finds that gradually increasing the number of βi until Nβ=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 βi reaches Nb=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:(26)MxMyMxyQxQyT=Pβ,where P is matrix of stress interpolation functions as shown in Table 1, and β is the stress parameters.Table 1
Assumed stress function.
i
M
x
M
y
M
x
y
Q
x
Q
y
1
1
0
0
0
0
2
0
1
0
0
0
3
0
0
1
0
0
4
x
0
0
1
0
5
y
0
0
0
0
6
0
x
0
0
0
7
0
y
0
0
1
8
0
0
x
0
1
9
0
0
y
1
0
10
x
2
0
-
x
y
x
-
y
11
x
y
0
0
y
0
12
y
2
0
0
0
0
13
0
x
2
0
0
0
14
0
x
y
0
0
x
15
0
y
2
-
x
y
-
x
y
16
0
0
x
2
0
2
x
17
0
0
y
2
2
y
0
18
x
3
0
-
1.5
x
2
y
1.5
x
2
-
3
x
y
19
x
2
y
0
-
0.5
x
y
2
x
y
-
0.5
y
2
20
x
y
2
0
0
y
2
0
21
y
3
0
0
0
0
22
0
x
3
0
0
0
23
0
x
2
y
0
0
x
2
24
0
x
y
2
-
0.5
x
2
y
-
0.5
x
2
x
y
25
0
y
3
-
1.5
x
y
2
-
3
x
y
1.5
y
2
26
0
0
x
3
0
3
x
2
27
0
0
y
3
3
y
2
0
28
x
4
-
6
x
2
y
2
0
4
x
3
-
12
x
2
y
29
x
3
y
0
-
0.75
x
2
y
2
1.5
x
2
y
-
1.5
x
y
2
30
x
y
3
0
0
y
3
0
31
y
4
0
0
0
0
32
0
x
4
0
0
0
33
0
x
3
y
0
0
x
3
34
0
x
y
3
-
0.75
x
2
y
2
-
1.5
x
2
y
1.5
x
y
2
35
-
6
x
2
y
2
y
4
0
-
12
x
y
2
4
y
3
36
0
0
x
4
0
4
x
3
37
0
-
3
x
2
y
2
x
3
y
x
3
-
3
x
2
y
38
-
3
x
2
y
2
0
x
y
3
-
3
x
y
2
y
3
39
0
0
y
4
4
y
3
0From (12), the boundary force T can be expressed as(27)T=RPβ.Then, along the boundary 1-2,G can be obtained: (28)G12=∫∂Ω1-2RPTL1-2ds.Similarly, along the other boundaries, G23, G34, G41 can also be obtained by cyclic permutation.


## 2.1. Fundamental Equations of Mindlin Plate

Consider a plate referred to as a Cartesian coordinate frameo,x,y,z, with the origin o on the mid-surface Ω and the z-axis in the thickness direction, -h/2≤z≤h/2, where h is the plate thickness. Let ∂Ω be the boundary of Ω. The Reissner-Mindlin theory, that is, the first-order shear deformable theory, is employed. Thus it is assumed that (1)ux,y,z=zθy,vx,y,z=-zθx,wx,y,z=wx,y,where u, v, w are displacements along the x-, y-, and z-axes, respectively, θx, θy are the rotations of the transverse normal about the x- and y-axes, and w is the transverse displacement field.The geometric equations can be written as follows:(2)χ=Bbθ,γ=Bsw+I^θ,where θ, χ, and γ are, respectively, the rotations, the curvatures, and the shear strains: (3)θ=θxθy,χ=χxχyχxy,γ=γxzγyz,and operators Bb, Bs and I^ are given by the following:(4)Bb=0∂∂x-∂∂y0-∂∂x∂∂y,Bs=∂∂x∂∂y,Ι^=01-10.The equilibrium equations can be obtained from the strain energy in the following form:(5)-BbTM+I^TS=0,BsTS=0,where vectors M and S are, respectively, the moment and shear resultants:(6)M=MxMyMxy,S=SxSy.The boundary forces can be written as follows:(7)T=RMS,R=0-sin⁡α-cos⁡α00cos⁡α0sin⁡α00000cos⁡αsin⁡α,where α is the angle between the normal of edge and the local x-axis of element.For a linearly elastic material, the constitutive equations can be written as follows:(8)M=Dbχ,S=Dsγ,where Db and Ds are the elasticity matrices of bending and transverse shear moduli. In the isotropic case, the elasticity matrices specialize as(9)Db=Eh3121-μ21μ0μ10001-μ2,Ds=κGh1001,where E is Young’s modulus, G the shear modulus, μ Poisson’s ration, and κ=5/6 a correction factor to account for nonuniform distribution of shear stresses through the thickness.


## 2.2. 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–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(10)Πe=∫Ωe12σTD-1σdΩ-∫∂ΩeΤTu-dS,where σ is the stress vector, D is the elasticity matrices, Τ is the vector of boundary force, u-=w~θ~xθ~yT is the boundary displacement vector, w~ is the transverse displacement, and θ~x, θ~y are the rotations of the transverse normal about the x- and y-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:(11)MxMyMxySxSyT=Pβ,where P is matrix of stress interpolation functions and β is the unknown stress parameters.The boundary forceT can be represented as follows:(12)T=RPβ,where R is the combination of direction cosine for the boundary normal.The boundary displacement field is described by(13)u-=w~θ~xθ~yT=Lq,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),(14)Πe=12βTHβ-βTGq,where(15)H=∫ΩePTD-1PdΩG=∫∂ΩeRPTLdS.The form of (15) is directly amenable to numerical integration (i.e., Gauss quadrature).Then the internal strain energy can be expressed as follows:(16)U=∫Ωe12σTD-1σdΩ=12βTHβ.By means of∂Πe/∂β, we obtained(17)Hβ=Gq.Consequently,(18)β=H-1Gq.Substitution ofβ in (16), the internal strain energy reduces to(19)U=12qTGTH-1Gq.Compared with U=1/2qTKq, the element stiffness matrix can be taken as(20)K=GTH-1G.The solution of the system yields the unknown nodal displacementq. After q is determined, element stress or internal forces can be recovered by use of (18) and (11). Thus(21)σ=PH-1Gq.


## 2.3. The Displacement Interpolation Function of QH8 Element

Euler-Bernoulli beam function has been successfully employed in the construction of refined thin plate elements. It is well-known that when constructing Mindlin plate element, both thick and thin plates should be taken into account, and it is necessary to eliminate shear-locking phenomenon. To seek out such element displacement function is definitely very difficult. Note that a closed form solution for both thick and thin beams exists in the form of the Timoshenko beam function, and it is possible to use it to derive more efficient Mindlin plate elements [41]. However, the use of Timoshenko beam function is capable of solving the problem of shear 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:(22)w=∑i=1nIiwi-Ln∏j=1nNj∑i=1n-1i-1n-1i-1θiθ=∑i=1nIiθi,where L is the beam length, wi and θi are the values of the displacements and the rotations at the n nodes equidistantly spaced between the beam ends, Ii are the standard Lagrange polynomials of order n-1, and Nj=x/L for j=1 and Nj=1-n-1/j-1x/L otherwise, in which x is the coordinate along the beam.An 8-node quadrilateral element was designed as given in Figure1. If any quadrilateral side is taken as a beam element, take the 1-2 boundary as example and the deflection w~ and rotations θ~x, θ~y can be derived as follows:(23)w~=I1w1+I2w2+I5w5-I0L-θn,1+2θn,5-θn,2=I1w1+I2w2+I5w5-I0-θx,1+2θx,5-θx,2b+-θy,1+2θy,5-θy,2aθ~x=I1θx,1+I2θx,2+I5θx,5θ~y=I1θy,1+I2θy,2+I5θy,5,where I1=L12L1-1, I2=L22L2-1, I5=4L1L2, I0=1/3L1L2L2-L1, L1=1-s/L, L2=s/L, L is the length of 1-2 boundary and s is the coordinate along the 1-2 edge.Figure 1
Eight-node quadrilateral plate element.The displacement componentsu- along the 1-2 boundary was given as follows:(24)u-1-2=w~θ~xθ~yT=L1-2q1-2,where(25)L1-2=I1I0bI0a0I1000I1I2I0bI0a0I2000I2I5-2I0b-2I0a0I5000I5and q1-2=q1q2q5T, qi=wiθxiθyiT(i=1,2,5).


## 2.4. 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 Nβ≥nd-nr, while the optimal selection is Nβ=nd-nr, nd, nr 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 Nβ=21. The finds that gradually increasing the number of βi until Nβ=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 βi reaches Nb=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:(26)MxMyMxyQxQyT=Pβ,where P is matrix of stress interpolation functions as shown in Table 1, and β is the stress parameters.Table 1
Assumed stress function.
i
M
x
M
y
M
x
y
Q
x
Q
y
1
1
0
0
0
0
2
0
1
0
0
0
3
0
0
1
0
0
4
x
0
0
1
0
5
y
0
0
0
0
6
0
x
0
0
0
7
0
y
0
0
1
8
0
0
x
0
1
9
0
0
y
1
0
10
x
2
0
-
x
y
x
-
y
11
x
y
0
0
y
0
12
y
2
0
0
0
0
13
0
x
2
0
0
0
14
0
x
y
0
0
x
15
0
y
2
-
x
y
-
x
y
16
0
0
x
2
0
2
x
17
0
0
y
2
2
y
0
18
x
3
0
-
1.5
x
2
y
1.5
x
2
-
3
x
y
19
x
2
y
0
-
0.5
x
y
2
x
y
-
0.5
y
2
20
x
y
2
0
0
y
2
0
21
y
3
0
0
0
0
22
0
x
3
0
0
0
23
0
x
2
y
0
0
x
2
24
0
x
y
2
-
0.5
x
2
y
-
0.5
x
2
x
y
25
0
y
3
-
1.5
x
y
2
-
3
x
y
1.5
y
2
26
0
0
x
3
0
3
x
2
27
0
0
y
3
3
y
2
0
28
x
4
-
6
x
2
y
2
0
4
x
3
-
12
x
2
y
29
x
3
y
0
-
0.75
x
2
y
2
1.5
x
2
y
-
1.5
x
y
2
30
x
y
3
0
0
y
3
0
31
y
4
0
0
0
0
32
0
x
4
0
0
0
33
0
x
3
y
0
0
x
3
34
0
x
y
3
-
0.75
x
2
y
2
-
1.5
x
2
y
1.5
x
y
2
35
-
6
x
2
y
2
y
4
0
-
12
x
y
2
4
y
3
36
0
0
x
4
0
4
x
3
37
0
-
3
x
2
y
2
x
3
y
x
3
-
3
x
2
y
38
-
3
x
2
y
2
0
x
y
3
-
3
x
y
2
y
3
39
0
0
y
4
4
y
3
0From (12), the boundary force T can be expressed as(27)T=RPβ.Then, along the boundary 1-2,G can be obtained: (28)G12=∫∂Ω1-2RPTL1-2ds.Similarly, along the other boundaries, G23, G34, G41 can also be obtained by cyclic permutation.


## 3. 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.


### 3.1. Free Vibration Analysis

It is well-known that the equation of free vibration can be expressed as follows:(29)K-ω2MΦ=0,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 matrixM. Two types of the mass matrix are of interest: lumped mass matrix and consistent mass matrix.


#### 3.1.1. 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: (30)Mlumpede=me36diag⁡100100100100800800800800,where me is the element mass.Figure 2
Mass lumping at the nodes.


#### 3.1.2. 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 componentsw, θx, θy as(31)w=∑i=18Niwi,θx=∑i=18Niθxi,θy=∑i=18Niθyi,where (32)Ni=141+ξiξ1+ηiηξiξ+ηiη-1,i=1,2,3,4Ni=121-ξ21+ηiη,i=5,7Ni=121-η21+ξiξ,i=6,8.Another 4-node shape functions are used to interpolate the deflection part of displacement function as(33)wb=∑i=14N-iwi,where (34)N-i=141+ξiξ1+ηiη,i=1,2,3,4.So the new deflection function partw∗ can be described as (35)w∗=w+αw-wb,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:(36)w∗θxθyT=w+αw-wbθxθyT=Nqand the element consistent mass matrix Mce can be expressed as(37)Mce=∬ΩNTmNdxdy,where m is the matrix containing the mass density of the material ρ and thickness t as follows:(38)m=ρt000t312000t312.


### 3.2. 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(39)K-λKGΦ=0,where K is the global stiffness matrix, KG 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 KGe can be written as(40)KGe=KG1e+KG2e+KG3e,where the terms of (40) are(41)KG1e=t∫Ωe∂w∂x∂w∂yσ0T∂w∂x∂w∂ydxdyKG2e=t312∫Ωe∂θx∂x∂θx∂yσ0T∂θx∂x∂θx∂ydxdyKG3e=t312∫Ωe∂θy∂x∂θy∂yσ0T∂θy∂x∂θy∂ydxdyand vector σ0=σx0τxy0τxy0σy0 stands for the initial stress of plate, σx0, σy0 are unit tractions in the direction of x and y, and τxy0 is the unit shear stress.Similarly, we can obtain the refined geometric stiffness matrix from the following relation:(42)∂w-∂x∂w-∂y=∂w∂x∂w∂y+α0∂w∂x∂w∂y-θy-θx+α0∂w∂x∂w∂y-∂wb∂x∂wb∂y,where α0 is an adjustable constant used to improve the accuracy of the stability analysis, the displacement interpolation functions w, θx, θy are the same as 8-node shape function (31), and wb is in terms of (33).Based on (42), we only modify the interpolation function w and geometric stiffness matrix part KG1e as(43)K-G1e=t∫Ωe∂w-∂x∂w-∂yσ0T∂w-∂x∂w-∂ydxdy.The refined geometric stiffness matrix can be written as(44)K-Ge=K-G1e+KG2e+KG3e.


## 3.1. Free Vibration Analysis

It is well-known that the equation of free vibration can be expressed as follows:(29)K-ω2MΦ=0,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 matrixM. Two types of the mass matrix are of interest: lumped mass matrix and consistent mass matrix.


### 3.1.1. 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: (30)Mlumpede=me36diag⁡100100100100800800800800,where me is the element mass.Figure 2
Mass lumping at the nodes.


### 3.1.2. 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 componentsw, θx, θy as(31)w=∑i=18Niwi,θx=∑i=18Niθxi,θy=∑i=18Niθyi,where (32)Ni=141+ξiξ1+ηiηξiξ+ηiη-1,i=1,2,3,4Ni=121-ξ21+ηiη,i=5,7Ni=121-η21+ξiξ,i=6,8.Another 4-node shape functions are used to interpolate the deflection part of displacement function as(33)wb=∑i=14N-iwi,where (34)N-i=141+ξiξ1+ηiη,i=1,2,3,4.So the new deflection function partw∗ can be described as (35)w∗=w+αw-wb,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:(36)w∗θxθyT=w+αw-wbθxθyT=Nqand the element consistent mass matrix Mce can be expressed as(37)Mce=∬ΩNTmNdxdy,where m is the matrix containing the mass density of the material ρ and thickness t as follows:(38)m=ρt000t312000t312.


## 3.1.1. 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: (30)Mlumpede=me36diag⁡100100100100800800800800,where me is the element mass.Figure 2
Mass lumping at the nodes.


## 3.1.2. 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 componentsw, θx, θy as(31)w=∑i=18Niwi,θx=∑i=18Niθxi,θy=∑i=18Niθyi,where (32)Ni=141+ξiξ1+ηiηξiξ+ηiη-1,i=1,2,3,4Ni=121-ξ21+ηiη,i=5,7Ni=121-η21+ξiξ,i=6,8.Another 4-node shape functions are used to interpolate the deflection part of displacement function as(33)wb=∑i=14N-iwi,where (34)N-i=141+ξiξ1+ηiη,i=1,2,3,4.So the new deflection function partw∗ can be described as (35)w∗=w+αw-wb,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:(36)w∗θxθyT=w+αw-wbθxθyT=Nqand the element consistent mass matrix Mce can be expressed as(37)Mce=∬ΩNTmNdxdy,where m is the matrix containing the mass density of the material ρ and thickness t as follows:(38)m=ρt000t312000t312.


## 3.2. 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(39)K-λKGΦ=0,where K is the global stiffness matrix, KG 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 KGe can be written as(40)KGe=KG1e+KG2e+KG3e,where the terms of (40) are(41)KG1e=t∫Ωe∂w∂x∂w∂yσ0T∂w∂x∂w∂ydxdyKG2e=t312∫Ωe∂θx∂x∂θx∂yσ0T∂θx∂x∂θx∂ydxdyKG3e=t312∫Ωe∂θy∂x∂θy∂yσ0T∂θy∂x∂θy∂ydxdyand vector σ0=σx0τxy0τxy0σy0 stands for the initial stress of plate, σx0, σy0 are unit tractions in the direction of x and y, and τxy0 is the unit shear stress.Similarly, we can obtain the refined geometric stiffness matrix from the following relation:(42)∂w-∂x∂w-∂y=∂w∂x∂w∂y+α0∂w∂x∂w∂y-θy-θx+α0∂w∂x∂w∂y-∂wb∂x∂wb∂y,where α0 is an adjustable constant used to improve the accuracy of the stability analysis, the displacement interpolation functions w, θx, θy are the same as 8-node shape function (31), and wb is in terms of (33).Based on (42), we only modify the interpolation function w and geometric stiffness matrix part KG1e as(43)K-G1e=t∫Ωe∂w-∂x∂w-∂yσ0T∂w-∂x∂w-∂ydxdy.The refined geometric stiffness matrix can be written as(44)K-Ge=K-G1e+KG2e+KG3e.


## 4. 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), soft-type 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 x=0, a free edge along y=0, a clamped edge along x=a, and a free edge along y=b.


### 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 Figure3. The size of the domain is 0.24×0.12. The mechanical properties of the plate are chosen as E=103, μ=0.25, κ=5/6.Figure 3
Constant stress patch test.For zero shear deformation patch test, the test function [20, 23] can be written as follows:(45)w=a0+a1x+a2y+a3x2+a4xy+a5y2θx=a2+a4x+2a5yθy=-a1-2a3x-a4y.For nonzero shear deformation patch test, the test function can be expressed as follows:(46)w=a0+a1x+a2y+a3x2+a4xy+a5y2+a6x3+a7x2y+a8xy2+a9y3θx=h23k1-μa7+3a9+a2+a4x+2a5y+a7x2+2a8xy+3a9y2θy=-h23k1-μ3a6+a8-a1-2a3x-a4y-3a6x2-2a7xy-a8y2,where ai, (i=0,…,9) are arbitrary constants. In the present paper, we assume that ai=i+1.The mesh of patch test is presented in Figure3. 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 Sx, Sy (Table 2) and the strict patch test with nonzero constant shear stresses Sx, Sy (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.Table 2
Numerical results of constant bending patch test (at node(0.04,0.02)).
w
θ
x
θ
y
M
x
M
y
M
x
y
S
x
S
y
t
=
0.05
Mesh A
1.1528
3.440
-
2.420
-
0.12222
-
0.15556
-
0.041667
0
0
Mesh B
1.1528
3.440
-
2.420
-
0.12222
-
0.15556
-
0.041667
0
0
Mesh C
1.1528
3.440
-
2.420
-
0.12222
-
0.15556
-
0.041667
0
0
Exact
1.1528
3.440
-
2.420
-
0.12222
-
0.15556
-
0.041667
0
0
t
=
0.005
Mesh A
1.1528
3.440
-
2.420
-
0.0001222
-
0.0001556
-
0.000041667
0
0
Mesh B
1.1528
3.440
-
2.420
-
0.0001222
-
0.0001556
-
0.000041667
0
0
Mesh C
1.1528
3.440
-
2.420
-
0.0001222
-
0.0001556
-
0.000041667
0
0
Exact
1.1528
3.440
-
2.420
-
0.0001222
-
0.0001556
-
0.000041667
0
0Table 3
Numerical results of nonzero constant shear stress patch test (at node(0.04,0.02)).
w
θ
x
θ
y
M
x
M
y
M
x
y
S
x
S
y
t
=
0.05
Mesh A
1.1537
3.529867
-
2.51
-
0.156667
-
0.196667
-
0.053
-
0.66667
-
0.84444
Mesh B
1.1537
3.529867
-
2.51
-
0.156667
-
0.196667
-
0.053
-
0.66667
-
0.84444
Mesh C
1.1537
3.529867
-
2.51
-
0.156667
-
0.196667
-
0.053
-
0.66667
-
0.84444
Exact
1.1537
3.529867
-
2.51
-
0.156667
-
0.196667
-
0.053
-
0.66667
-
0.84444
t
=
0.005
Mesh A
1.1537
3.479707
-
2.4704
-
0.0001567
-
0.0001967
-
0.000053
-
0.00067
-
0.000844
Mesh B
1.1537
3.479707
-
2.4704
-
0.0001567
-
0.0001967
-
0.000053
-
0.00067
-
0.000844
Mesh C
1.1537
3.479707
-
2.4704
-
0.0001567
-
0.0001967
-
0.000053
-
0.00067
-
0.000844
Exact
1.1537
3.479707
-
2.4704
-
0.0001567
-
0.0001967
-
0.000053
-
0.00067
-
0.000844


### 4.2. 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 modulusE=2.0×1011, 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/b2D/ρ are obtained. b is the width of the edge, D=Et3/12(1-μ2). Unless otherwise specified, the actual computations are done using dimensionless units.


#### 4.2.1. 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 a and the width of the plate b: a/b=1) is divided into mesh of 16×16 elements. The thickness-to-width aspect ratios t/b=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.Table 4
The six lowest nondimensional frequency parameters of a CCCC thin square plate(t/L=0.001).
Element (α)
Mode sequence number
1
2
3
4
5
6
QH8-RCV(-1)
3.6742
7.5726
7.5726
11.2791
13.8025
13.8626
QH8-RCV(-0.5)
3.6602
7.5052
7.5052
11.1242
13.5707
13.6325
QH8-RCV(-0.25)
3.6532
7.4711
7.4711
11.0458
13.4531
13.5157
QH8-RCV(-0.08)
3.6484
7.4478
7.4478
10.9922
13.3726
13.4356
QH8-RCV(0)
3.6461
7.4368
7.4368
10.9669
13.3345
13.3978
QH8-LV
3.6460
7.4363
7.4363
10.9644
13.3318
13.3949
Reference [27]
3.646
7.436
7.436
10.964
13.332
13.395Table 5
The six lowest nondimensional frequency parameters of a CCCC thick square plate(t/L=0.1).
Element (α)
Mode sequence number
1
2
3
4
5
6
QH8-RCV(-1)
3.3118
6.3402
6.3402
8.9003
10.4844
10.5761
QH8-RCV(-0.5)
3.3066
6.3332
6.3332
8.9090
10.5178
10.6160
QH8-RCV(-0.25)
3.3002
6.3038
6.3038
8.8443
10.4237
10.5216
QH8-RCV(-0.08)
3.2943
6.2743
6.2743
8.7752
10.3198
10.4162
QH8-RCV(0)
3.2912
6.2578
6.2578
8.7360
10.2604
10.3558
QH8-LV
3.3181
6.3688
6.3688
8.9622
10.5759
10.6670
Reference [27]
3.292
6.276
6.276
8.792
10.356
10.455Table 6
The six lowest nondimensional frequency parameters of a SSSS thin square plate(t/L=0.001).
Element (α)
Mode sequence number
1
2
3
4
5
6
QH8-RCV(-1)
2.0129
5.0805
5.0805
8.2063
10.3226
10.3226
QH8-RCV(-0.5)
2.0064
5.0403
5.0403
8.1034
10.1616
10.1616
QH8-RCV(-0.25)
2.0032
5.0202
5.0202
8.0522
10.0813
10.0813
QH8-RCV(-0.08)
2.0010
5.0066
5.0066
8.0174
10.0268
10.0268
QH8-RCV(0)
2.0000
5.0002
5.0002
8.0010
10.0012
10.0012
QH8-LV
2.0000
5.0000
5.0000
7.9998
9.9999
9.9999
Reference [27]
2.000
5.000
5.000
8.000
10.000
10.000Table 7
The six lowest nondimensional frequency parameters of a SSSS thick square plate(t/L=0.1).
Element (α)
Mode sequence number
1
2
3
4
5
6
QH8-RCV(-1)
1.9424
4.6599
4.6599
7.1742
8.7513
8.7513
QH8-RCV(-0.5)
1.9376
4.6420
4.6420
7.1505
8.7325
8.7325
QH8-RCV(-0.25)
1.9345
4.6229
4.6229
7.1032
8.6600
8.6600
QH8-RCV(-0.08)
1.9321
4.6062
4.6062
7.0581
8.5881
8.5881
QH8-RCV(0)
1.9308
4.5973
4.5973
7.0335
8.5482
8.5482
QH8-LV
1.9457
4.6784
4.6784
7.2173
8.8157
8.8157
Reference [27]
1.931
4.605
4.605
7.064
8.605
8.605Figure 4
Error of six lowest nondimensional frequency parameters.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 N×N uniform meshes with N=4,8,16. The results agree well with [27] when N=8,16.Table 8
Convergence of six lowest nondimensional frequency parameters (t/h=0.001, 0.1 and SSSS, CCCC).
t
/
L
Boundary
Mesh
Mode sequence number
1
2
3
4
5
6
(a) QH8 with modified consistent mass matrix (QH8-RCV,α=-0.08)
0.001
SSSS
4
×
4
2.0208
5.1600
5.1600
8.5457
10.7782
10.7782
8
×
8
2.0044
5.0290
5.0290
8.0832
10.1251
10.1251
16
×
16
2.0010
5.0066
5.0066
8.0174
10.0268
10.0268
Ref. sol
2.000
5.000
5.000
8.000
10.000
10.000
CCCC
4
×
4
3.6986
7.7478
7.7478
11.8207
14.7803
14.9043
8
×
8
3.6563
7.4895
7.4895
11.1031
13.5317
13.5943
16
×
16
3.6484
7.4478
7.4478
10.9922
13.3726
13.4356
Ref. sol
3.646
7.436
7.436
10.964
13.332
13.395
0.1
SSSS
4
×
4
1.9440
4.6822
4.6822
7.3572
8.9607
8.9607
8
×
8
1.9335
4.6062
4.6062
7.0516
8.5663
8.5663
16
×
16
1.9321
4.6062
4.6062
7.0581
8.5881
8.5881
Ref. sol
1.931
4.605
4.605
7.064
8.605
8.605
CCCC
4
×
4
3.3008
6.3401
6.3401
9.0667
10.7861
10.8677
8
×
8
3.2925
6.2555
6.2555
8.7327
10.2492
10.3400
16
×
16
3.2943
6.2743
6.2743
8.7752
10.3198
10.4162
Ref. sol
3.292
6.276
6.276
8.792
10.356
10.455
(b) QH8 with lumped mass matrix (QH8-LV)
0.001
SSSS
4
×
4
1.9997
4.9970
4.9970
7.9691
9.9694
9.9694
8
×
8
2.0000
4.9999
4.9999
7.9986
10.0001
10.0001
16
×
16
2.0000
5.0000
5.0000
7.9998
9.9999
9.9999
Ref. sol
2.000
5.000
5.000
8.000
10.000
10.000
CCCC
4
×
4
3.6452
7.4121
7.4121
10.8713
12.9885
13.0468
8
×
8
3.6461
7.4365
7.4365
10.9637
13.3312
13.3935
16
×
16
3.6460
7.4363
7.4363
10.9644
13.3318
13.3949
Ref. sol
3.646
7.436
7.436
10.964
13.332
13.395
0.1
SSSS
4
×
4
1.9427
4.6385
4.6385
7.0569
8.5511
8.5511
8
×
8
1.9452
4.6710
4.6710
7.1895
8.7670
8.7670
16
×
16
1.9457
4.6784
4.6784
7.2173
8.8157
8.8157
Ref. sol
1.931
4.605
4.605
7.064
8.605
8.605
CCCC
4
×
4
3.3019
6.2524
6.2524
8.6322
10.0169
10.0917
8
×
8
3.3155
6.3496
6.3496
8.9075
10.4886
10.5773
16
×
16
3.3181
6.3688
6.3688
8.9622
10.5759
10.6670
Ref. sol
3.292
6.276
6.276
8.792
10.356
10.455Figure 5
Convergence of six lowest nondimensional frequency parameters (t/h=0.001, 0.1, and SSSS, CCCC).Then, natural frequencies analyses are further carried on the thin(t/b=0.001) and thick (t/b=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.Table 9
Frequency parameters for square Mindlin plates with different boundary conditions.
Boundary
t
/
b
Element
Mode sequence number
1
2
3
4
5
6
SSFF
0.001
QH8-RCV
0.3407
1.7530
1.9541
3.8675
5.1729
5.4231
QH8-LV
0.3405
1.7485
1.9504
3.8550
5.1471
5.4000
Ref. sol
0.341
1.754
1.955
3.872
5.171
5.419
0.1
QH8-RCV
0.3329
1.6780
1.8751
3.5579
4.7179
4.9445
QH8-LV
0.3335
1.6944
1.8922
3.6122
4.8227
5.0496
Ref. sol
0.333
1.677
1.874
3.557
4.718
4.945
SFSF
0.001
QH8-RCV
0.9759
1.6329
3.7170
3.9486
4.7357
7.1656
QH8-LV
0.9755
1.6303
3.7009
3.9431
4.7240
7.1313
Ref. sol
0.976
1.635
3.721
3.946
4.735
7.167
0.1
QH8-RCV
0.9567
1.5594
3.4297
3.6829
4.3331
6.2859
QH8-LV
0.9601
1.5703
3.4900
3.7305
4.4011
6.4372
Ref. sol
0.956
1.559
3.429
3.681
4.332
6.290
SCSC
0.001
QH8-RCV
2.9350
5.5542
7.0348
9.6060
10.3849
13.1192
QH8-LV
2.9333
5.5466
7.0242
9.5833
10.3566
13.0799
Ref. sol
2.933
5.547
7.024
9.583
10.356
13.080
0.1
QH8-RCV
2.7019
4.9725
5.9898
7.9621
8.7685
10.2141
QH8-LV
2.7212
5.0501
6.0804
8.1367
8.9997
10.4651
Ref. sol
2.700
4.971
5.990
7.972
8.787
10.249
CFFF
0.001
QH8-RCV
0.3516
0.8603
2.1563
2.7540
3.1315
5.4804
QH8-LV
0.3514
0.8590
2.1499
2.7427
3.1211
5.4505
Ref. sol
0.352
0.862
2.157
2.756
3.136
5.490
0.1
QH8-RCV
0.3476
0.8167
2.0353
2.5832
2.8610
4.8108
QH8-LV
0.3481
0.8205
2.0567
2.6220
2.8963
4.9106
Ref. sol
0.348
0.816
2.034
2.582
2.860
4.811
CCFF
0.001
QH8-RCV
0.7003
2.4197
2.6929
4.8226
6.3563
6.6444
QH8-LV
0.6997
2.4117
2.6870
4.8039
6.3216
6.6129
Ref. sol
0.701
2.422
2.694
4.828
6.353
6.640
0.1
QH8-RCV
0.6761
2.2434
2.5046
4.2523
5.5551
5.8093
QH8-LV
0.6777
2.2687
2.5294
4.3213
5.6786
5.9314
Ref. sol
0.676
2.242
2.503
4.251
5.557
5.812
CFCF
0.001
QH8-RCV
2.2459
2.6732
4.4127
6.2028
6.8074
8.0892
QH8-LV
2.2448
2.6689
4.3935
6.1928
6.7890
8.0240
Ref. sol
2.246
2.675
4.417
6.198
6.806
8.087
0.1
QH8-RCV
2.0899
2.4334
3.9032
5.3315
5.7707
6.9214
QH8-LV
2.0979
2.4508
3.9716
5.3982
5.8585
7.1420
Ref. sol
2.088
2.431
3.901
5.331
5.771
6.929Figure 6
Shape of the six lowest modes of SSSS thin square plate obtained from QH8-RCV (α=-0.08, t/h=0.001).


#### 4.2.2. Free Vibration of Circular Plate

Free vibration of a circular plate (Figure7) of R=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 E=2.0×1011, and Poisson’s ratio μ=0.3.Figure 7
Circular plate model and its discretization for QH8 element (Mesh A: 205 nodes, 60 elements; Mesh B: 321 nodes, 96 elements).60 proposed elements and 205 nodes as shown in Figure7, 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 (t/2R=0.01) and thick (t/2R=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.Table 10
Ten lowest parameterized natural frequenciesω-=ωR2ρt/D of a clamped circular plate with t/2R=0.01.
Mode
NS-ES-FEM
S8R
MIN3
ES-FEM+DSG3
ANS4
QH8-RCV
QH8-LV
Ref. sol [42]
1
10.2700
10.3095
10.4082
10.3109
10.2572
10.3613
10.3669
10.216
2
21.1535
21.4233
22.2198
21.6702
21.4981
21.6355
21.5152
21.260
3
21.1617
21.4283
22.2444
21.6900
21.4981
21.6356
21.5167
21.260
4
34.4423
35.2809
37.7461
36.3124
35.3941
35.6621
35.1192
34.880
5
34.4573
35.3685
37.7816
36.3816
35.5173
35.7071
35.2812
34.880
6
39.1472
40.1039
43.0344
41.3801
40.8975
40.7676
40.1565
39.771
7
49.6819
52.2399
57.8881
54.7796
52.2054
52.6337
51.3212
51.040
8
49.7157
52.2597
58.0836
54.8922
52.2054
52.6344
51.3358
51.040
9
58.8517
62.2179
68.7260
64.6300
63.2397
63.0358
61.0974
60.820
10
59.0184
62.2543
69.2354
65.1330
63.2397
63.0368
61.1017
60.820Table 11
Ten lowest parameterized natural frequenciesω-=ωR2ρt/D of a clamped circular plate with t/2R=0.1.
Mode
NS-ES-FEM
S8R
MIN3
ES-FEM+DSG3
ANS4
QH8-RCV
QH8-LV
Ref. sol [43]
1
9.2460
9.2539
9.9682
9.3262
9.2605
9.3380
9.4444
9.240
2
17.6424
17.8115
20.306
18.0461
17.9469
17.8602
18.1926
17.834
3
17.6461
17.8115
20.3256
18.0673
17.9469
17.8603
18.1938
17.834
4
26.5699
27.0923
32.7800
27.8438
27.0345
26.9746
27.5602
27.214
5
26.5740
27.1650
32.7905
27.8856
27.6566
26.9836
27.6854
27.214
6
29.7654
30.4101
37.1591
31.1280
30.3221
30.0846
30.8631
30.211
7
35.5595
33.8815
47.5437
38.6936
37.2579
36.5419
37.3990
37.109
8
35.5999
33.8848
47.5993
38.7042
37.2579
36.5431
37.4096
37.109
9
41.0560
37.0290
55.3530
44.1130
43.2702
41.8999
42.8445
42.409
10
41.1226
37.0324
55.6608
44.3286
43.2702
41.9011
42.8465
42.409Figure 8
Ten lowest parameterized natural frequencies of a clamped circular plate.
(a)
Thin circular plate(t/2R=0.01)
(b)
Thick circular plate(t/2R=0.1)First six lowest modes of this clamped circular obtained from QH8-RCV using the mesh of Figure7 Mesh B are shown in Figure 9, which exactly describes the real physical modes of the circular plate.Figure 9
Shape of the six lowest modes of a clamped thick circular plate obtained from QH8-RCV(t/2R=0.1).


#### 4.2.3. 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(a). The material parameters of the triangular plate are given as follows: the density mass ρ=8000, Young’s modulus E=2.0×1011, and Poisson’s ratio μ=0.3. The plate is discretized by the meshes as shown in Figure 10 for QH8. The convergence study of nondimensional frequency parameter ω-=ωL2ρt/D/π2 for thin (t/L=0.001) and thick (t/L=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.Table 12
The lowest six parameterized natural frequencies of triangular plates witht/L=0.001.
Skew angle
Method
Mode sequence number
1
2
3
4
5
6
0°
DSG3
0.6288
2.4399
3.4796
5.9913
8.3737
11.2081
MIN3
0.6272
2.4242
3.4344
5.9458
8.2460
11.0178
ES-FEM+DSG3
0.6241
2.3970
3.3631
5.7976
7.9696
10.5356
NS+ES-FEM
0.6253
2.3335
3.1888
5.4772
7.4350
9.2790
Rayleigh-Ritz
0.6240
2.3770
3.3080
5.6890
7.7430
—
ASN4
0.6240
2.3790
3.3170
5.7240
7.7940
10.2000
QH8-RCV
0.6232
2.3709
3.2980
5.6766
7.7627
10.0574
QH8-LV
0.6165
2.3143
3.1843
5.4810
7.5749
9.6271
15°
DSG3
0.5915
2.2360
3.6467
5.5686
7.8875
11.1572
MIN3
0.5886
2.2229
3.5717
5.5286
7.7447
10.9544
ES-FEM+DSG3
0.5843
2.1979
3.4937
5.3893
7.4745
10.4532
NS+ES-FEM
0.5854
2.1449
3.2938
5.0391
7.0313
9.1652
Rayleigh-Ritz
0.5840
2.1810
3.4090
5.2800
7.2640
—
ASN4
0.5830
2.1810
3.4130
5.3030
7.2890
10.0950
QH8-RCV
0.5829
2.1759
3.3995
5.2687
7.2983
9.9691
QH8-LV
0.5750
2.1291
3.2332
5.0577
7.0972
9.4300
30°
DSG3
0.5928
2.2293
3.9729
5.8523
7.8772
11.9930
MIN3
0.5866
2.2164
3.8656
5.8059
7.6697
11.5791
ES-FEM+DSG3
0.5795
2.1881
3.7620
5.6464
7.3608
10.9609
NS+ES-FEM
0.5801
2.1433
3.5402
5.2019
6.8932
9.6373
Rayleigh-Ritz
0.5760
2.1740
3.6390
5.5110
7.1080
—
ASN4
0.5750
2.1740
3.6380
5.5340
7.1390
10.4770
QH8-RCV
0.5763
2.1685
3.6379
5.4957
7.1296
10.3946
QH8-LV
0.5714
2.1368
3.5544
5.3345
6.9689
9.9234
45°
DSG3
0.6265
2.4298
4.5709
6.9888
8.7055
13.4984
MIN3
0.6125
2.4038
4.4237
6.9194
8.3447
12.7848
ES-FEM+DSG3
0.6001
2.3557
4.2652
6.6442
7.9210
11.8923
NS+ES-FEM
0.6020
2.3189
4.0352
6.0299
7.3322
10.5988
Rayleigh-Ritz
0.5900
2.3290
4.1370
6.3810
7.6020
—
ASN4
0.5880
2.3240
4.1260
6.3810
7.6140
11.2240
QH8-RCV
0.5925
2.3290
4.1374
6.4133
7.6019
11.2205
QH8-LV
0.5858
2.3042
3.9784
6.1554
7.4240
10.8429
60°
DSG3
0.6969
2.9282
6.0785
9.3379
11.8340
17.8590
MIN3
0.6631
2.8106
5.8343
8.9727
11.3197
16.3597
ES-FEM+DSG3
0.6371
2.6804
5.4909
8.1353
10.6516
14.5225
NS+ES-FEM
0.6546
2.6792
5.3330
7.3384
9.6825
12.9900
Rayleigh-Ritz
0.6170
2.5760
5.3760
7.5240
10.2850
—
ASN4
0.6130
2.5640
5.3530
7.4600
10.3060
12.9420
QH8-RCV
0.6234
2.5978
5.3502
7.5829
10.2623
13.0926
QH8-LV
0.6164
2.5712
5.1891
7.3141
9.8462
12.8085Table 13
The lowest six parameterized natural frequencies of triangular plates witht/L=0.2.
Skew angle
Method
Mode sequence number
1
2
3
4
5
6
0°
DSG3
0.5872
1.9529
2.4678
4.1305
5.2243
6.2886
MIN3
0.6042
2.1109
2.7737
4.6581
6.1311
7.3208
ES-FEM+DSG3
0.5856
1.9376
2.4411
4.0678
5.1289
6.1359
Rayleigh-Ritz
0.5820
1.9000
2.4080
3.9360
—
—
ASN4
0.5820
1.9150
2.4280
3.9840
5.0180
5.9440
NS+ES-FEM
0.5821
1.8710
2.3600
3.8330
4.8217
5.5683
QH8-RCV
0.5821
1.8907
2.3931
3.8677
4.8701
5.6897
QH8-LV
0.5826
1.9628
2.5320
4.1115
5.2710
6.4029
15°
DSG3
0.5478
1.8160
2.4429
3.8111
5.0253
6.0033
MIN3
0.5656
1.9555
2.7772
4.2851
5.8857
6.9977
ES-FEM+DSG3
0.5457
1.8018
2.4121
3.7529
4.9317
5.8327
Rayleigh-Ritz
0.5440
1.7710
2.3860
3.6280
—
—
ASN4
0.5450
1.7640
2.4200
3.6080
4.8200
5.4310
NS+ES-FEM
0.5442
1.7466
2.3391
3.5289
4.6721
5.1901
QH8-RCV
0.5441
1.7644
2.3687
3.5813
4.7024
5.3504
QH8-LV
0.5427
1.8251
2.5244
3.8353
5.0144
6.1823
30°
DSG3
0.5365
1.8126
2.4954
3.7549
4.9796
5.9219
MIN3
0.5582
1.9551
2.8595
4.2605
5.8273
6.9945
ES-FEM+DSG3
0.5329
1.7963
2.4552
3.6902
4.8717
5.6942
Rayleigh-Ritz
0.5330
1.7720
2.4190
3.5650
—
—
ASN4
0.5320
1.7730
2.4370
3.5910
4.7650
5.3230
NS+ES-FEM
0.5336
1.7485
2.3730
3.4615
4.5927
4.9630
QH8-RCV
0.5331
1.7663
2.4045
3.5167
4.6384
5.1512
QH8-LV
0.5342
1.8391
2.7165
3.9146
4.9778
6.4454
45°
DSG3
0.5448
1.9269
2.6246
3.9196
5.1474
5.9838
MIN3
0.5724
2.1044
3.0165
4.5535
6.0525
7.2736
ES-FEM+DSG3
0.5380
1.9002
2.5682
3.8239
4.9806
5.6000
Rayleigh-Ritz
0.5400
1.8850
2.4890
3.6740
—
—
ASN4
0.5410
1.8840
2.5180
3.7480
4.7400
5.2920
NS+ES-FEM
0.5408
1.8634
2.4323
3.5593
4.4299
4.9328
QH8-RCV
0.5391
1.8765
2.4729
3.6337
4.5699
5.0583
QH8-LV
0.5389
1.9562
2.9904
4.2186
5.2159
6.9240
60°
DSG3
0.5691
2.1281
2.7093
4.3427
5.3437
6.5278
MIN3
0.6017
2.3742
3.2377
5.1762
6.5578
8.1003
ES-FEM+DSG3
0.5569
2.0692
2.5955
4.1109
4.7967
6.0192
Rayleigh-Ritz
0.5590
2.0590
2.3960
3.5900
—
—
ASN4
0.5590
2.0950
2.4830
3.9100
4.5170
5.7630
NS+ES-FEM
0.5611
2.0446
2.3276
3.5969
4.2312
5.3615
QH8-RCV
0.5574
2.0455
2.3807
3.6046
4.2131
5.2640
QH8-LV
0.5573
2.1170
3.8244
4.6656
6.2099
7.7794Figure 10
Triangular plate model with various shape geometries and its discretization meshes (θ=0°, 15°, 30°, 45°, 60°).
(a)
Triangular plate
(b)
0°, 28 elements and 106 nodes
(c)
15°, 26 elements and 103 nodes
(d)
30°, 27 elements and 106 nodes
(e)
45°, 26 elements and 105 nodes
(f)
60°, 28 elements and 107 nodesFigure 11
Convergence study of the lowest four parameterized natural frequencies of triangular plates witht/L=0.001.Figure 12
Convergence study of the lowest four parameterized natural frequencies of triangular plates witht/L=0.2.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.Figure 13
The six lowest modes of a CFF thin triangular plate obtained from QH8-RCV (θ=15°).


### 4.3. 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.(47)Kb=λb2π2D,with b being the edge width of the plate that is used throughout the paper and D=Et3/12(1-μ2).


#### 4.3.1. Rectangular Plates Subjected to an In-Plane Compressive Load

The parameterα0 in (42) for the present element QH8-R (8-node hybrid element with refined geometric stiffness matrix) is determined by comparison using numerical examples. A square plate (a/b=1) compressed in one direction or two perpendicular directions are divided into different meshes 2×2, 4×4, 6×6, 8×8, 12×12, 16×16. Thickness-to-width aspect ratios t/b=0.01, 0.1, α0=-1,-0.75,-0.5,-0.3,-0.17,0, and boundary conditions SSSS, CCCC are chosen. The computational results of the critical load are given in Tables 14 and 15 and Figure 14, in which the exact solutions can be found in [44–46]. After analysis of numerical examples, the present 8-node hybrid stress element QH8-R(α0) has higher accuracy when -0.4≤α0≤-0.1, and we advise α0=-0.17 for following numerical examples.Table 14
The factor of uniaxial buckling loads along thex-axis of rectangular plates with length-to-width aspect ratios a/b=1 and thickness-to-width aspect ratios t/b=0.01.
Boundary condition
Element (α0)
Mesh
2
×
2
4
×
4
6
×
6
8
×
8
12
×
12
16
×
16
SSSS
QH8-R(-1)
3.1427
3.4982
3.7191
3.8183
3.8982
3.9293
(
-
0.75
)
4.0364
3.8502
3.9045
3.9338
3.9588
3.9691
(
-
0.5
)
4.7314
4.0776
4.0190
4.0043
3.9961
3.9940
(
-
0.3
)
4.8659
4.1389
4.0521
4.0258
4.0084
4.0027
(
-
0.17
)
4.7097
4.1161
4.0444
4.0225
4.0079
4.0031
(
0
)
4.2694
4.0139
4.0002
3.9980
3.9972
3.9971
Reference [44]
4
CCCC
QH8-R(-1)
5.3071
7.4051
8.3520
8.9070
9.4027
9.6049
(
-
0.75
)
8.2670
9.1692
9.4445
9.6292
9.7956
9.8651
(
-
0.5
)
10.9935
10.6774
10.2199
10.1074
10.0459
10.0304
(
-
0.3
)
12.9592
11.1920
10.4649
10.2587
10.1293
10.0885
(
-
0.17
)
13.6765
11.0371
10.4113
10.2349
10.1247
10.0896
(
0
)
13.6316
10.2573
10.0875
10.0576
10.0482
10.0472
Reference [44]
10.0738
Note. κ=5/6.Table 15
The factor of uniaxial buckling loads along thex-axis of rectangular plates with length-to-width aspect ratios a/b=1 and thickness-to-width aspect ratios t/b=0.1.
Boundary
Element (α0)
Mesh
2
×
2
4
×
4
6
×
6
8
×
8
12
×
12
16
×
16
SSSS
QH8-R(-1)
2.3811
2.4760
2.7861
3.0392
3.3424
3.4914
(
-
0.75
)
3.3240
3.0601
3.2413
3.3983
3.5815
3.6689
(
-
0.5
)
4.1951
3.5827
3.6110
3.6676
3.7403
3.7756
(
-
0.3
)
4.4995
3.8245
3.7755
3.7778
3.7917
3.8002
(
-
0.17
)
4.3933
3.8501
3.7970
3.7869
3.7845
3.7850
(
0
)
3.9119
3.6966
3.7056
3.7146
3.7232
3.7266
Reference [45]
3.7864
CCCC
QH8-R(-1)
3.4293
4.0301
4.6232
5.2042
6.1652
6.7926
(
-
0.75
)
5.9999
5.8082
6.0193
6.4602
7.1856
7.6362
(
-
0.5
)
8.7108
7.7026
7.4918
7.6519
7.9938
8.2101
(
-
0.3
)
10.4843
8.8128
8.3089
8.2536
8.3146
8.3741
(
-
0.17
)
10.9492
8.8625
8.3903
8.2963
8.2811
8.2936
(
0
)
9.9368
7.5669
7.5101
7.5613
7.6923
7.7874
Reference [46]
8.2917
Note. κ=5/6.Figure 14
The factor of uniaxial buckling loads (σx0) along the x-axis of rectangular plates for QH8-R (α0) 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 are investigated. Several thickness-span aspect ratiost/b=0.1, 0.2 (moderately thick plate) and t/b=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 nodes) and MK, Shear-MK use a regular pattern of 21×21 elements (441 nodes). A shear correction factor κ=π2/12 is used here for a proper comparison.Table 16
The factor of uniaxial buckling loads (σx0) along the x-axis of rectangular plates (a/b=1) with various boundary conditions and various thickness-to-width ratios (t/b).
Boundary
Method
t
/
b
0.001
0.05
0.1
0.2
SSSS
SSM
4.000
3.944
3.784
3.256
Shear-MK
4.0243
3.9153
3.7895
3.2748
MK
5.467
5.0111
3.7905
3.2799
Exact [45]
4.000
3.9437
3.7838
3.2558
S8R
3.997
3.9445
3.7866
1.7025
QH8-R
4.0132
3.9478
3.7827
3.2464
SFSC
SSM
1.653
1.620
1.556
1.370
Shear-MK
1.6552
1.6211
1.5658
1.3765
MK
2.5876
2.1778
1.5699
1.3612
Exact [45]
1.6522
1.6197
1.5558
1.3701
S8R
1.6560
1.6202
1.5570
1.2532
QH8-R
1.6522
1.6209
1.5549
1.3669
SFSS
SSM
1.402
1.378
1.327
1.173
Shear-MK
1.4211
1.3877
1.3795
1.2700
MK
2.2456
2.0323
1.3791
1.2755
Exact [45]
1.4014
1.3813
1.3707
1.2138
S8R
1.4017
1.3816
1.3415
1.2158
QH8-R
1.4022
1.3827
1.3410
1.2133
SFSF
SSM
0.9523
0.9412
0.9146
0.8274
Shear-MK
0.9676
0.9455
0.9231
0.8575
MK
1.7644
1.6578
0.9344
0.8654
Exact [45]
0.95225
0.94314
0.92187
0.85011
S8R
0.9524
0.9433
0.9223
0.8325
QH8-R
0.9532
0.9442
0.9227
0.8519
Note. QH8-R (α0=-0.17). Shear correction factor κ=π2/12.Figure 15
Normalized uniaxial buckling load factors along thex-axis (σx0) versus the thickness-span aspect ratio t/b=0.001,0.05,0.1,0.2 with various boundaries.Through the achieved results shown in Table16 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, t/b=0.001, 0.05, but conversely matches well with thick plates, that is, t/b=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.The results of uniaxial buckling loads(σy0) along they-axis and the biaxial in-plane compression load (σx0=σy0) are shown in Tables 17 and 18 and Figures 16 and 17 associated with five different densities of regular elements such as 2×2 (21 nodes), 4×4 (65 nodes), 6×6 (133 nodes), 8×8 (225 nodes), and 10×10 (341 nodes).Table 17
The factor of uniaxial buckling loads (σy0) along the y-axis of rectangular plates (a/b=1) with various boundary conditions and various thickness-to-width ratios (t/b).
Boundary
Element
t
/
b
0.001
0.05
0.1
0.2
SSSS
2
×
2
4.7132
4.6261
4.3905
3.6771
4
×
4
4.1193
4.0431
3.8479
3.2455
6
×
6
4.0476
3.9751
3.7949
3.2231
8
×
8
4.0257
3.9564
3.785
3.2321
10
×
10
4.0161
3.9495
3.7829
3.2422
Exact [45]
4.0000
3.9437
3.7838
3.2558
SFSC
2
×
2
2.5834
2.5216
2.3678
1.9583
4
×
4
2.3993
2.3194
2.1530
1.7580
6
×
6
2.3821
2.2898
2.1191
1.7288
8
×
8
2.3797
2.2794
2.1079
1.7213
10
×
10
2.3799
2.2743
2.1031
1.7190
Exact [45]
2.3901
2.2667
2.1010
1.7200
SFSS
2
×
2
2.5626
2.502
2.3511
1.9480
4
×
4
2.3741
2.2965
2.1346
1.7483
6
×
6
2.3567
2.2672
2.1010
1.7195
8
×
8
2.3541
2.2569
2.0901
1.7122
10
×
10
2.3542
2.2518
2.0853
1.7099
Exact [45]
2.3639
2.2442
2.0829
1.7105
SFSF
2
×
2
2.1000
2.058
1.9553
1.6775
4
×
4
2.0327
1.9737
1.8516
1.5567
6
×
6
2.0286
1.9589
1.8315
1.5358
8
×
8
2.0295
1.9531
1.8244
1.5301
10
×
10
2.0309
1.9500
1.8212
1.5282
Exact [45]
2.0413
1.9457
1.8216
1.5333
Note. QH8-R (α0=-0.17). Shear correction factor κ=π2/12.Table 18
The factor of biaxial buckling loads (σx0=σy0) for rectangular plates (a/b=1) with various boundary conditions and various thickness-to-width ratios (t/b).
Boundary
Element
t
/
b
0.001
0.005
0.05
0.1
SSSS
2
×
2
2.3774
2.3770
2.3377
2.2277
4
×
4
2.0671
2.0667
2.0334
1.9432
6
×
6
2.0275
2.0272
1.9954
1.9108
8
×
8
2.015
2.0147
1.984
1.9023
10
×
10
2.0095
2.0091
1.9793
1.8992
Exact [45]a
2.0000
1.9997
1.9718
1.8919
SFSC
2
×
2
1.1820
1.1818
1.1652
1.1205
4
×
4
1.1474
1.1472
1.1261
1.0767
6
×
6
1.1427
1.1424
1.1184
1.0679
8
×
8
1.1416
1.1413
1.1155
1.0649
10
×
10
1.1414
1.1410
1.1140
1.0636
Exact [45]b
1.1431
1.1412
1.1119
1.0641
SFSS
2
×
2
1.0925
1.0923
1.0794
1.0439
4
×
4
1.0606
1.0605
1.0445
1.0060
6
×
6
1.0558
1.0556
1.0378
0.9986
8
×
8
1.0546
1.0543
1.0352
0.9962
10
×
10
1.0541
1.0538
1.0340
0.9952
Exact [45]b
1.0548
1.0535
1.0322
0.99541
SFSF
2
×
2
0.9671
0.9670
0.9583
0.9334
4
×
4
0.9400
0.9399
0.9302
0.9053
6
×
6
0.9351
0.9349
0.9247
0.9002
8
×
8
0.9335
0.9333
0.9228
0.8986
10
×
10
0.9328
0.9326
0.9219
0.8979
Exact [45]b
0.93209
0.9316
0.92071
0.89774
Note. QH8-R (α0=-0.17). aShear correction factor κ=0.823045. bShear correction factor κ=5/6.Figure 16
Normalized axial buckling load factors along they-axis (σy0) versus the thickness-span aspect ratio t/b=0.001, 0.05, 0.1, 0.2 with various boundaries.Figure 17
Normalized biaxial buckling load factors (σx0=σy0) versus the thickness-span aspect ratio t/b=0.001, 0.05, 0.1, 0.2 with various boundaries.We also consider simply supported plates (SSSS) with various thickness-to-width ratios,t/b = 0.05, 0.1, 0.2, and length-to-width ratios, a/b = 0.5, 1.0, 1.5, 2.0, 2.5.Table19 and Figure 18 show the buckling factors using the element QH8-R with regular mesh of 9×9 rectangular elements (280 nodes) when a/b=1.0. The QH8-R results are also compared to the DSG3 and ES-DSG3 [47] (mesh of 16×16 elements and 289 nodes), S8R (mesh of 9×9 elements and 280 nodes), mesh-free method [48], and the pb-2 Ritz [46]. It is seen that the QH8-R exhibits a good agreement with the pb-2 Ritz method, and S8R is not suitable for plate t/b = 0.2. Figure 19 also depicts the uniaxial buckling modes of simply supported rectangular plates with thickness-to- width ratios t/b = 0.01 and various length-to-width ratios, a/b = 1.0, 1.5, 2.0, 2.5.Table 19
The factor of uniaxial buckling loads along thex-axis of rectangular plates (SSSS) with various length-to-width ratios and various thickness-to-width ratios.
a
/
b
t
/
b
S8R
ES-DSG3 [47]
Mesh-free [48]
QH8-R
pb-2 Ritz [46]
0.5
0.05
6.0376
5.9873
6.0405
6.0822
6.0372
0.1
5.4779
5.3064
5.3116
5.5071
5.4777
0.2
1.9525
3.720
3.7157
3.9754
3.9963
1
0.05
3.9441
3.9412
3.9293
3.9528
3.9444
0.1
3.7861
3.7402
3.7270
3.7855
3.7865
0.2
2.0444
3.1263
3.1471
3.2427
3.2637
1.5
0.05
4.2576
4.2852
4.2116
4.2745
4.2570
0.1
4.0255
3.9844
3.8982
4.0332
4.0250
0.2
2.0132
3.1461
3.1032
3.2966
3.3048
2
0.05
3.9440
3.9811
3.8657
3.9528
3.9444
0.1
3.7860
3.7711
3.6797
3.7855
3.7865
0.2
1.9486
3.1415
3.0783
3.2427
3.2637
2.5
0.05
4.0673
4.1691
3.9600
4.0795
4.0645
0.1
3.8700
3.8924
3.7311
3.8737
3.8683
0.2
1.8697
3.1234
3.0306
3.2338
3.2421Figure 18
The factor of axial buckling loads of SSSS plates with various length-to-width ratios (a/b) and various thickness-to-width ratios (t/b).Figure 19
Axial buckling modes of simply supported rectangular plates with thickness-to-width ratiost/b = 0.01 and various length-to-width ratios a/b.
(a)
a
/
b
=
1.0
(b)
a
/
b
=
1.5
(c)
a
/
b
=
2.0
(d)
a
/
b
=
2.5


#### 4.3.2. 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, t/b, 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.Figure 20
Skew plate under uniaxial in-plane loading.Numerical results for thick plates with thickness-to-width ratios,t/b = 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.Table 20
Buckling load intensity factor of skew plates havingS∗S∗S∗S∗, CCCC, CFCF, and S∗FS∗F boundary conditions subject to uniaxial loads (σx0).
Boundary
t
/
b
Skew angle
Mesh-free [48]
6
×
6
9
×
9
16
×
16
pb-2 Ritz [46]
S
∗
S
∗
S
∗
S
∗
0.05
0°
—
3.8057
3.7905
3.7815
3.7835
15°
—
4.1583
4.1380
4.1276
4.1314
30°
5.3874
5.4755
5.4304
5.4109
5.4182
45°
8.6926
8.9757
8.8002
8.7247
8.7328
0.1
0°
—
3.5111
3.5007
3.4965
3.4950
15°
—
3.8234
3.8094
3.8047
3.8027
30°
4.8534
4.9813
4.9475
4.9370
4.9324
45°
7.6106
7.9388
7.7931
7.7430
7.7236
CCCC
0.05
0°
9.5418
9.8091
9.6238
9.5647
9.5588
15°
10.2123
10.5549
10.3209
10.2430
10.2312
30°
12.5225
13.2413
12.7803
12.6114
12.5711
45°
17.8211
19.9900
18.6615
18.1211
17.9652
0.1
0°
8.1849
8.3903
8.2836
8.2936
8.2917
15°
8.7560
8.9338
8.7863
8.7845
8.7741
30°
10.3273
10.8235
10.4883
10.4238
10.3760
45°
13.2881
15.2116
14.2092
13.8638
13.6909
CFCF
0.05
0°
3.7984
3.8724
3.8281
3.8082
3.8007
15°
4.1310
4.2202
4.1688
4.1460
4.1387
30°
5.3510
5.4976
5.4102
5.3720
5.3660
45°
7.3658
7.9413
7.6341
7.4986
7.4670
0.1
0°
3.5063
3.5523
3.5262
3.5185
3.5077
15°
3.7684
3.8441
3.8135
3.8049
3.7937
30°
4.7773
4.8862
4.8316
4.8168
4.8043
45°
6.3092
6.6886
6.4547
6.3680
6.3311
S
∗
F
S
∗
F
0.05
0°
—
0.9459
0.9436
0.9424
0.9421
15°
—
1.0535
1.0513
1.0501
1.0501
30°
—
1.4612
1.4594
1.4583
1.4603
45°
—
2.5333
2.5317
2.5311
2.5450
0.1
0°
—
0.9227
0.9210
0.9202
0.9199
15°
—
1.0215
1.0199
1.0191
1.0189
30°
—
1.3892
1.3878
1.3873
1.3887
45°
—
2.3192
2.3187
2.3190
2.3316
Note. 6×6, 9×9, 16×16: the meshes for element QH8-R(-0.17).


#### 4.3.3. 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 lengtha=b=10 and the thickness t=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 x-direction, σx0, pure shear load, and biaxial compressive loads in x and y directions, where σx0=σy0, respectively.Figure 21
Square plate with a circular hole and the distributions (308 nodes, 84 elements).The results of the normalized buckling load factor calculated for those loads with different boundary conditions are shown in Table21. The results agree very well with those computed by the EFG [57] involving both the third-order 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.Table 21
The uniaxial, biaxial, and pure shear buckling load factor for a square plate with a circular hole and different boundary conditions.
Buckling load
Boundary
EFG (FSDT)
EFG (TSDT)
Shear-MK
HQ8-R
Uniaxial
CCCC
7.995
8.097
8.0777
8.0543
SSSS
1.986
1.969
1.9544
2.0580
Biaxial
CCCC
4.781
4.860
4.9878
4.8477
SSSS
1.032
1.021
1.0511
1.0696
Pure shear
CCCC
12.669
13.357
13.5252
12.6976
SSSS
7.867
7.873
8.1544
7.8971Figure 22
The first uniaxial, biaxial, and pear shear buckling modes of the plate with a circular hole associated with different boundaries as SSSS, CCCC (thicknesst=1.0).
(a)
Uniaxial SSSS
(b)
Uniaxial CCCC
(c)
Biaxial SSSS
(d)
Biaxial CCCC
(e)
Pure shear SSSS
(f)
Pure shear CCCC


## 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 Figure3. The size of the domain is 0.24×0.12. The mechanical properties of the plate are chosen as E=103, μ=0.25, κ=5/6.Figure 3
Constant stress patch test.For zero shear deformation patch test, the test function [20, 23] can be written as follows:(45)w=a0+a1x+a2y+a3x2+a4xy+a5y2θx=a2+a4x+2a5yθy=-a1-2a3x-a4y.For nonzero shear deformation patch test, the test function can be expressed as follows:(46)w=a0+a1x+a2y+a3x2+a4xy+a5y2+a6x3+a7x2y+a8xy2+a9y3θx=h23k1-μa7+3a9+a2+a4x+2a5y+a7x2+2a8xy+3a9y2θy=-h23k1-μ3a6+a8-a1-2a3x-a4y-3a6x2-2a7xy-a8y2,where ai, (i=0,…,9) are arbitrary constants. In the present paper, we assume that ai=i+1.The mesh of patch test is presented in Figure3. 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 Sx, Sy (Table 2) and the strict patch test with nonzero constant shear stresses Sx, Sy (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.Table 2
Numerical results of constant bending patch test (at node(0.04,0.02)).
w
θ
x
θ
y
M
x
M
y
M
x
y
S
x
S
y
t
=
0.05
Mesh A
1.1528
3.440
-
2.420
-
0.12222
-
0.15556
-
0.041667
0
0
Mesh B
1.1528
3.440
-
2.420
-
0.12222
-
0.15556
-
0.041667
0
0
Mesh C
1.1528
3.440
-
2.420
-
0.12222
-
0.15556
-
0.041667
0
0
Exact
1.1528
3.440
-
2.420
-
0.12222
-
0.15556
-
0.041667
0
0
t
=
0.005
Mesh A
1.1528
3.440
-
2.420
-
0.0001222
-
0.0001556
-
0.000041667
0
0
Mesh B
1.1528
3.440
-
2.420
-
0.0001222
-
0.0001556
-
0.000041667
0
0
Mesh C
1.1528
3.440
-
2.420
-
0.0001222
-
0.0001556
-
0.000041667
0
0
Exact
1.1528
3.440
-
2.420
-
0.0001222
-
0.0001556
-
0.000041667
0
0Table 3
Numerical results of nonzero constant shear stress patch test (at node(0.04,0.02)).
w
θ
x
θ
y
M
x
M
y
M
x
y
S
x
S
y
t
=
0.05
Mesh A
1.1537
3.529867
-
2.51
-
0.156667
-
0.196667
-
0.053
-
0.66667
-
0.84444
Mesh B
1.1537
3.529867
-
2.51
-
0.156667
-
0.196667
-
0.053
-
0.66667
-
0.84444
Mesh C
1.1537
3.529867
-
2.51
-
0.156667
-
0.196667
-
0.053
-
0.66667
-
0.84444
Exact
1.1537
3.529867
-
2.51
-
0.156667
-
0.196667
-
0.053
-
0.66667
-
0.84444
t
=
0.005
Mesh A
1.1537
3.479707
-
2.4704
-
0.0001567
-
0.0001967
-
0.000053
-
0.00067
-
0.000844
Mesh B
1.1537
3.479707
-
2.4704
-
0.0001567
-
0.0001967
-
0.000053
-
0.00067
-
0.000844
Mesh C
1.1537
3.479707
-
2.4704
-
0.0001567
-
0.0001967
-
0.000053
-
0.00067
-
0.000844
Exact
1.1537
3.479707
-
2.4704
-
0.0001567
-
0.0001967
-
0.000053
-
0.00067
-
0.000844


## 4.2. 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 modulusE=2.0×1011, 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/b2D/ρ are obtained. b is the width of the edge, D=Et3/12(1-μ2). Unless otherwise specified, the actual computations are done using dimensionless units.


### 4.2.1. 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 a and the width of the plate b: a/b=1) is divided into mesh of 16×16 elements. The thickness-to-width aspect ratios t/b=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.Table 4
The six lowest nondimensional frequency parameters of a CCCC thin square plate(t/L=0.001).
Element (α)
Mode sequence number
1
2
3
4
5
6
QH8-RCV(-1)
3.6742
7.5726
7.5726
11.2791
13.8025
13.8626
QH8-RCV(-0.5)
3.6602
7.5052
7.5052
11.1242
13.5707
13.6325
QH8-RCV(-0.25)
3.6532
7.4711
7.4711
11.0458
13.4531
13.5157
QH8-RCV(-0.08)
3.6484
7.4478
7.4478
10.9922
13.3726
13.4356
QH8-RCV(0)
3.6461
7.4368
7.4368
10.9669
13.3345
13.3978
QH8-LV
3.6460
7.4363
7.4363
10.9644
13.3318
13.3949
Reference [27]
3.646
7.436
7.436
10.964
13.332
13.395Table 5
The six lowest nondimensional frequency parameters of a CCCC thick square plate(t/L=0.1).
Element (α)
Mode sequence number
1
2
3
4
5
6
QH8-RCV(-1)
3.3118
6.3402
6.3402
8.9003
10.4844
10.5761
QH8-RCV(-0.5)
3.3066
6.3332
6.3332
8.9090
10.5178
10.6160
QH8-RCV(-0.25)
3.3002
6.3038
6.3038
8.8443
10.4237
10.5216
QH8-RCV(-0.08)
3.2943
6.2743
6.2743
8.7752
10.3198
10.4162
QH8-RCV(0)
3.2912
6.2578
6.2578
8.7360
10.2604
10.3558
QH8-LV
3.3181
6.3688
6.3688
8.9622
10.5759
10.6670
Reference [27]
3.292
6.276
6.276
8.792
10.356
10.455Table 6
The six lowest nondimensional frequency parameters of a SSSS thin square plate(t/L=0.001).
Element (α)
Mode sequence number
1
2
3
4
5
6
QH8-RCV(-1)
2.0129
5.0805
5.0805
8.2063
10.3226
10.3226
QH8-RCV(-0.5)
2.0064
5.0403
5.0403
8.1034
10.1616
10.1616
QH8-RCV(-0.25)
2.0032
5.0202
5.0202
8.0522
10.0813
10.0813
QH8-RCV(-0.08)
2.0010
5.0066
5.0066
8.0174
10.0268
10.0268
QH8-RCV(0)
2.0000
5.0002
5.0002
8.0010
10.0012
10.0012
QH8-LV
2.0000
5.0000
5.0000
7.9998
9.9999
9.9999
Reference [27]
2.000
5.000
5.000
8.000
10.000
10.000Table 7
The six lowest nondimensional frequency parameters of a SSSS thick square plate(t/L=0.1).
Element (α)
Mode sequence number
1
2
3
4
5
6
QH8-RCV(-1)
1.9424
4.6599
4.6599
7.1742
8.7513
8.7513
QH8-RCV(-0.5)
1.9376
4.6420
4.6420
7.1505
8.7325
8.7325
QH8-RCV(-0.25)
1.9345
4.6229
4.6229
7.1032
8.6600
8.6600
QH8-RCV(-0.08)
1.9321
4.6062
4.6062
7.0581
8.5881
8.5881
QH8-RCV(0)
1.9308
4.5973
4.5973
7.0335
8.5482
8.5482
QH8-LV
1.9457
4.6784
4.6784
7.2173
8.8157
8.8157
Reference [27]
1.931
4.605
4.605
7.064
8.605
8.605Figure 4
Error of six lowest nondimensional frequency parameters.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 N×N uniform meshes with N=4,8,16. The results agree well with [27] when N=8,16.Table 8
Convergence of six lowest nondimensional frequency parameters (t/h=0.001, 0.1 and SSSS, CCCC).
t
/
L
Boundary
Mesh
Mode sequence number
1
2
3
4
5
6
(a) QH8 with modified consistent mass matrix (QH8-RCV,α=-0.08)
0.001
SSSS
4
×
4
2.0208
5.1600
5.1600
8.5457
10.7782
10.7782
8
×
8
2.0044
5.0290
5.0290
8.0832
10.1251
10.1251
16
×
16
2.0010
5.0066
5.0066
8.0174
10.0268
10.0268
Ref. sol
2.000
5.000
5.000
8.000
10.000
10.000
CCCC
4
×
4
3.6986
7.7478
7.7478
11.8207
14.7803
14.9043
8
×
8
3.6563
7.4895
7.4895
11.1031
13.5317
13.5943
16
×
16
3.6484
7.4478
7.4478
10.9922
13.3726
13.4356
Ref. sol
3.646
7.436
7.436
10.964
13.332
13.395
0.1
SSSS
4
×
4
1.9440
4.6822
4.6822
7.3572
8.9607
8.9607
8
×
8
1.9335
4.6062
4.6062
7.0516
8.5663
8.5663
16
×
16
1.9321
4.6062
4.6062
7.0581
8.5881
8.5881
Ref. sol
1.931
4.605
4.605
7.064
8.605
8.605
CCCC
4
×
4
3.3008
6.3401
6.3401
9.0667
10.7861
10.8677
8
×
8
3.2925
6.2555
6.2555
8.7327
10.2492
10.3400
16
×
16
3.2943
6.2743
6.2743
8.7752
10.3198
10.4162
Ref. sol
3.292
6.276
6.276
8.792
10.356
10.455
(b) QH8 with lumped mass matrix (QH8-LV)
0.001
SSSS
4
×
4
1.9997
4.9970
4.9970
7.9691
9.9694
9.9694
8
×
8
2.0000
4.9999
4.9999
7.9986
10.0001
10.0001
16
×
16
2.0000
5.0000
5.0000
7.9998
9.9999
9.9999
Ref. sol
2.000
5.000
5.000
8.000
10.000
10.000
CCCC
4
×
4
3.6452
7.4121
7.4121
10.8713
12.9885
13.0468
8
×
8
3.6461
7.4365
7.4365
10.9637
13.3312
13.3935
16
×
16
3.6460
7.4363
7.4363
10.9644
13.3318
13.3949
Ref. sol
3.646
7.436
7.436
10.964
13.332
13.395
0.1
SSSS
4
×
4
1.9427
4.6385
4.6385
7.0569
8.5511
8.5511
8
×
8
1.9452
4.6710
4.6710
7.1895
8.7670
8.7670
16
×
16
1.9457
4.6784
4.6784
7.2173
8.8157
8.8157
Ref. sol
1.931
4.605
4.605
7.064
8.605
8.605
CCCC
4
×
4
3.3019
6.2524
6.2524
8.6322
10.0169
10.0917
8
×
8
3.3155
6.3496
6.3496
8.9075
10.4886
10.5773
16
×
16
3.3181
6.3688
6.3688
8.9622
10.5759
10.6670
Ref. sol
3.292
6.276
6.276
8.792
10.356
10.455Figure 5
Convergence of six lowest nondimensional frequency parameters (t/h=0.001, 0.1, and SSSS, CCCC).Then, natural frequencies analyses are further carried on the thin(t/b=0.001) and thick (t/b=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.Table 9
Frequency parameters for square Mindlin plates with different boundary conditions.
Boundary
t
/
b
Element
Mode sequence number
1
2
3
4
5
6
SSFF
0.001
QH8-RCV
0.3407
1.7530
1.9541
3.8675
5.1729
5.4231
QH8-LV
0.3405
1.7485
1.9504
3.8550
5.1471
5.4000
Ref. sol
0.341
1.754
1.955
3.872
5.171
5.419
0.1
QH8-RCV
0.3329
1.6780
1.8751
3.5579
4.7179
4.9445
QH8-LV
0.3335
1.6944
1.8922
3.6122
4.8227
5.0496
Ref. sol
0.333
1.677
1.874
3.557
4.718
4.945
SFSF
0.001
QH8-RCV
0.9759
1.6329
3.7170
3.9486
4.7357
7.1656
QH8-LV
0.9755
1.6303
3.7009
3.9431
4.7240
7.1313
Ref. sol
0.976
1.635
3.721
3.946
4.735
7.167
0.1
QH8-RCV
0.9567
1.5594
3.4297
3.6829
4.3331
6.2859
QH8-LV
0.9601
1.5703
3.4900
3.7305
4.4011
6.4372
Ref. sol
0.956
1.559
3.429
3.681
4.332
6.290
SCSC
0.001
QH8-RCV
2.9350
5.5542
7.0348
9.6060
10.3849
13.1192
QH8-LV
2.9333
5.5466
7.0242
9.5833
10.3566
13.0799
Ref. sol
2.933
5.547
7.024
9.583
10.356
13.080
0.1
QH8-RCV
2.7019
4.9725
5.9898
7.9621
8.7685
10.2141
QH8-LV
2.7212
5.0501
6.0804
8.1367
8.9997
10.4651
Ref. sol
2.700
4.971
5.990
7.972
8.787
10.249
CFFF
0.001
QH8-RCV
0.3516
0.8603
2.1563
2.7540
3.1315
5.4804
QH8-LV
0.3514
0.8590
2.1499
2.7427
3.1211
5.4505
Ref. sol
0.352
0.862
2.157
2.756
3.136
5.490
0.1
QH8-RCV
0.3476
0.8167
2.0353
2.5832
2.8610
4.8108
QH8-LV
0.3481
0.8205
2.0567
2.6220
2.8963
4.9106
Ref. sol
0.348
0.816
2.034
2.582
2.860
4.811
CCFF
0.001
QH8-RCV
0.7003
2.4197
2.6929
4.8226
6.3563
6.6444
QH8-LV
0.6997
2.4117
2.6870
4.8039
6.3216
6.6129
Ref. sol
0.701
2.422
2.694
4.828
6.353
6.640
0.1
QH8-RCV
0.6761
2.2434
2.5046
4.2523
5.5551
5.8093
QH8-LV
0.6777
2.2687
2.5294
4.3213
5.6786
5.9314
Ref. sol
0.676
2.242
2.503
4.251
5.557
5.812
CFCF
0.001
QH8-RCV
2.2459
2.6732
4.4127
6.2028
6.8074
8.0892
QH8-LV
2.2448
2.6689
4.3935
6.1928
6.7890
8.0240
Ref. sol
2.246
2.675
4.417
6.198
6.806
8.087
0.1
QH8-RCV
2.0899
2.4334
3.9032
5.3315
5.7707
6.9214
QH8-LV
2.0979
2.4508
3.9716
5.3982
5.8585
7.1420
Ref. sol
2.088
2.431
3.901
5.331
5.771
6.929Figure 6
Shape of the six lowest modes of SSSS thin square plate obtained from QH8-RCV (α=-0.08, t/h=0.001).


### 4.2.2. Free Vibration of Circular Plate

Free vibration of a circular plate (Figure7) of R=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 E=2.0×1011, and Poisson’s ratio μ=0.3.Figure 7
Circular plate model and its discretization for QH8 element (Mesh A: 205 nodes, 60 elements; Mesh B: 321 nodes, 96 elements).60 proposed elements and 205 nodes as shown in Figure7, 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 (t/2R=0.01) and thick (t/2R=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.Table 10
Ten lowest parameterized natural frequenciesω-=ωR2ρt/D of a clamped circular plate with t/2R=0.01.
Mode
NS-ES-FEM
S8R
MIN3
ES-FEM+DSG3
ANS4
QH8-RCV
QH8-LV
Ref. sol [42]
1
10.2700
10.3095
10.4082
10.3109
10.2572
10.3613
10.3669
10.216
2
21.1535
21.4233
22.2198
21.6702
21.4981
21.6355
21.5152
21.260
3
21.1617
21.4283
22.2444
21.6900
21.4981
21.6356
21.5167
21.260
4
34.4423
35.2809
37.7461
36.3124
35.3941
35.6621
35.1192
34.880
5
34.4573
35.3685
37.7816
36.3816
35.5173
35.7071
35.2812
34.880
6
39.1472
40.1039
43.0344
41.3801
40.8975
40.7676
40.1565
39.771
7
49.6819
52.2399
57.8881
54.7796
52.2054
52.6337
51.3212
51.040
8
49.7157
52.2597
58.0836
54.8922
52.2054
52.6344
51.3358
51.040
9
58.8517
62.2179
68.7260
64.6300
63.2397
63.0358
61.0974
60.820
10
59.0184
62.2543
69.2354
65.1330
63.2397
63.0368
61.1017
60.820Table 11
Ten lowest parameterized natural frequenciesω-=ωR2ρt/D of a clamped circular plate with t/2R=0.1.
Mode
NS-ES-FEM
S8R
MIN3
ES-FEM+DSG3
ANS4
QH8-RCV
QH8-LV
Ref. sol [43]
1
9.2460
9.2539
9.9682
9.3262
9.2605
9.3380
9.4444
9.240
2
17.6424
17.8115
20.306
18.0461
17.9469
17.8602
18.1926
17.834
3
17.6461
17.8115
20.3256
18.0673
17.9469
17.8603
18.1938
17.834
4
26.5699
27.0923
32.7800
27.8438
27.0345
26.9746
27.5602
27.214
5
26.5740
27.1650
32.7905
27.8856
27.6566
26.9836
27.6854
27.214
6
29.7654
30.4101
37.1591
31.1280
30.3221
30.0846
30.8631
30.211
7
35.5595
33.8815
47.5437
38.6936
37.2579
36.5419
37.3990
37.109
8
35.5999
33.8848
47.5993
38.7042
37.2579
36.5431
37.4096
37.109
9
41.0560
37.0290
55.3530
44.1130
43.2702
41.8999
42.8445
42.409
10
41.1226
37.0324
55.6608
44.3286
43.2702
41.9011
42.8465
42.409Figure 8
Ten lowest parameterized natural frequencies of a clamped circular plate.
(a)
Thin circular plate(t/2R=0.01)
(b)
Thick circular plate(t/2R=0.1)First six lowest modes of this clamped circular obtained from QH8-RCV using the mesh of Figure7 Mesh B are shown in Figure 9, which exactly describes the real physical modes of the circular plate.Figure 9
Shape of the six lowest modes of a clamped thick circular plate obtained from QH8-RCV(t/2R=0.1).


### 4.2.3. 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(a). The material parameters of the triangular plate are given as follows: the density mass ρ=8000, Young’s modulus E=2.0×1011, and Poisson’s ratio μ=0.3. The plate is discretized by the meshes as shown in Figure 10 for QH8. The convergence study of nondimensional frequency parameter ω-=ωL2ρt/D/π2 for thin (t/L=0.001) and thick (t/L=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.Table 12
The lowest six parameterized natural frequencies of triangular plates witht/L=0.001.
Skew angle
Method
Mode sequence number
1
2
3
4
5
6
0°
DSG3
0.6288
2.4399
3.4796
5.9913
8.3737
11.2081
MIN3
0.6272
2.4242
3.4344
5.9458
8.2460
11.0178
ES-FEM+DSG3
0.6241
2.3970
3.3631
5.7976
7.9696
10.5356
NS+ES-FEM
0.6253
2.3335
3.1888
5.4772
7.4350
9.2790
Rayleigh-Ritz
0.6240
2.3770
3.3080
5.6890
7.7430
—
ASN4
0.6240
2.3790
3.3170
5.7240
7.7940
10.2000
QH8-RCV
0.6232
2.3709
3.2980
5.6766
7.7627
10.0574
QH8-LV
0.6165
2.3143
3.1843
5.4810
7.5749
9.6271
15°
DSG3
0.5915
2.2360
3.6467
5.5686
7.8875
11.1572
MIN3
0.5886
2.2229
3.5717
5.5286
7.7447
10.9544
ES-FEM+DSG3
0.5843
2.1979
3.4937
5.3893
7.4745
10.4532
NS+ES-FEM
0.5854
2.1449
3.2938
5.0391
7.0313
9.1652
Rayleigh-Ritz
0.5840
2.1810
3.4090
5.2800
7.2640
—
ASN4
0.5830
2.1810
3.4130
5.3030
7.2890
10.0950
QH8-RCV
0.5829
2.1759
3.3995
5.2687
7.2983
9.9691
QH8-LV
0.5750
2.1291
3.2332
5.0577
7.0972
9.4300
30°
DSG3
0.5928
2.2293
3.9729
5.8523
7.8772
11.9930
MIN3
0.5866
2.2164
3.8656
5.8059
7.6697
11.5791
ES-FEM+DSG3
0.5795
2.1881
3.7620
5.6464
7.3608
10.9609
NS+ES-FEM
0.5801
2.1433
3.5402
5.2019
6.8932
9.6373
Rayleigh-Ritz
0.5760
2.1740
3.6390
5.5110
7.1080
—
ASN4
0.5750
2.1740
3.6380
5.5340
7.1390
10.4770
QH8-RCV
0.5763
2.1685
3.6379
5.4957
7.1296
10.3946
QH8-LV
0.5714
2.1368
3.5544
5.3345
6.9689
9.9234
45°
DSG3
0.6265
2.4298
4.5709
6.9888
8.7055
13.4984
MIN3
0.6125
2.4038
4.4237
6.9194
8.3447
12.7848
ES-FEM+DSG3
0.6001
2.3557
4.2652
6.6442
7.9210
11.8923
NS+ES-FEM
0.6020
2.3189
4.0352
6.0299
7.3322
10.5988
Rayleigh-Ritz
0.5900
2.3290
4.1370
6.3810
7.6020
—
ASN4
0.5880
2.3240
4.1260
6.3810
7.6140
11.2240
QH8-RCV
0.5925
2.3290
4.1374
6.4133
7.6019
11.2205
QH8-LV
0.5858
2.3042
3.9784
6.1554
7.4240
10.8429
60°
DSG3
0.6969
2.9282
6.0785
9.3379
11.8340
17.8590
MIN3
0.6631
2.8106
5.8343
8.9727
11.3197
16.3597
ES-FEM+DSG3
0.6371
2.6804
5.4909
8.1353
10.6516
14.5225
NS+ES-FEM
0.6546
2.6792
5.3330
7.3384
9.6825
12.9900
Rayleigh-Ritz
0.6170
2.5760
5.3760
7.5240
10.2850
—
ASN4
0.6130
2.5640
5.3530
7.4600
10.3060
12.9420
QH8-RCV
0.6234
2.5978
5.3502
7.5829
10.2623
13.0926
QH8-LV
0.6164
2.5712
5.1891
7.3141
9.8462
12.8085Table 13
The lowest six parameterized natural frequencies of triangular plates witht/L=0.2.
Skew angle
Method
Mode sequence number
1
2
3
4
5
6
0°
DSG3
0.5872
1.9529
2.4678
4.1305
5.2243
6.2886
MIN3
0.6042
2.1109
2.7737
4.6581
6.1311
7.3208
ES-FEM+DSG3
0.5856
1.9376
2.4411
4.0678
5.1289
6.1359
Rayleigh-Ritz
0.5820
1.9000
2.4080
3.9360
—
—
ASN4
0.5820
1.9150
2.4280
3.9840
5.0180
5.9440
NS+ES-FEM
0.5821
1.8710
2.3600
3.8330
4.8217
5.5683
QH8-RCV
0.5821
1.8907
2.3931
3.8677
4.8701
5.6897
QH8-LV
0.5826
1.9628
2.5320
4.1115
5.2710
6.4029
15°
DSG3
0.5478
1.8160
2.4429
3.8111
5.0253
6.0033
MIN3
0.5656
1.9555
2.7772
4.2851
5.8857
6.9977
ES-FEM+DSG3
0.5457
1.8018
2.4121
3.7529
4.9317
5.8327
Rayleigh-Ritz
0.5440
1.7710
2.3860
3.6280
—
—
ASN4
0.5450
1.7640
2.4200
3.6080
4.8200
5.4310
NS+ES-FEM
0.5442
1.7466
2.3391
3.5289
4.6721
5.1901
QH8-RCV
0.5441
1.7644
2.3687
3.5813
4.7024
5.3504
QH8-LV
0.5427
1.8251
2.5244
3.8353
5.0144
6.1823
30°
DSG3
0.5365
1.8126
2.4954
3.7549
4.9796
5.9219
MIN3
0.5582
1.9551
2.8595
4.2605
5.8273
6.9945
ES-FEM+DSG3
0.5329
1.7963
2.4552
3.6902
4.8717
5.6942
Rayleigh-Ritz
0.5330
1.7720
2.4190
3.5650
—
—
ASN4
0.5320
1.7730
2.4370
3.5910
4.7650
5.3230
NS+ES-FEM
0.5336
1.7485
2.3730
3.4615
4.5927
4.9630
QH8-RCV
0.5331
1.7663
2.4045
3.5167
4.6384
5.1512
QH8-LV
0.5342
1.8391
2.7165
3.9146
4.9778
6.4454
45°
DSG3
0.5448
1.9269
2.6246
3.9196
5.1474
5.9838
MIN3
0.5724
2.1044
3.0165
4.5535
6.0525
7.2736
ES-FEM+DSG3
0.5380
1.9002
2.5682
3.8239
4.9806
5.6000
Rayleigh-Ritz
0.5400
1.8850
2.4890
3.6740
—
—
ASN4
0.5410
1.8840
2.5180
3.7480
4.7400
5.2920
NS+ES-FEM
0.5408
1.8634
2.4323
3.5593
4.4299
4.9328
QH8-RCV
0.5391
1.8765
2.4729
3.6337
4.5699
5.0583
QH8-LV
0.5389
1.9562
2.9904
4.2186
5.2159
6.9240
60°
DSG3
0.5691
2.1281
2.7093
4.3427
5.3437
6.5278
MIN3
0.6017
2.3742
3.2377
5.1762
6.5578
8.1003
ES-FEM+DSG3
0.5569
2.0692
2.5955
4.1109
4.7967
6.0192
Rayleigh-Ritz
0.5590
2.0590
2.3960
3.5900
—
—
ASN4
0.5590
2.0950
2.4830
3.9100
4.5170
5.7630
NS+ES-FEM
0.5611
2.0446
2.3276
3.5969
4.2312
5.3615
QH8-RCV
0.5574
2.0455
2.3807
3.6046
4.2131
5.2640
QH8-LV
0.5573
2.1170
3.8244
4.6656
6.2099
7.7794Figure 10
Triangular plate model with various shape geometries and its discretization meshes (θ=0°, 15°, 30°, 45°, 60°).
(a)
Triangular plate
(b)
0°, 28 elements and 106 nodes
(c)
15°, 26 elements and 103 nodes
(d)
30°, 27 elements and 106 nodes
(e)
45°, 26 elements and 105 nodes
(f)
60°, 28 elements and 107 nodesFigure 11
Convergence study of the lowest four parameterized natural frequencies of triangular plates witht/L=0.001.Figure 12
Convergence study of the lowest four parameterized natural frequencies of triangular plates witht/L=0.2.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.Figure 13
The six lowest modes of a CFF thin triangular plate obtained from QH8-RCV (θ=15°).


## 4.2.1. 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 a and the width of the plate b: a/b=1) is divided into mesh of 16×16 elements. The thickness-to-width aspect ratios t/b=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.Table 4
The six lowest nondimensional frequency parameters of a CCCC thin square plate(t/L=0.001).
Element (α)
Mode sequence number
1
2
3
4
5
6
QH8-RCV(-1)
3.6742
7.5726
7.5726
11.2791
13.8025
13.8626
QH8-RCV(-0.5)
3.6602
7.5052
7.5052
11.1242
13.5707
13.6325
QH8-RCV(-0.25)
3.6532
7.4711
7.4711
11.0458
13.4531
13.5157
QH8-RCV(-0.08)
3.6484
7.4478
7.4478
10.9922
13.3726
13.4356
QH8-RCV(0)
3.6461
7.4368
7.4368
10.9669
13.3345
13.3978
QH8-LV
3.6460
7.4363
7.4363
10.9644
13.3318
13.3949
Reference [27]
3.646
7.436
7.436
10.964
13.332
13.395Table 5
The six lowest nondimensional frequency parameters of a CCCC thick square plate(t/L=0.1).
Element (α)
Mode sequence number
1
2
3
4
5
6
QH8-RCV(-1)
3.3118
6.3402
6.3402
8.9003
10.4844
10.5761
QH8-RCV(-0.5)
3.3066
6.3332
6.3332
8.9090
10.5178
10.6160
QH8-RCV(-0.25)
3.3002
6.3038
6.3038
8.8443
10.4237
10.5216
QH8-RCV(-0.08)
3.2943
6.2743
6.2743
8.7752
10.3198
10.4162
QH8-RCV(0)
3.2912
6.2578
6.2578
8.7360
10.2604
10.3558
QH8-LV
3.3181
6.3688
6.3688
8.9622
10.5759
10.6670
Reference [27]
3.292
6.276
6.276
8.792
10.356
10.455Table 6
The six lowest nondimensional frequency parameters of a SSSS thin square plate(t/L=0.001).
Element (α)
Mode sequence number
1
2
3
4
5
6
QH8-RCV(-1)
2.0129
5.0805
5.0805
8.2063
10.3226
10.3226
QH8-RCV(-0.5)
2.0064
5.0403
5.0403
8.1034
10.1616
10.1616
QH8-RCV(-0.25)
2.0032
5.0202
5.0202
8.0522
10.0813
10.0813
QH8-RCV(-0.08)
2.0010
5.0066
5.0066
8.0174
10.0268
10.0268
QH8-RCV(0)
2.0000
5.0002
5.0002
8.0010
10.0012
10.0012
QH8-LV
2.0000
5.0000
5.0000
7.9998
9.9999
9.9999
Reference [27]
2.000
5.000
5.000
8.000
10.000
10.000Table 7
The six lowest nondimensional frequency parameters of a SSSS thick square plate(t/L=0.1).
Element (α)
Mode sequence number
1
2
3
4
5
6
QH8-RCV(-1)
1.9424
4.6599
4.6599
7.1742
8.7513
8.7513
QH8-RCV(-0.5)
1.9376
4.6420
4.6420
7.1505
8.7325
8.7325
QH8-RCV(-0.25)
1.9345
4.6229
4.6229
7.1032
8.6600
8.6600
QH8-RCV(-0.08)
1.9321
4.6062
4.6062
7.0581
8.5881
8.5881
QH8-RCV(0)
1.9308
4.5973
4.5973
7.0335
8.5482
8.5482
QH8-LV
1.9457
4.6784
4.6784
7.2173
8.8157
8.8157
Reference [27]
1.931
4.605
4.605
7.064
8.605
8.605Figure 4
Error of six lowest nondimensional frequency parameters.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 N×N uniform meshes with N=4,8,16. The results agree well with [27] when N=8,16.Table 8
Convergence of six lowest nondimensional frequency parameters (t/h=0.001, 0.1 and SSSS, CCCC).
t
/
L
Boundary
Mesh
Mode sequence number
1
2
3
4
5
6
(a) QH8 with modified consistent mass matrix (QH8-RCV,α=-0.08)
0.001
SSSS
4
×
4
2.0208
5.1600
5.1600
8.5457
10.7782
10.7782
8
×
8
2.0044
5.0290
5.0290
8.0832
10.1251
10.1251
16
×
16
2.0010
5.0066
5.0066
8.0174
10.0268
10.0268
Ref. sol
2.000
5.000
5.000
8.000
10.000
10.000
CCCC
4
×
4
3.6986
7.7478
7.7478
11.8207
14.7803
14.9043
8
×
8
3.6563
7.4895
7.4895
11.1031
13.5317
13.5943
16
×
16
3.6484
7.4478
7.4478
10.9922
13.3726
13.4356
Ref. sol
3.646
7.436
7.436
10.964
13.332
13.395
0.1
SSSS
4
×
4
1.9440
4.6822
4.6822
7.3572
8.9607
8.9607
8
×
8
1.9335
4.6062
4.6062
7.0516
8.5663
8.5663
16
×
16
1.9321
4.6062
4.6062
7.0581
8.5881
8.5881
Ref. sol
1.931
4.605
4.605
7.064
8.605
8.605
CCCC
4
×
4
3.3008
6.3401
6.3401
9.0667
10.7861
10.8677
8
×
8
3.2925
6.2555
6.2555
8.7327
10.2492
10.3400
16
×
16
3.2943
6.2743
6.2743
8.7752
10.3198
10.4162
Ref. sol
3.292
6.276
6.276
8.792
10.356
10.455
(b) QH8 with lumped mass matrix (QH8-LV)
0.001
SSSS
4
×
4
1.9997
4.9970
4.9970
7.9691
9.9694
9.9694
8
×
8
2.0000
4.9999
4.9999
7.9986
10.0001
10.0001
16
×
16
2.0000
5.0000
5.0000
7.9998
9.9999
9.9999
Ref. sol
2.000
5.000
5.000
8.000
10.000
10.000
CCCC
4
×
4
3.6452
7.4121
7.4121
10.8713
12.9885
13.0468
8
×
8
3.6461
7.4365
7.4365
10.9637
13.3312
13.3935
16
×
16
3.6460
7.4363
7.4363
10.9644
13.3318
13.3949
Ref. sol
3.646
7.436
7.436
10.964
13.332
13.395
0.1
SSSS
4
×
4
1.9427
4.6385
4.6385
7.0569
8.5511
8.5511
8
×
8
1.9452
4.6710
4.6710
7.1895
8.7670
8.7670
16
×
16
1.9457
4.6784
4.6784
7.2173
8.8157
8.8157
Ref. sol
1.931
4.605
4.605
7.064
8.605
8.605
CCCC
4
×
4
3.3019
6.2524
6.2524
8.6322
10.0169
10.0917
8
×
8
3.3155
6.3496
6.3496
8.9075
10.4886
10.5773
16
×
16
3.3181
6.3688
6.3688
8.9622
10.5759
10.6670
Ref. sol
3.292
6.276
6.276
8.792
10.356
10.455Figure 5
Convergence of six lowest nondimensional frequency parameters (t/h=0.001, 0.1, and SSSS, CCCC).Then, natural frequencies analyses are further carried on the thin(t/b=0.001) and thick (t/b=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.Table 9
Frequency parameters for square Mindlin plates with different boundary conditions.
Boundary
t
/
b
Element
Mode sequence number
1
2
3
4
5
6
SSFF
0.001
QH8-RCV
0.3407
1.7530
1.9541
3.8675
5.1729
5.4231
QH8-LV
0.3405
1.7485
1.9504
3.8550
5.1471
5.4000
Ref. sol
0.341
1.754
1.955
3.872
5.171
5.419
0.1
QH8-RCV
0.3329
1.6780
1.8751
3.5579
4.7179
4.9445
QH8-LV
0.3335
1.6944
1.8922
3.6122
4.8227
5.0496
Ref. sol
0.333
1.677
1.874
3.557
4.718
4.945
SFSF
0.001
QH8-RCV
0.9759
1.6329
3.7170
3.9486
4.7357
7.1656
QH8-LV
0.9755
1.6303
3.7009
3.9431
4.7240
7.1313
Ref. sol
0.976
1.635
3.721
3.946
4.735
7.167
0.1
QH8-RCV
0.9567
1.5594
3.4297
3.6829
4.3331
6.2859
QH8-LV
0.9601
1.5703
3.4900
3.7305
4.4011
6.4372
Ref. sol
0.956
1.559
3.429
3.681
4.332
6.290
SCSC
0.001
QH8-RCV
2.9350
5.5542
7.0348
9.6060
10.3849
13.1192
QH8-LV
2.9333
5.5466
7.0242
9.5833
10.3566
13.0799
Ref. sol
2.933
5.547
7.024
9.583
10.356
13.080
0.1
QH8-RCV
2.7019
4.9725
5.9898
7.9621
8.7685
10.2141
QH8-LV
2.7212
5.0501
6.0804
8.1367
8.9997
10.4651
Ref. sol
2.700
4.971
5.990
7.972
8.787
10.249
CFFF
0.001
QH8-RCV
0.3516
0.8603
2.1563
2.7540
3.1315
5.4804
QH8-LV
0.3514
0.8590
2.1499
2.7427
3.1211
5.4505
Ref. sol
0.352
0.862
2.157
2.756
3.136
5.490
0.1
QH8-RCV
0.3476
0.8167
2.0353
2.5832
2.8610
4.8108
QH8-LV
0.3481
0.8205
2.0567
2.6220
2.8963
4.9106
Ref. sol
0.348
0.816
2.034
2.582
2.860
4.811
CCFF
0.001
QH8-RCV
0.7003
2.4197
2.6929
4.8226
6.3563
6.6444
QH8-LV
0.6997
2.4117
2.6870
4.8039
6.3216
6.6129
Ref. sol
0.701
2.422
2.694
4.828
6.353
6.640
0.1
QH8-RCV
0.6761
2.2434
2.5046
4.2523
5.5551
5.8093
QH8-LV
0.6777
2.2687
2.5294
4.3213
5.6786
5.9314
Ref. sol
0.676
2.242
2.503
4.251
5.557
5.812
CFCF
0.001
QH8-RCV
2.2459
2.6732
4.4127
6.2028
6.8074
8.0892
QH8-LV
2.2448
2.6689
4.3935
6.1928
6.7890
8.0240
Ref. sol
2.246
2.675
4.417
6.198
6.806
8.087
0.1
QH8-RCV
2.0899
2.4334
3.9032
5.3315
5.7707
6.9214
QH8-LV
2.0979
2.4508
3.9716
5.3982
5.8585
7.1420
Ref. sol
2.088
2.431
3.901
5.331
5.771
6.929Figure 6
Shape of the six lowest modes of SSSS thin square plate obtained from QH8-RCV (α=-0.08, t/h=0.001).


## 4.2.2. Free Vibration of Circular Plate

Free vibration of a circular plate (Figure7) of R=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 E=2.0×1011, and Poisson’s ratio μ=0.3.Figure 7
Circular plate model and its discretization for QH8 element (Mesh A: 205 nodes, 60 elements; Mesh B: 321 nodes, 96 elements).60 proposed elements and 205 nodes as shown in Figure7, 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 (t/2R=0.01) and thick (t/2R=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.Table 10
Ten lowest parameterized natural frequenciesω-=ωR2ρt/D of a clamped circular plate with t/2R=0.01.
Mode
NS-ES-FEM
S8R
MIN3
ES-FEM+DSG3
ANS4
QH8-RCV
QH8-LV
Ref. sol [42]
1
10.2700
10.3095
10.4082
10.3109
10.2572
10.3613
10.3669
10.216
2
21.1535
21.4233
22.2198
21.6702
21.4981
21.6355
21.5152
21.260
3
21.1617
21.4283
22.2444
21.6900
21.4981
21.6356
21.5167
21.260
4
34.4423
35.2809
37.7461
36.3124
35.3941
35.6621
35.1192
34.880
5
34.4573
35.3685
37.7816
36.3816
35.5173
35.7071
35.2812
34.880
6
39.1472
40.1039
43.0344
41.3801
40.8975
40.7676
40.1565
39.771
7
49.6819
52.2399
57.8881
54.7796
52.2054
52.6337
51.3212
51.040
8
49.7157
52.2597
58.0836
54.8922
52.2054
52.6344
51.3358
51.040
9
58.8517
62.2179
68.7260
64.6300
63.2397
63.0358
61.0974
60.820
10
59.0184
62.2543
69.2354
65.1330
63.2397
63.0368
61.1017
60.820Table 11
Ten lowest parameterized natural frequenciesω-=ωR2ρt/D of a clamped circular plate with t/2R=0.1.
Mode
NS-ES-FEM
S8R
MIN3
ES-FEM+DSG3
ANS4
QH8-RCV
QH8-LV
Ref. sol [43]
1
9.2460
9.2539
9.9682
9.3262
9.2605
9.3380
9.4444
9.240
2
17.6424
17.8115
20.306
18.0461
17.9469
17.8602
18.1926
17.834
3
17.6461
17.8115
20.3256
18.0673
17.9469
17.8603
18.1938
17.834
4
26.5699
27.0923
32.7800
27.8438
27.0345
26.9746
27.5602
27.214
5
26.5740
27.1650
32.7905
27.8856
27.6566
26.9836
27.6854
27.214
6
29.7654
30.4101
37.1591
31.1280
30.3221
30.0846
30.8631
30.211
7
35.5595
33.8815
47.5437
38.6936
37.2579
36.5419
37.3990
37.109
8
35.5999
33.8848
47.5993
38.7042
37.2579
36.5431
37.4096
37.109
9
41.0560
37.0290
55.3530
44.1130
43.2702
41.8999
42.8445
42.409
10
41.1226
37.0324
55.6608
44.3286
43.2702
41.9011
42.8465
42.409Figure 8
Ten lowest parameterized natural frequencies of a clamped circular plate.
(a)
Thin circular plate(t/2R=0.01)
(b)
Thick circular plate(t/2R=0.1)First six lowest modes of this clamped circular obtained from QH8-RCV using the mesh of Figure7 Mesh B are shown in Figure 9, which exactly describes the real physical modes of the circular plate.Figure 9
Shape of the six lowest modes of a clamped thick circular plate obtained from QH8-RCV(t/2R=0.1).


## 4.2.3. 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(a). The material parameters of the triangular plate are given as follows: the density mass ρ=8000, Young’s modulus E=2.0×1011, and Poisson’s ratio μ=0.3. The plate is discretized by the meshes as shown in Figure 10 for QH8. The convergence study of nondimensional frequency parameter ω-=ωL2ρt/D/π2 for thin (t/L=0.001) and thick (t/L=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.Table 12
The lowest six parameterized natural frequencies of triangular plates witht/L=0.001.
Skew angle
Method
Mode sequence number
1
2
3
4
5
6
0°
DSG3
0.6288
2.4399
3.4796
5.9913
8.3737
11.2081
MIN3
0.6272
2.4242
3.4344
5.9458
8.2460
11.0178
ES-FEM+DSG3
0.6241
2.3970
3.3631
5.7976
7.9696
10.5356
NS+ES-FEM
0.6253
2.3335
3.1888
5.4772
7.4350
9.2790
Rayleigh-Ritz
0.6240
2.3770
3.3080
5.6890
7.7430
—
ASN4
0.6240
2.3790
3.3170
5.7240
7.7940
10.2000
QH8-RCV
0.6232
2.3709
3.2980
5.6766
7.7627
10.0574
QH8-LV
0.6165
2.3143
3.1843
5.4810
7.5749
9.6271
15°
DSG3
0.5915
2.2360
3.6467
5.5686
7.8875
11.1572
MIN3
0.5886
2.2229
3.5717
5.5286
7.7447
10.9544
ES-FEM+DSG3
0.5843
2.1979
3.4937
5.3893
7.4745
10.4532
NS+ES-FEM
0.5854
2.1449
3.2938
5.0391
7.0313
9.1652
Rayleigh-Ritz
0.5840
2.1810
3.4090
5.2800
7.2640
—
ASN4
0.5830
2.1810
3.4130
5.3030
7.2890
10.0950
QH8-RCV
0.5829
2.1759
3.3995
5.2687
7.2983
9.9691
QH8-LV
0.5750
2.1291
3.2332
5.0577
7.0972
9.4300
30°
DSG3
0.5928
2.2293
3.9729
5.8523
7.8772
11.9930
MIN3
0.5866
2.2164
3.8656
5.8059
7.6697
11.5791
ES-FEM+DSG3
0.5795
2.1881
3.7620
5.6464
7.3608
10.9609
NS+ES-FEM
0.5801
2.1433
3.5402
5.2019
6.8932
9.6373
Rayleigh-Ritz
0.5760
2.1740
3.6390
5.5110
7.1080
—
ASN4
0.5750
2.1740
3.6380
5.5340
7.1390
10.4770
QH8-RCV
0.5763
2.1685
3.6379
5.4957
7.1296
10.3946
QH8-LV
0.5714
2.1368
3.5544
5.3345
6.9689
9.9234
45°
DSG3
0.6265
2.4298
4.5709
6.9888
8.7055
13.4984
MIN3
0.6125
2.4038
4.4237
6.9194
8.3447
12.7848
ES-FEM+DSG3
0.6001
2.3557
4.2652
6.6442
7.9210
11.8923
NS+ES-FEM
0.6020
2.3189
4.0352
6.0299
7.3322
10.5988
Rayleigh-Ritz
0.5900
2.3290
4.1370
6.3810
7.6020
—
ASN4
0.5880
2.3240
4.1260
6.3810
7.6140
11.2240
QH8-RCV
0.5925
2.3290
4.1374
6.4133
7.6019
11.2205
QH8-LV
0.5858
2.3042
3.9784
6.1554
7.4240
10.8429
60°
DSG3
0.6969
2.9282
6.0785
9.3379
11.8340
17.8590
MIN3
0.6631
2.8106
5.8343
8.9727
11.3197
16.3597
ES-FEM+DSG3
0.6371
2.6804
5.4909
8.1353
10.6516
14.5225
NS+ES-FEM
0.6546
2.6792
5.3330
7.3384
9.6825
12.9900
Rayleigh-Ritz
0.6170
2.5760
5.3760
7.5240
10.2850
—
ASN4
0.6130
2.5640
5.3530
7.4600
10.3060
12.9420
QH8-RCV
0.6234
2.5978
5.3502
7.5829
10.2623
13.0926
QH8-LV
0.6164
2.5712
5.1891
7.3141
9.8462
12.8085Table 13
The lowest six parameterized natural frequencies of triangular plates witht/L=0.2.
Skew angle
Method
Mode sequence number
1
2
3
4
5
6
0°
DSG3
0.5872
1.9529
2.4678
4.1305
5.2243
6.2886
MIN3
0.6042
2.1109
2.7737
4.6581
6.1311
7.3208
ES-FEM+DSG3
0.5856
1.9376
2.4411
4.0678
5.1289
6.1359
Rayleigh-Ritz
0.5820
1.9000
2.4080
3.9360
—
—
ASN4
0.5820
1.9150
2.4280
3.9840
5.0180
5.9440
NS+ES-FEM
0.5821
1.8710
2.3600
3.8330
4.8217
5.5683
QH8-RCV
0.5821
1.8907
2.3931
3.8677
4.8701
5.6897
QH8-LV
0.5826
1.9628
2.5320
4.1115
5.2710
6.4029
15°
DSG3
0.5478
1.8160
2.4429
3.8111
5.0253
6.0033
MIN3
0.5656
1.9555
2.7772
4.2851
5.8857
6.9977
ES-FEM+DSG3
0.5457
1.8018
2.4121
3.7529
4.9317
5.8327
Rayleigh-Ritz
0.5440
1.7710
2.3860
3.6280
—
—
ASN4
0.5450
1.7640
2.4200
3.6080
4.8200
5.4310
NS+ES-FEM
0.5442
1.7466
2.3391
3.5289
4.6721
5.1901
QH8-RCV
0.5441
1.7644
2.3687
3.5813
4.7024
5.3504
QH8-LV
0.5427
1.8251
2.5244
3.8353
5.0144
6.1823
30°
DSG3
0.5365
1.8126
2.4954
3.7549
4.9796
5.9219
MIN3
0.5582
1.9551
2.8595
4.2605
5.8273
6.9945
ES-FEM+DSG3
0.5329
1.7963
2.4552
3.6902
4.8717
5.6942
Rayleigh-Ritz
0.5330
1.7720
2.4190
3.5650
—
—
ASN4
0.5320
1.7730
2.4370
3.5910
4.7650
5.3230
NS+ES-FEM
0.5336
1.7485
2.3730
3.4615
4.5927
4.9630
QH8-RCV
0.5331
1.7663
2.4045
3.5167
4.6384
5.1512
QH8-LV
0.5342
1.8391
2.7165
3.9146
4.9778
6.4454
45°
DSG3
0.5448
1.9269
2.6246
3.9196
5.1474
5.9838
MIN3
0.5724
2.1044
3.0165
4.5535
6.0525
7.2736
ES-FEM+DSG3
0.5380
1.9002
2.5682
3.8239
4.9806
5.6000
Rayleigh-Ritz
0.5400
1.8850
2.4890
3.6740
—
—
ASN4
0.5410
1.8840
2.5180
3.7480
4.7400
5.2920
NS+ES-FEM
0.5408
1.8634
2.4323
3.5593
4.4299
4.9328
QH8-RCV
0.5391
1.8765
2.4729
3.6337
4.5699
5.0583
QH8-LV
0.5389
1.9562
2.9904
4.2186
5.2159
6.9240
60°
DSG3
0.5691
2.1281
2.7093
4.3427
5.3437
6.5278
MIN3
0.6017
2.3742
3.2377
5.1762
6.5578
8.1003
ES-FEM+DSG3
0.5569
2.0692
2.5955
4.1109
4.7967
6.0192
Rayleigh-Ritz
0.5590
2.0590
2.3960
3.5900
—
—
ASN4
0.5590
2.0950
2.4830
3.9100
4.5170
5.7630
NS+ES-FEM
0.5611
2.0446
2.3276
3.5969
4.2312
5.3615
QH8-RCV
0.5574
2.0455
2.3807
3.6046
4.2131
5.2640
QH8-LV
0.5573
2.1170
3.8244
4.6656
6.2099
7.7794Figure 10
Triangular plate model with various shape geometries and its discretization meshes (θ=0°, 15°, 30°, 45°, 60°).
(a)
Triangular plate
(b)
0°, 28 elements and 106 nodes
(c)
15°, 26 elements and 103 nodes
(d)
30°, 27 elements and 106 nodes
(e)
45°, 26 elements and 105 nodes
(f)
60°, 28 elements and 107 nodesFigure 11
Convergence study of the lowest four parameterized natural frequencies of triangular plates witht/L=0.001.Figure 12
Convergence study of the lowest four parameterized natural frequencies of triangular plates witht/L=0.2.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.Figure 13
The six lowest modes of a CFF thin triangular plate obtained from QH8-RCV (θ=15°).


## 4.3. 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.(47)Kb=λb2π2D,with b being the edge width of the plate that is used throughout the paper and D=Et3/12(1-μ2).


### 4.3.1. Rectangular Plates Subjected to an In-Plane Compressive Load

The parameterα0 in (42) for the present element QH8-R (8-node hybrid element with refined geometric stiffness matrix) is determined by comparison using numerical examples. A square plate (a/b=1) compressed in one direction or two perpendicular directions are divided into different meshes 2×2, 4×4, 6×6, 8×8, 12×12, 16×16. Thickness-to-width aspect ratios t/b=0.01, 0.1, α0=-1,-0.75,-0.5,-0.3,-0.17,0, and boundary conditions SSSS, CCCC are chosen. The computational results of the critical load are given in Tables 14 and 15 and Figure 14, in which the exact solutions can be found in [44–46]. After analysis of numerical examples, the present 8-node hybrid stress element QH8-R(α0) has higher accuracy when -0.4≤α0≤-0.1, and we advise α0=-0.17 for following numerical examples.Table 14
The factor of uniaxial buckling loads along thex-axis of rectangular plates with length-to-width aspect ratios a/b=1 and thickness-to-width aspect ratios t/b=0.01.
Boundary condition
Element (α0)
Mesh
2
×
2
4
×
4
6
×
6
8
×
8
12
×
12
16
×
16
SSSS
QH8-R(-1)
3.1427
3.4982
3.7191
3.8183
3.8982
3.9293
(
-
0.75
)
4.0364
3.8502
3.9045
3.9338
3.9588
3.9691
(
-
0.5
)
4.7314
4.0776
4.0190
4.0043
3.9961
3.9940
(
-
0.3
)
4.8659
4.1389
4.0521
4.0258
4.0084
4.0027
(
-
0.17
)
4.7097
4.1161
4.0444
4.0225
4.0079
4.0031
(
0
)
4.2694
4.0139
4.0002
3.9980
3.9972
3.9971
Reference [44]
4
CCCC
QH8-R(-1)
5.3071
7.4051
8.3520
8.9070
9.4027
9.6049
(
-
0.75
)
8.2670
9.1692
9.4445
9.6292
9.7956
9.8651
(
-
0.5
)
10.9935
10.6774
10.2199
10.1074
10.0459
10.0304
(
-
0.3
)
12.9592
11.1920
10.4649
10.2587
10.1293
10.0885
(
-
0.17
)
13.6765
11.0371
10.4113
10.2349
10.1247
10.0896
(
0
)
13.6316
10.2573
10.0875
10.0576
10.0482
10.0472
Reference [44]
10.0738
Note. κ=5/6.Table 15
The factor of uniaxial buckling loads along thex-axis of rectangular plates with length-to-width aspect ratios a/b=1 and thickness-to-width aspect ratios t/b=0.1.
Boundary
Element (α0)
Mesh
2
×
2
4
×
4
6
×
6
8
×
8
12
×
12
16
×
16
SSSS
QH8-R(-1)
2.3811
2.4760
2.7861
3.0392
3.3424
3.4914
(
-
0.75
)
3.3240
3.0601
3.2413
3.3983
3.5815
3.6689
(
-
0.5
)
4.1951
3.5827
3.6110
3.6676
3.7403
3.7756
(
-
0.3
)
4.4995
3.8245
3.7755
3.7778
3.7917
3.8002
(
-
0.17
)
4.3933
3.8501
3.7970
3.7869
3.7845
3.7850
(
0
)
3.9119
3.6966
3.7056
3.7146
3.7232
3.7266
Reference [45]
3.7864
CCCC
QH8-R(-1)
3.4293
4.0301
4.6232
5.2042
6.1652
6.7926
(
-
0.75
)
5.9999
5.8082
6.0193
6.4602
7.1856
7.6362
(
-
0.5
)
8.7108
7.7026
7.4918
7.6519
7.9938
8.2101
(
-
0.3
)
10.4843
8.8128
8.3089
8.2536
8.3146
8.3741
(
-
0.17
)
10.9492
8.8625
8.3903
8.2963
8.2811
8.2936
(
0
)
9.9368
7.5669
7.5101
7.5613
7.6923
7.7874
Reference [46]
8.2917
Note. κ=5/6.Figure 14
The factor of uniaxial buckling loads (σx0) along the x-axis of rectangular plates for QH8-R (α0) 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 are investigated. Several thickness-span aspect ratiost/b=0.1, 0.2 (moderately thick plate) and t/b=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 nodes) and MK, Shear-MK use a regular pattern of 21×21 elements (441 nodes). A shear correction factor κ=π2/12 is used here for a proper comparison.Table 16
The factor of uniaxial buckling loads (σx0) along the x-axis of rectangular plates (a/b=1) with various boundary conditions and various thickness-to-width ratios (t/b).
Boundary
Method
t
/
b
0.001
0.05
0.1
0.2
SSSS
SSM
4.000
3.944
3.784
3.256
Shear-MK
4.0243
3.9153
3.7895
3.2748
MK
5.467
5.0111
3.7905
3.2799
Exact [45]
4.000
3.9437
3.7838
3.2558
S8R
3.997
3.9445
3.7866
1.7025
QH8-R
4.0132
3.9478
3.7827
3.2464
SFSC
SSM
1.653
1.620
1.556
1.370
Shear-MK
1.6552
1.6211
1.5658
1.3765
MK
2.5876
2.1778
1.5699
1.3612
Exact [45]
1.6522
1.6197
1.5558
1.3701
S8R
1.6560
1.6202
1.5570
1.2532
QH8-R
1.6522
1.6209
1.5549
1.3669
SFSS
SSM
1.402
1.378
1.327
1.173
Shear-MK
1.4211
1.3877
1.3795
1.2700
MK
2.2456
2.0323
1.3791
1.2755
Exact [45]
1.4014
1.3813
1.3707
1.2138
S8R
1.4017
1.3816
1.3415
1.2158
QH8-R
1.4022
1.3827
1.3410
1.2133
SFSF
SSM
0.9523
0.9412
0.9146
0.8274
Shear-MK
0.9676
0.9455
0.9231
0.8575
MK
1.7644
1.6578
0.9344
0.8654
Exact [45]
0.95225
0.94314
0.92187
0.85011
S8R
0.9524
0.9433
0.9223
0.8325
QH8-R
0.9532
0.9442
0.9227
0.8519
Note. QH8-R (α0=-0.17). Shear correction factor κ=π2/12.Figure 15
Normalized uniaxial buckling load factors along thex-axis (σx0) versus the thickness-span aspect ratio t/b=0.001,0.05,0.1,0.2 with various boundaries.Through the achieved results shown in Table16 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, t/b=0.001, 0.05, but conversely matches well with thick plates, that is, t/b=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.The results of uniaxial buckling loads(σy0) along they-axis and the biaxial in-plane compression load (σx0=σy0) are shown in Tables 17 and 18 and Figures 16 and 17 associated with five different densities of regular elements such as 2×2 (21 nodes), 4×4 (65 nodes), 6×6 (133 nodes), 8×8 (225 nodes), and 10×10 (341 nodes).Table 17
The factor of uniaxial buckling loads (σy0) along the y-axis of rectangular plates (a/b=1) with various boundary conditions and various thickness-to-width ratios (t/b).
Boundary
Element
t
/
b
0.001
0.05
0.1
0.2
SSSS
2
×
2
4.7132
4.6261
4.3905
3.6771
4
×
4
4.1193
4.0431
3.8479
3.2455
6
×
6
4.0476
3.9751
3.7949
3.2231
8
×
8
4.0257
3.9564
3.785
3.2321
10
×
10
4.0161
3.9495
3.7829
3.2422
Exact [45]
4.0000
3.9437
3.7838
3.2558
SFSC
2
×
2
2.5834
2.5216
2.3678
1.9583
4
×
4
2.3993
2.3194
2.1530
1.7580
6
×
6
2.3821
2.2898
2.1191
1.7288
8
×
8
2.3797
2.2794
2.1079
1.7213
10
×
10
2.3799
2.2743
2.1031
1.7190
Exact [45]
2.3901
2.2667
2.1010
1.7200
SFSS
2
×
2
2.5626
2.502
2.3511
1.9480
4
×
4
2.3741
2.2965
2.1346
1.7483
6
×
6
2.3567
2.2672
2.1010
1.7195
8
×
8
2.3541
2.2569
2.0901
1.7122
10
×
10
2.3542
2.2518
2.0853
1.7099
Exact [45]
2.3639
2.2442
2.0829
1.7105
SFSF
2
×
2
2.1000
2.058
1.9553
1.6775
4
×
4
2.0327
1.9737
1.8516
1.5567
6
×
6
2.0286
1.9589
1.8315
1.5358
8
×
8
2.0295
1.9531
1.8244
1.5301
10
×
10
2.0309
1.9500
1.8212
1.5282
Exact [45]
2.0413
1.9457
1.8216
1.5333
Note. QH8-R (α0=-0.17). Shear correction factor κ=π2/12.Table 18
The factor of biaxial buckling loads (σx0=σy0) for rectangular plates (a/b=1) with various boundary conditions and various thickness-to-width ratios (t/b).
Boundary
Element
t
/
b
0.001
0.005
0.05
0.1
SSSS
2
×
2
2.3774
2.3770
2.3377
2.2277
4
×
4
2.0671
2.0667
2.0334
1.9432
6
×
6
2.0275
2.0272
1.9954
1.9108
8
×
8
2.015
2.0147
1.984
1.9023
10
×
10
2.0095
2.0091
1.9793
1.8992
Exact [45]a
2.0000
1.9997
1.9718
1.8919
SFSC
2
×
2
1.1820
1.1818
1.1652
1.1205
4
×
4
1.1474
1.1472
1.1261
1.0767
6
×
6
1.1427
1.1424
1.1184
1.0679
8
×
8
1.1416
1.1413
1.1155
1.0649
10
×
10
1.1414
1.1410
1.1140
1.0636
Exact [45]b
1.1431
1.1412
1.1119
1.0641
SFSS
2
×
2
1.0925
1.0923
1.0794
1.0439
4
×
4
1.0606
1.0605
1.0445
1.0060
6
×
6
1.0558
1.0556
1.0378
0.9986
8
×
8
1.0546
1.0543
1.0352
0.9962
10
×
10
1.0541
1.0538
1.0340
0.9952
Exact [45]b
1.0548
1.0535
1.0322
0.99541
SFSF
2
×
2
0.9671
0.9670
0.9583
0.9334
4
×
4
0.9400
0.9399
0.9302
0.9053
6
×
6
0.9351
0.9349
0.9247
0.9002
8
×
8
0.9335
0.9333
0.9228
0.8986
10
×
10
0.9328
0.9326
0.9219
0.8979
Exact [45]b
0.93209
0.9316
0.92071
0.89774
Note. QH8-R (α0=-0.17). aShear correction factor κ=0.823045. bShear correction factor κ=5/6.Figure 16
Normalized axial buckling load factors along they-axis (σy0) versus the thickness-span aspect ratio t/b=0.001, 0.05, 0.1, 0.2 with various boundaries.Figure 17
Normalized biaxial buckling load factors (σx0=σy0) versus the thickness-span aspect ratio t/b=0.001, 0.05, 0.1, 0.2 with various boundaries.We also consider simply supported plates (SSSS) with various thickness-to-width ratios,t/b = 0.05, 0.1, 0.2, and length-to-width ratios, a/b = 0.5, 1.0, 1.5, 2.0, 2.5.Table19 and Figure 18 show the buckling factors using the element QH8-R with regular mesh of 9×9 rectangular elements (280 nodes) when a/b=1.0. The QH8-R results are also compared to the DSG3 and ES-DSG3 [47] (mesh of 16×16 elements and 289 nodes), S8R (mesh of 9×9 elements and 280 nodes), mesh-free method [48], and the pb-2 Ritz [46]. It is seen that the QH8-R exhibits a good agreement with the pb-2 Ritz method, and S8R is not suitable for plate t/b = 0.2. Figure 19 also depicts the uniaxial buckling modes of simply supported rectangular plates with thickness-to- width ratios t/b = 0.01 and various length-to-width ratios, a/b = 1.0, 1.5, 2.0, 2.5.Table 19
The factor of uniaxial buckling loads along thex-axis of rectangular plates (SSSS) with various length-to-width ratios and various thickness-to-width ratios.
a
/
b
t
/
b
S8R
ES-DSG3 [47]
Mesh-free [48]
QH8-R
pb-2 Ritz [46]
0.5
0.05
6.0376
5.9873
6.0405
6.0822
6.0372
0.1
5.4779
5.3064
5.3116
5.5071
5.4777
0.2
1.9525
3.720
3.7157
3.9754
3.9963
1
0.05
3.9441
3.9412
3.9293
3.9528
3.9444
0.1
3.7861
3.7402
3.7270
3.7855
3.7865
0.2
2.0444
3.1263
3.1471
3.2427
3.2637
1.5
0.05
4.2576
4.2852
4.2116
4.2745
4.2570
0.1
4.0255
3.9844
3.8982
4.0332
4.0250
0.2
2.0132
3.1461
3.1032
3.2966
3.3048
2
0.05
3.9440
3.9811
3.8657
3.9528
3.9444
0.1
3.7860
3.7711
3.6797
3.7855
3.7865
0.2
1.9486
3.1415
3.0783
3.2427
3.2637
2.5
0.05
4.0673
4.1691
3.9600
4.0795
4.0645
0.1
3.8700
3.8924
3.7311
3.8737
3.8683
0.2
1.8697
3.1234
3.0306
3.2338
3.2421Figure 18
The factor of axial buckling loads of SSSS plates with various length-to-width ratios (a/b) and various thickness-to-width ratios (t/b).Figure 19
Axial buckling modes of simply supported rectangular plates with thickness-to-width ratiost/b = 0.01 and various length-to-width ratios a/b.
(a)
a
/
b
=
1.0
(b)
a
/
b
=
1.5
(c)
a
/
b
=
2.0
(d)
a
/
b
=
2.5


### 4.3.2. 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, t/b, 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.Figure 20
Skew plate under uniaxial in-plane loading.Numerical results for thick plates with thickness-to-width ratios,t/b = 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.Table 20
Buckling load intensity factor of skew plates havingS∗S∗S∗S∗, CCCC, CFCF, and S∗FS∗F boundary conditions subject to uniaxial loads (σx0).
Boundary
t
/
b
Skew angle
Mesh-free [48]
6
×
6
9
×
9
16
×
16
pb-2 Ritz [46]
S
∗
S
∗
S
∗
S
∗
0.05
0°
—
3.8057
3.7905
3.7815
3.7835
15°
—
4.1583
4.1380
4.1276
4.1314
30°
5.3874
5.4755
5.4304
5.4109
5.4182
45°
8.6926
8.9757
8.8002
8.7247
8.7328
0.1
0°
—
3.5111
3.5007
3.4965
3.4950
15°
—
3.8234
3.8094
3.8047
3.8027
30°
4.8534
4.9813
4.9475
4.9370
4.9324
45°
7.6106
7.9388
7.7931
7.7430
7.7236
CCCC
0.05
0°
9.5418
9.8091
9.6238
9.5647
9.5588
15°
10.2123
10.5549
10.3209
10.2430
10.2312
30°
12.5225
13.2413
12.7803
12.6114
12.5711
45°
17.8211
19.9900
18.6615
18.1211
17.9652
0.1
0°
8.1849
8.3903
8.2836
8.2936
8.2917
15°
8.7560
8.9338
8.7863
8.7845
8.7741
30°
10.3273
10.8235
10.4883
10.4238
10.3760
45°
13.2881
15.2116
14.2092
13.8638
13.6909
CFCF
0.05
0°
3.7984
3.8724
3.8281
3.8082
3.8007
15°
4.1310
4.2202
4.1688
4.1460
4.1387
30°
5.3510
5.4976
5.4102
5.3720
5.3660
45°
7.3658
7.9413
7.6341
7.4986
7.4670
0.1
0°
3.5063
3.5523
3.5262
3.5185
3.5077
15°
3.7684
3.8441
3.8135
3.8049
3.7937
30°
4.7773
4.8862
4.8316
4.8168
4.8043
45°
6.3092
6.6886
6.4547
6.3680
6.3311
S
∗
F
S
∗
F
0.05
0°
—
0.9459
0.9436
0.9424
0.9421
15°
—
1.0535
1.0513
1.0501
1.0501
30°
—
1.4612
1.4594
1.4583
1.4603
45°
—
2.5333
2.5317
2.5311
2.5450
0.1
0°
—
0.9227
0.9210
0.9202
0.9199
15°
—
1.0215
1.0199
1.0191
1.0189
30°
—
1.3892
1.3878
1.3873
1.3887
45°
—
2.3192
2.3187
2.3190
2.3316
Note. 6×6, 9×9, 16×16: the meshes for element QH8-R(-0.17).


### 4.3.3. 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 lengtha=b=10 and the thickness t=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 x-direction, σx0, pure shear load, and biaxial compressive loads in x and y directions, where σx0=σy0, respectively.Figure 21
Square plate with a circular hole and the distributions (308 nodes, 84 elements).The results of the normalized buckling load factor calculated for those loads with different boundary conditions are shown in Table21. The results agree very well with those computed by the EFG [57] involving both the third-order 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.Table 21
The uniaxial, biaxial, and pure shear buckling load factor for a square plate with a circular hole and different boundary conditions.
Buckling load
Boundary
EFG (FSDT)
EFG (TSDT)
Shear-MK
HQ8-R
Uniaxial
CCCC
7.995
8.097
8.0777
8.0543
SSSS
1.986
1.969
1.9544
2.0580
Biaxial
CCCC
4.781
4.860
4.9878
4.8477
SSSS
1.032
1.021
1.0511
1.0696
Pure shear
CCCC
12.669
13.357
13.5252
12.6976
SSSS
7.867
7.873
8.1544
7.8971Figure 22
The first uniaxial, biaxial, and pear shear buckling modes of the plate with a circular hole associated with different boundaries as SSSS, CCCC (thicknesst=1.0).
(a)
Uniaxial SSSS
(b)
Uniaxial CCCC
(c)
Biaxial SSSS
(d)
Biaxial CCCC
(e)
Pure shear SSSS
(f)
Pure shear CCCC


## 4.3.1. Rectangular Plates Subjected to an In-Plane Compressive Load

The parameterα0 in (42) for the present element QH8-R (8-node hybrid element with refined geometric stiffness matrix) is determined by comparison using numerical examples. A square plate (a/b=1) compressed in one direction or two perpendicular directions are divided into different meshes 2×2, 4×4, 6×6, 8×8, 12×12, 16×16. Thickness-to-width aspect ratios t/b=0.01, 0.1, α0=-1,-0.75,-0.5,-0.3,-0.17,0, and boundary conditions SSSS, CCCC are chosen. The computational results of the critical load are given in Tables 14 and 15 and Figure 14, in which the exact solutions can be found in [44–46]. After analysis of numerical examples, the present 8-node hybrid stress element QH8-R(α0) has higher accuracy when -0.4≤α0≤-0.1, and we advise α0=-0.17 for following numerical examples.Table 14
The factor of uniaxial buckling loads along thex-axis of rectangular plates with length-to-width aspect ratios a/b=1 and thickness-to-width aspect ratios t/b=0.01.
Boundary condition
Element (α0)
Mesh
2
×
2
4
×
4
6
×
6
8
×
8
12
×
12
16
×
16
SSSS
QH8-R(-1)
3.1427
3.4982
3.7191
3.8183
3.8982
3.9293
(
-
0.75
)
4.0364
3.8502
3.9045
3.9338
3.9588
3.9691
(
-
0.5
)
4.7314
4.0776
4.0190
4.0043
3.9961
3.9940
(
-
0.3
)
4.8659
4.1389
4.0521
4.0258
4.0084
4.0027
(
-
0.17
)
4.7097
4.1161
4.0444
4.0225
4.0079
4.0031
(
0
)
4.2694
4.0139
4.0002
3.9980
3.9972
3.9971
Reference [44]
4
CCCC
QH8-R(-1)
5.3071
7.4051
8.3520
8.9070
9.4027
9.6049
(
-
0.75
)
8.2670
9.1692
9.4445
9.6292
9.7956
9.8651
(
-
0.5
)
10.9935
10.6774
10.2199
10.1074
10.0459
10.0304
(
-
0.3
)
12.9592
11.1920
10.4649
10.2587
10.1293
10.0885
(
-
0.17
)
13.6765
11.0371
10.4113
10.2349
10.1247
10.0896
(
0
)
13.6316
10.2573
10.0875
10.0576
10.0482
10.0472
Reference [44]
10.0738
Note. κ=5/6.Table 15
The factor of uniaxial buckling loads along thex-axis of rectangular plates with length-to-width aspect ratios a/b=1 and thickness-to-width aspect ratios t/b=0.1.
Boundary
Element (α0)
Mesh
2
×
2
4
×
4
6
×
6
8
×
8
12
×
12
16
×
16
SSSS
QH8-R(-1)
2.3811
2.4760
2.7861
3.0392
3.3424
3.4914
(
-
0.75
)
3.3240
3.0601
3.2413
3.3983
3.5815
3.6689
(
-
0.5
)
4.1951
3.5827
3.6110
3.6676
3.7403
3.7756
(
-
0.3
)
4.4995
3.8245
3.7755
3.7778
3.7917
3.8002
(
-
0.17
)
4.3933
3.8501
3.7970
3.7869
3.7845
3.7850
(
0
)
3.9119
3.6966
3.7056
3.7146
3.7232
3.7266
Reference [45]
3.7864
CCCC
QH8-R(-1)
3.4293
4.0301
4.6232
5.2042
6.1652
6.7926
(
-
0.75
)
5.9999
5.8082
6.0193
6.4602
7.1856
7.6362
(
-
0.5
)
8.7108
7.7026
7.4918
7.6519
7.9938
8.2101
(
-
0.3
)
10.4843
8.8128
8.3089
8.2536
8.3146
8.3741
(
-
0.17
)
10.9492
8.8625
8.3903
8.2963
8.2811
8.2936
(
0
)
9.9368
7.5669
7.5101
7.5613
7.6923
7.7874
Reference [46]
8.2917
Note. κ=5/6.Figure 14
The factor of uniaxial buckling loads (σx0) along the x-axis of rectangular plates for QH8-R (α0) 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 are investigated. Several thickness-span aspect ratiost/b=0.1, 0.2 (moderately thick plate) and t/b=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 nodes) and MK, Shear-MK use a regular pattern of 21×21 elements (441 nodes). A shear correction factor κ=π2/12 is used here for a proper comparison.Table 16
The factor of uniaxial buckling loads (σx0) along the x-axis of rectangular plates (a/b=1) with various boundary conditions and various thickness-to-width ratios (t/b).
Boundary
Method
t
/
b
0.001
0.05
0.1
0.2
SSSS
SSM
4.000
3.944
3.784
3.256
Shear-MK
4.0243
3.9153
3.7895
3.2748
MK
5.467
5.0111
3.7905
3.2799
Exact [45]
4.000
3.9437
3.7838
3.2558
S8R
3.997
3.9445
3.7866
1.7025
QH8-R
4.0132
3.9478
3.7827
3.2464
SFSC
SSM
1.653
1.620
1.556
1.370
Shear-MK
1.6552
1.6211
1.5658
1.3765
MK
2.5876
2.1778
1.5699
1.3612
Exact [45]
1.6522
1.6197
1.5558
1.3701
S8R
1.6560
1.6202
1.5570
1.2532
QH8-R
1.6522
1.6209
1.5549
1.3669
SFSS
SSM
1.402
1.378
1.327
1.173
Shear-MK
1.4211
1.3877
1.3795
1.2700
MK
2.2456
2.0323
1.3791
1.2755
Exact [45]
1.4014
1.3813
1.3707
1.2138
S8R
1.4017
1.3816
1.3415
1.2158
QH8-R
1.4022
1.3827
1.3410
1.2133
SFSF
SSM
0.9523
0.9412
0.9146
0.8274
Shear-MK
0.9676
0.9455
0.9231
0.8575
MK
1.7644
1.6578
0.9344
0.8654
Exact [45]
0.95225
0.94314
0.92187
0.85011
S8R
0.9524
0.9433
0.9223
0.8325
QH8-R
0.9532
0.9442
0.9227
0.8519
Note. QH8-R (α0=-0.17). Shear correction factor κ=π2/12.Figure 15
Normalized uniaxial buckling load factors along thex-axis (σx0) versus the thickness-span aspect ratio t/b=0.001,0.05,0.1,0.2 with various boundaries.Through the achieved results shown in Table16 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, t/b=0.001, 0.05, but conversely matches well with thick plates, that is, t/b=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.The results of uniaxial buckling loads(σy0) along they-axis and the biaxial in-plane compression load (σx0=σy0) are shown in Tables 17 and 18 and Figures 16 and 17 associated with five different densities of regular elements such as 2×2 (21 nodes), 4×4 (65 nodes), 6×6 (133 nodes), 8×8 (225 nodes), and 10×10 (341 nodes).Table 17
The factor of uniaxial buckling loads (σy0) along the y-axis of rectangular plates (a/b=1) with various boundary conditions and various thickness-to-width ratios (t/b).
Boundary
Element
t
/
b
0.001
0.05
0.1
0.2
SSSS
2
×
2
4.7132
4.6261
4.3905
3.6771
4
×
4
4.1193
4.0431
3.8479
3.2455
6
×
6
4.0476
3.9751
3.7949
3.2231
8
×
8
4.0257
3.9564
3.785
3.2321
10
×
10
4.0161
3.9495
3.7829
3.2422
Exact [45]
4.0000
3.9437
3.7838
3.2558
SFSC
2
×
2
2.5834
2.5216
2.3678
1.9583
4
×
4
2.3993
2.3194
2.1530
1.7580
6
×
6
2.3821
2.2898
2.1191
1.7288
8
×
8
2.3797
2.2794
2.1079
1.7213
10
×
10
2.3799
2.2743
2.1031
1.7190
Exact [45]
2.3901
2.2667
2.1010
1.7200
SFSS
2
×
2
2.5626
2.502
2.3511
1.9480
4
×
4
2.3741
2.2965
2.1346
1.7483
6
×
6
2.3567
2.2672
2.1010
1.7195
8
×
8
2.3541
2.2569
2.0901
1.7122
10
×
10
2.3542
2.2518
2.0853
1.7099
Exact [45]
2.3639
2.2442
2.0829
1.7105
SFSF
2
×
2
2.1000
2.058
1.9553
1.6775
4
×
4
2.0327
1.9737
1.8516
1.5567
6
×
6
2.0286
1.9589
1.8315
1.5358
8
×
8
2.0295
1.9531
1.8244
1.5301
10
×
10
2.0309
1.9500
1.8212
1.5282
Exact [45]
2.0413
1.9457
1.8216
1.5333
Note. QH8-R (α0=-0.17). Shear correction factor κ=π2/12.Table 18
The factor of biaxial buckling loads (σx0=σy0) for rectangular plates (a/b=1) with various boundary conditions and various thickness-to-width ratios (t/b).
Boundary
Element
t
/
b
0.001
0.005
0.05
0.1
SSSS
2
×
2
2.3774
2.3770
2.3377
2.2277
4
×
4
2.0671
2.0667
2.0334
1.9432
6
×
6
2.0275
2.0272
1.9954
1.9108
8
×
8
2.015
2.0147
1.984
1.9023
10
×
10
2.0095
2.0091
1.9793
1.8992
Exact [45]a
2.0000
1.9997
1.9718
1.8919
SFSC
2
×
2
1.1820
1.1818
1.1652
1.1205
4
×
4
1.1474
1.1472
1.1261
1.0767
6
×
6
1.1427
1.1424
1.1184
1.0679
8
×
8
1.1416
1.1413
1.1155
1.0649
10
×
10
1.1414
1.1410
1.1140
1.0636
Exact [45]b
1.1431
1.1412
1.1119
1.0641
SFSS
2
×
2
1.0925
1.0923
1.0794
1.0439
4
×
4
1.0606
1.0605
1.0445
1.0060
6
×
6
1.0558
1.0556
1.0378
0.9986
8
×
8
1.0546
1.0543
1.0352
0.9962
10
×
10
1.0541
1.0538
1.0340
0.9952
Exact [45]b
1.0548
1.0535
1.0322
0.99541
SFSF
2
×
2
0.9671
0.9670
0.9583
0.9334
4
×
4
0.9400
0.9399
0.9302
0.9053
6
×
6
0.9351
0.9349
0.9247
0.9002
8
×
8
0.9335
0.9333
0.9228
0.8986
10
×
10
0.9328
0.9326
0.9219
0.8979
Exact [45]b
0.93209
0.9316
0.92071
0.89774
Note. QH8-R (α0=-0.17). aShear correction factor κ=0.823045. bShear correction factor κ=5/6.Figure 16
Normalized axial buckling load factors along they-axis (σy0) versus the thickness-span aspect ratio t/b=0.001, 0.05, 0.1, 0.2 with various boundaries.Figure 17
Normalized biaxial buckling load factors (σx0=σy0) versus the thickness-span aspect ratio t/b=0.001, 0.05, 0.1, 0.2 with various boundaries.We also consider simply supported plates (SSSS) with various thickness-to-width ratios,t/b = 0.05, 0.1, 0.2, and length-to-width ratios, a/b = 0.5, 1.0, 1.5, 2.0, 2.5.Table19 and Figure 18 show the buckling factors using the element QH8-R with regular mesh of 9×9 rectangular elements (280 nodes) when a/b=1.0. The QH8-R results are also compared to the DSG3 and ES-DSG3 [47] (mesh of 16×16 elements and 289 nodes), S8R (mesh of 9×9 elements and 280 nodes), mesh-free method [48], and the pb-2 Ritz [46]. It is seen that the QH8-R exhibits a good agreement with the pb-2 Ritz method, and S8R is not suitable for plate t/b = 0.2. Figure 19 also depicts the uniaxial buckling modes of simply supported rectangular plates with thickness-to- width ratios t/b = 0.01 and various length-to-width ratios, a/b = 1.0, 1.5, 2.0, 2.5.Table 19
The factor of uniaxial buckling loads along thex-axis of rectangular plates (SSSS) with various length-to-width ratios and various thickness-to-width ratios.
a
/
b
t
/
b
S8R
ES-DSG3 [47]
Mesh-free [48]
QH8-R
pb-2 Ritz [46]
0.5
0.05
6.0376
5.9873
6.0405
6.0822
6.0372
0.1
5.4779
5.3064
5.3116
5.5071
5.4777
0.2
1.9525
3.720
3.7157
3.9754
3.9963
1
0.05
3.9441
3.9412
3.9293
3.9528
3.9444
0.1
3.7861
3.7402
3.7270
3.7855
3.7865
0.2
2.0444
3.1263
3.1471
3.2427
3.2637
1.5
0.05
4.2576
4.2852
4.2116
4.2745
4.2570
0.1
4.0255
3.9844
3.8982
4.0332
4.0250
0.2
2.0132
3.1461
3.1032
3.2966
3.3048
2
0.05
3.9440
3.9811
3.8657
3.9528
3.9444
0.1
3.7860
3.7711
3.6797
3.7855
3.7865
0.2
1.9486
3.1415
3.0783
3.2427
3.2637
2.5
0.05
4.0673
4.1691
3.9600
4.0795
4.0645
0.1
3.8700
3.8924
3.7311
3.8737
3.8683
0.2
1.8697
3.1234
3.0306
3.2338
3.2421Figure 18
The factor of axial buckling loads of SSSS plates with various length-to-width ratios (a/b) and various thickness-to-width ratios (t/b).Figure 19
Axial buckling modes of simply supported rectangular plates with thickness-to-width ratiost/b = 0.01 and various length-to-width ratios a/b.
(a)
a
/
b
=
1.0
(b)
a
/
b
=
1.5
(c)
a
/
b
=
2.0
(d)
a
/
b
=
2.5


## 4.3.2. 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, t/b, 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.Figure 20
Skew plate under uniaxial in-plane loading.Numerical results for thick plates with thickness-to-width ratios,t/b = 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.Table 20
Buckling load intensity factor of skew plates havingS∗S∗S∗S∗, CCCC, CFCF, and S∗FS∗F boundary conditions subject to uniaxial loads (σx0).
Boundary
t
/
b
Skew angle
Mesh-free [48]
6
×
6
9
×
9
16
×
16
pb-2 Ritz [46]
S
∗
S
∗
S
∗
S
∗
0.05
0°
—
3.8057
3.7905
3.7815
3.7835
15°
—
4.1583
4.1380
4.1276
4.1314
30°
5.3874
5.4755
5.4304
5.4109
5.4182
45°
8.6926
8.9757
8.8002
8.7247
8.7328
0.1
0°
—
3.5111
3.5007
3.4965
3.4950
15°
—
3.8234
3.8094
3.8047
3.8027
30°
4.8534
4.9813
4.9475
4.9370
4.9324
45°
7.6106
7.9388
7.7931
7.7430
7.7236
CCCC
0.05
0°
9.5418
9.8091
9.6238
9.5647
9.5588
15°
10.2123
10.5549
10.3209
10.2430
10.2312
30°
12.5225
13.2413
12.7803
12.6114
12.5711
45°
17.8211
19.9900
18.6615
18.1211
17.9652
0.1
0°
8.1849
8.3903
8.2836
8.2936
8.2917
15°
8.7560
8.9338
8.7863
8.7845
8.7741
30°
10.3273
10.8235
10.4883
10.4238
10.3760
45°
13.2881
15.2116
14.2092
13.8638
13.6909
CFCF
0.05
0°
3.7984
3.8724
3.8281
3.8082
3.8007
15°
4.1310
4.2202
4.1688
4.1460
4.1387
30°
5.3510
5.4976
5.4102
5.3720
5.3660
45°
7.3658
7.9413
7.6341
7.4986
7.4670
0.1
0°
3.5063
3.5523
3.5262
3.5185
3.5077
15°
3.7684
3.8441
3.8135
3.8049
3.7937
30°
4.7773
4.8862
4.8316
4.8168
4.8043
45°
6.3092
6.6886
6.4547
6.3680
6.3311
S
∗
F
S
∗
F
0.05
0°
—
0.9459
0.9436
0.9424
0.9421
15°
—
1.0535
1.0513
1.0501
1.0501
30°
—
1.4612
1.4594
1.4583
1.4603
45°
—
2.5333
2.5317
2.5311
2.5450
0.1
0°
—
0.9227
0.9210
0.9202
0.9199
15°
—
1.0215
1.0199
1.0191
1.0189
30°
—
1.3892
1.3878
1.3873
1.3887
45°
—
2.3192
2.3187
2.3190
2.3316
Note. 6×6, 9×9, 16×16: the meshes for element QH8-R(-0.17).


## 4.3.3. 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 lengtha=b=10 and the thickness t=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 x-direction, σx0, pure shear load, and biaxial compressive loads in x and y directions, where σx0=σy0, respectively.Figure 21
Square plate with a circular hole and the distributions (308 nodes, 84 elements).The results of the normalized buckling load factor calculated for those loads with different boundary conditions are shown in Table21. The results agree very well with those computed by the EFG [57] involving both the third-order 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.Table 21
The uniaxial, biaxial, and pure shear buckling load factor for a square plate with a circular hole and different boundary conditions.
Buckling load
Boundary
EFG (FSDT)
EFG (TSDT)
Shear-MK
HQ8-R
Uniaxial
CCCC
7.995
8.097
8.0777
8.0543
SSSS
1.986
1.969
1.9544
2.0580
Biaxial
CCCC
4.781
4.860
4.9878
4.8477
SSSS
1.032
1.021
1.0511
1.0696
Pure shear
CCCC
12.669
13.357
13.5252
12.6976
SSSS
7.867
7.873
8.1544
7.8971Figure 22
The first uniaxial, biaxial, and pear shear buckling modes of the plate with a circular hole associated with different boundaries as SSSS, CCCC (thicknesst=1.0).
(a)
Uniaxial SSSS
(b)
Uniaxial CCCC
(c)
Biaxial SSSS
(d)
Biaxial CCCC
(e)
Pure shear SSSS
(f)
Pure shear CCCC


## 5. 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.



---

*Source: 6492081-2017-11-27.xml*