Mathematical Problems in Engineering
Volume 2017 (2017), Article ID 6492081, 33 pages
https://doi.org/10.1155/2017/6492081

Research Article

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

Tan Li,¹ Xu Ma,¹ Wanji Chen,² and Xili Jing¹

¹School of Science, Yanshan University, Qinhuangdao 066004, China
²Key Laboratory of Liaoning Province for Composite Structural Analysis of Aerocraft and Simulation, Shenyang Aerospace University, Shenyang 110136, China

Correspondence should be addressed to Xu Ma; maxu824@ysu.edu.cn

Received 3 July 2017; Revised 20 September 2017; Accepted 6 November 2017; Published 27 November 2017

Academic Editor: Giovanni Garcea

Copyright © 2017 Tan Li et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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 the 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 .

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 ). 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 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 ). 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 frame , with the origin on the mid-surface and the -axis in the thickness direction, , where 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

where , , are displacements along the -, -, and -axes, respectively, , are the rotations of the transverse normal about the - and -axes, and is the transverse displacement field.

The geometric equations can be written as follows:

where , , and are, respectively, the rotations, the curvatures, and the shear strains:

and operators , and are given by the following:

The equilibrium equations can be obtained from the strain energy in the following form:

where vectors and are, respectively, the moment and shear resultants:

The boundary forces can be written as follows:

where is the angle between the normal of edge and the local -axis of element.

For a linearly elastic material, the constitutive equations can be written as follows:

where and are the elasticity matrices of bending and transverse shear moduli. In the isotropic case, the elasticity matrices specialize as

where is Young’s modulus, the shear modulus, Poisson’s ration, and 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

where is the stress vector, is the elasticity matrices, is the vector of boundary force, is the boundary displacement vector, is the transverse displacement, and , are the rotations of the transverse normal about the - and -axes.

The approximation for stress and boundary displacements can now be incorporated in the functional. The stress field is described in the interior of the element as follows:

where is matrix of stress interpolation functions and is the unknown stress parameters.

The boundary force can be represented as follows:

where is the combination of direction cosine for the boundary normal.

The boundary displacement field is described by

where are interpolation functions and is nodal displacement parameters.

Substituting the stress equation (11), boundary force equation (12), and displacement approximations equation (13) into the functional (10),

where

The form of (15) is directly amenable to numerical integration (i.e., Gauss quadrature).

Then the internal strain energy can be expressed as follows:

By means of , we obtained

Consequently,

Substitution of in (16), the internal strain energy reduces to

Compared with , the element stiffness matrix can be taken as

The solution of the system yields the unknown nodal displacement . After is determined, element stress or internal forces can be recovered by use of (18) and (11). Thus

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:

where is the beam length, and are the values of the displacements and the rotations at the nodes equidistantly spaced between the beam ends, are the standard Lagrange polynomials of order , and for and otherwise, in which is the coordinate along the beam.

An 8-node quadrilateral element was designed as given in Figure 1. If any quadrilateral side is taken as a beam element, take the 1-2 boundary as example and the deflection and rotations , can be derived as follows:

where , , , , , , is the length of 1-2 boundary and is the coordinate along the 1-2 edge.

Figure 1: Eight-node quadrilateral plate element.

The displacement components along the 1-2 boundary was given as follows:

where

and , .

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 . The number of stress parameters, which is the number of columns in , must be at least equal to the number of degrees of freedom of the element less the number of degrees of freedom necessary to prevent rigid body motion.

In a series of studies of hybrid element, there is a general consensus on the selection of parameter that is , while the optimal selection is , , being the degrees of freedom of the element and the number of allowed rigid body motions, respectively. The element has 24 DOF; therefore, a stress field with at least 21 parameters is needed to describe the stress field and without spurious zero-energy modes. Numerical results show that the stiffness matrix has spurious zero-energy modes and converges slowly when . The finds that gradually increasing the number of until , there is a proper rank for the stiffness matrix and the absence of spurious zero-energy modes. Moreover, the calculation results are accurate and converged faster. And they are better if the number of reaches (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, are chosen as the assumed stress field.

The stress field is described in the interior of the element as follows:

where is matrix of stress interpolation functions as shown in Table 1, and is the stress parameters.

Table 1: Assumed stress function.

1
1
0
0
0
0

2
0
1
0
0
0

3
0
0
1
0
0

4

0
0
1
0

5

0
0
0
0

6
0

0
0
0

7
0

0
0
1

8
0
0

0
1

9
0
0

1
0

10

0

11

0
0

0

12

0
0
0
0

13
0

0
0
0

14
0

0
0

15
0

16
0
0

0

17
0
0

0

18

0

19

0

20

0
0

0

21

0
0
0
0

22
0

0
0
0

23
0

0
0

24
0

25
0

26
0
0

0

27
0
0

0

28

0

29

0

30

0
0

0

31

0
0
0
0

32
0

0
0
0

33
0

0
0

34
0

35

0

36
0
0

0

37
0

38

0

39
0
0

0

From (12), the boundary force can be expressed as

Then, along the boundary 1-2, can be obtained:

Similarly, along the other boundaries, , , 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:

where is the global stiffness matrix, is the global mass matrix, is the mode shapes, and is the natural frequencies.

The major complication of the standard eigenvalue problem introduced by the finite element method is the mass matrix . 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:

where 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 components , , as

where

Another 4-node shape functions are used to interpolate the deflection part of displacement function as

where

So the new deflection function part can be described as

where is an adjustable constant which is used to improve the accuracy of the vibration analysis.

Based on (31), (33), and (35), the shape function of refined consistent mass matrix is of the following form:

and the element consistent mass matrix can be expressed as

where is the matrix containing the mass density of the material and thickness as follows:

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

where is the global stiffness matrix, 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 can be written as

where the terms of (40) are

and vector stands for the initial stress of plate, , are unit tractions in the direction of and , and is the unit shear stress.

Similarly, we can obtain the refined geometric stiffness matrix from the following relation:

where is an adjustable constant used to improve the accuracy of the stability analysis, the displacement interpolation functions , , are the same as 8-node shape function (31), and is in terms of (33).

Based on (42), we only modify the interpolation function and geometric stiffness matrix part as

The refined geometric stiffness matrix can be written as

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 (), hard-type simply supported (), soft-type simply supported (), and free (), 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 , a free edge along , a clamped edge along , and a free edge along .

4.1. Patch Test: Consistency Assessment

Consistency of the developed elements is tested for the constant strain and stress states on the patch example with five elements, covering a rectangular domain of a plate as shown in Figure 3. The size of the domain is . The mechanical properties of the plate are chosen as , , .

Figure 3: Constant stress patch test.

For zero shear deformation patch test, the test function [20, 23] can be written as follows:

For nonzero shear deformation patch test, the test function can be expressed as follows:

where , are arbitrary constants. In the present paper, we assume that .

The mesh of patch test is presented in Figure 3. The displacements at the nodes on the boundaries are imposed according to (45) and (46), given the displacements and rotations at the boundary nodes (8 displacements and 16 rotations), while all the internal nodal displacements and rotations are to be calculated by the finite element solution procedure.

Test results show that the element QH8-39β passes both the constant bending with zero shear stresses , (Table 2) and the strict patch test with nonzero constant shear stresses , (Table 3). No matter the inner element edges are straight or curved, and no matter the shapes of the elements are convex or concave, the exact results of the displacements and stresses at inner node obtained using the present QH8-39β element agree well with the test functions. This demonstrates that the new element passes the strict constant stress patch test, thus ensuring solution convergence.

Table 2: Numerical results of constant bending patch test (at node ).

Mesh A

Mesh B

Mesh C

Exact

Mesh A

Mesh B

Mesh C

Exact

Table 3: Numerical results of nonzero constant shear stress patch test (at node ).

Mesh A

Mesh B

Mesh C

Exact

Mesh A

Mesh B

Mesh C

Exact

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 modulus , Poisson’s ratio , the mass density of the material , and the shear correction factor is taken if not specified otherwise.

In the numerical calculation, the dimensionless solutions of natural frequencies with a multiplier are obtained. is the width of the edge, . 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 and the width of the plate : ) is divided into mesh of elements. The thickness-to-width aspect ratios , boundary conditions SSSS, CCCC, and 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 of the present QH8-RCV for following numerical examples.

Table 4: The six lowest nondimensional frequency parameters of a CCCC thin square plate .

Element ()
Mode sequence number

1
2
3
4
5
6

QH8-RCV
3.6742
7.5726
7.5726
11.2791
13.8025
13.8626

QH8-RCV
3.6602
7.5052
7.5052
11.1242
13.5707
13.6325

QH8-RCV
3.6532
7.4711
7.4711
11.0458
13.4531
13.5157

QH8-RCV
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.395

Table 5: The six lowest nondimensional frequency parameters of a CCCC thick square plate .

Element ()
Mode sequence number

1
2
3
4
5
6

QH8-RCV
3.3118
6.3402
6.3402
8.9003
10.4844
10.5761

QH8-RCV
3.3066
6.3332
6.3332
8.9090
10.5178
10.6160

QH8-RCV
3.3002
6.3038
6.3038
8.8443
10.4237
10.5216

QH8-RCV
3.2943
6.2743
6.2743
8.7752
10.3198
10.4162

QH8-RCV
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.455

Table 6: The six lowest nondimensional frequency parameters of a SSSS thin square plate .

Element ()
Mode sequence number

1
2
3
4
5
6

QH8-RCV
2.0129
5.0805
5.0805
8.2063
10.3226
10.3226

QH8-RCV
2.0064
5.0403
5.0403
8.1034
10.1616
10.1616

QH8-RCV
2.0032
5.0202
5.0202
8.0522
10.0813
10.0813

QH8-RCV
2.0010
5.0066
5.0066
8.0174
10.0268
10.0268

QH8-RCV
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.000

Table 7: The six lowest nondimensional frequency parameters of a SSSS thick square plate .

Element ()
Mode sequence number

1
2
3
4
5
6

QH8-RCV
1.9424
4.6599
4.6599
7.1742
8.7513
8.7513

QH8-RCV
1.9376
4.6420
4.6420
7.1505
8.7325
8.7325

QH8-RCV
1.9345
4.6229
4.6229
7.1032
8.6600
8.6600

QH8-RCV
1.9321
4.6062
4.6062
7.0581
8.5881
8.5881

QH8-RCV
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.605

Figure 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 and QH8-LV are shown in Table 8 and Figure 5. The square plate is discretized by uniform meshes with . The results agree well with [27] when .

Table 8: Convergence of six lowest nondimensional frequency parameters (, 0.1 and SSSS, CCCC).

Boundary
Mesh
Mode sequence number

1
2
3
4
5
6

(a) QH8 with modified consistent mass matrix (QH8-RCV, )

0.001
SSSS

2.0208
5.1600
5.1600
8.5457
10.7782
10.7782

2.0044
5.0290
5.0290
8.0832
10.1251
10.1251

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

3.6986
7.7478
7.7478
11.8207
14.7803
14.9043

3.6563
7.4895
7.4895
11.1031
13.5317
13.5943

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

1.9440
4.6822
4.6822
7.3572
8.9607
8.9607

1.9335
4.6062
4.6062
7.0516
8.5663
8.5663

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

3.3008
6.3401
6.3401
9.0667
10.7861
10.8677

3.2925
6.2555
6.2555
8.7327
10.2492
10.3400

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

1.9997
4.9970
4.9970
7.9691
9.9694
9.9694

2.0000
4.9999
4.9999
7.9986
10.0001
10.0001

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

3.6452
7.4121
7.4121
10.8713
12.9885
13.0468

3.6461
7.4365
7.4365
10.9637
13.3312
13.3935

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

1.9427
4.6385
4.6385
7.0569
8.5511
8.5511

1.9452
4.6710
4.6710
7.1895
8.7670
8.7670

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

3.3019
6.2524
6.2524
8.6322
10.0169
10.0917

3.3155
6.3496
6.3496
8.9075
10.4886
10.5773

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.455

Figure 5: Convergence of six lowest nondimensional frequency parameters (, 0.1, and SSSS, CCCC).

Then, natural frequencies analyses are further carried on the thin and thick square plate with the meshes of 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

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.929

Figure 6: Shape of the six lowest modes of SSSS thin square plate obtained from QH8-RCV (, ).

4.2.2. Free Vibration of Circular Plate

Free vibration of a circular plate (Figure 7) of with the clamped boundary is then studied in this subsection. The material parameters of the circular plate are given as follows: the density mass , Young’s modulus , and Poisson’s ratio .

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 Figure 7, Mesh A, are used to discretize the whole plate for S8R of Abaqus, QH8-RCV 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 () and thick () 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 of a clamped circular plate with .

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.820

Table 11: Ten lowest parameterized natural frequencies of a clamped circular plate with .

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.409

(a) Thin circular plate

(b) Thick circular plate

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

First six lowest modes of this clamped circular obtained from QH8-RCV using the mesh of Figure 7 Mesh B are shown in Figure 9, which exactly describes the real physical modes of the circular plate.

Figure 9: Shape of the six lowest modes of a clamped thick circular plate obtained from QH8-RCV .

4.2.3. Free Vibration of Triangular Plates

QH8-RCV 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 , Young’s modulus , and Poisson’s ratio . The plate is discretized by the meshes as shown in Figure 10 for QH8. The convergence study of nondimensional frequency parameter for thin () and thick () triangular plates against different skew angles , 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 with .

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.8085

Table 13: The lowest six parameterized natural frequencies of triangular plates with .

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.7794

(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 nodes

Figure 10: Triangular plate model with various shape geometries and its discretization meshes (, 15°, 30°, 45°, 60°).

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

Figure 12: Convergence study of the lowest four parameterized natural frequencies of triangular plates with .

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 ().

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.

with being the edge width of the plate that is used throughout the paper and .

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

The parameter 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 () compressed in one direction or two perpendicular directions are divided into different meshes , , , , , . Thickness-to-width aspect ratios , 0.1, , 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() has higher accuracy when , and we advise for following numerical examples.

Table 14: The factor of uniaxial buckling loads along the -axis of rectangular plates with length-to-width aspect ratios and thickness-to-width aspect ratios .

Boundary condition
Element ()
Mesh

SSSS
QH8-R
3.1427
3.4982
3.7191
3.8183
3.8982
3.9293

4.0364
3.8502
3.9045
3.9338
3.9588
3.9691

4.7314
4.0776
4.0190
4.0043
3.9961
3.9940

4.8659
4.1389
4.0521
4.0258
4.0084
4.0027

4.7097
4.1161
4.0444
4.0225
4.0079
4.0031

4.2694
4.0139
4.0002
3.9980
3.9972
3.9971

Reference [44]
4

CCCC
QH8-R
5.3071
7.4051
8.3520
8.9070
9.4027
9.6049

8.2670
9.1692
9.4445
9.6292
9.7956
9.8651

10.9935
10.6774
10.2199
10.1074
10.0459
10.0304

12.9592
11.1920
10.4649
10.2587
10.1293
10.0885

13.6765
11.0371
10.4113
10.2349
10.1247
10.0896

13.6316
10.2573
10.0875
10.0576
10.0482
10.0472

Reference [44]
10.0738

Note. .

Table 15: The factor of uniaxial buckling loads along the -axis of rectangular plates with length-to-width aspect ratios and thickness-to-width aspect ratios .

Boundary
Element ()
Mesh

SSSS
QH8-R
2.3811
2.4760
2.7861
3.0392
3.3424
3.4914

3.3240
3.0601
3.2413
3.3983
3.5815
3.6689

4.1951
3.5827
3.6110
3.6676
3.7403
3.7756

4.4995
3.8245
3.7755
3.7778
3.7917
3.8002

4.3933
3.8501
3.7970
3.7869
3.7845
3.7850

3.9119
3.6966
3.7056
3.7146
3.7232
3.7266

Reference [45]
3.7864

CCCC
QH8-R
3.4293
4.0301
4.6232
5.2042
6.1652
6.7926

5.9999
5.8082
6.0193
6.4602
7.1856
7.6362

8.7108
7.7026
7.4918
7.6519
7.9938
8.2101

10.4843
8.8128
8.3089
8.2536
8.3146
8.3741

10.9492
8.8625
8.3903
8.2963
8.2811
8.2936

9.9368
7.5669
7.5101
7.5613
7.6923
7.7874

Reference [46]
8.2917

Note. .

Figure 14: The factor of uniaxial buckling loads () along the -axis of rectangular plates for QH8-R 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 ratios , 0.2 (moderately thick plate) and , 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 element (408 nodes) and MK, Shear-MK use a regular pattern of elements (441 nodes). A shear correction factor is used here for a proper comparison.

Table 16: The factor of uniaxial buckling loads () along the x-axis of rectangular plates () with various boundary conditions and various thickness-to-width ratios ().

Boundary
Method

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 . Shear correction factor .

Figure 15: Normalized uniaxial buckling load factors along the -axis () versus the thickness-span aspect ratio with various boundaries.

Through the achieved results shown in Table 16 and Figure 15, the shear-locking phenomenon, on the one side, is found very clearly because of losing the accuracy significantly when the standard MK is applied to deal with thin plates, that is, , 0.05, but conversely matches well with thick plates, that is, , 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 along the y-axis and the biaxial in-plane compression load () are shown in Tables 17 and 18 and Figures 16 and 17 associated with five different densities of regular elements such as (21 nodes), (65 nodes), (133 nodes), (225 nodes), and (341 nodes).

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

Boundary
Element

0.001
0.05
0.1
0.2

SSSS

4.7132
4.6261
4.3905
3.6771

4.1193
4.0431
3.8479
3.2455

4.0476
3.9751
3.7949
3.2231

4.0257
3.9564
3.785
3.2321

4.0161
3.9495
3.7829
3.2422

Exact [45]
4.0000
3.9437
3.7838
3.2558

SFSC

2.5834
2.5216
2.3678
1.9583

2.3993
2.3194
2.1530
1.7580

2.3821
2.2898
2.1191
1.7288

2.3797
2.2794
2.1079
1.7213

2.3799
2.2743
2.1031
1.7190

Exact [45]
2.3901
2.2667
2.1010
1.7200

SFSS

2.5626
2.502
2.3511
1.9480

2.3741
2.2965
2.1346
1.7483

2.3567
2.2672
2.1010
1.7195

2.3541
2.2569
2.0901
1.7122

2.3542
2.2518
2.0853
1.7099

Exact [45]
2.3639
2.2442
2.0829
1.7105

SFSF

2.1000
2.058
1.9553
1.6775

2.0327
1.9737
1.8516
1.5567

2.0286
1.9589
1.8315
1.5358

2.0295
1.9531
1.8244
1.5301

2.0309
1.9500
1.8212
1.5282

Exact [45]
2.0413
1.9457
1.8216
1.5333

Note. QH8-R (). Shear correction factor .

Table 18: The factor of biaxial buckling loads () for rectangular plates () with various boundary conditions and various thickness-to-width ratios ().

Boundary
Element

0.001
0.005
0.05
0.1

SSSS

2.3774
2.3770
2.3377
2.2277

2.0671
2.0667
2.0334
1.9432

2.0275
2.0272
1.9954
1.9108

2.015
2.0147
1.984
1.9023

2.0095
2.0091
1.9793
1.8992

Exact [45]
2.0000
1.9997
1.9718
1.8919

SFSC

1.1820
1.1818
1.1652
1.1205

1.1474
1.1472
1.1261
1.0767

1.1427
1.1424
1.1184
1.0679

1.1416
1.1413
1.1155
1.0649

1.1414
1.1410
1.1140
1.0636

Exact [45]
1.1431
1.1412
1.1119
1.0641

SFSS

1.0925
1.0923
1.0794
1.0439

1.0606
1.0605
1.0445
1.0060

1.0558
1.0556
1.0378
0.9986

1.0546
1.0543
1.0352
0.9962

1.0541
1.0538
1.0340
0.9952

Exact [45]
1.0548
1.0535
1.0322
0.99541

SFSF

0.9671
0.9670
0.9583
0.9334

0.9400
0.9399
0.9302
0.9053

0.9351
0.9349
0.9247
0.9002

0.9335
0.9333
0.9228
0.8986

0.9328
0.9326
0.9219
0.8979

Exact [45]
0.93209
0.9316
0.92071
0.89774

Note. QH8-R (). Shear correction factor . Shear correction factor .

Figure 16: Normalized axial buckling load factors along the -axis () versus the thickness-span aspect ratio , 0.05, 0.1, 0.2 with various boundaries.

Figure 17: Normalized biaxial buckling load factors () versus the thickness-span aspect ratio , 0.05, 0.1, 0.2 with various boundaries.

We also consider simply supported plates (SSSS) with various thickness-to-width ratios, = 0.05, 0.1, 0.2, and length-to-width ratios, = 0.5, 1.0, 1.5, 2.0, 2.5.

Table 19 and Figure 18 show the buckling factors using the element QH8-R with regular mesh of rectangular elements (280 nodes) when . The QH8-R results are also compared to the DSG3 and ES-DSG3 [47] (mesh of elements and 289 nodes), S8R (mesh of 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 = 0.2. Figure 19 also depicts the uniaxial buckling modes of simply supported rectangular plates with thickness-to- width ratios = 0.01 and various length-to-width ratios, = 1.0, 1.5, 2.0, 2.5.

Table 19: The factor of uniaxial buckling loads along the -axis of rectangular plates (SSSS) with various length-to-width ratios and various thickness-to-width ratios.

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.2421

Figure 18: The factor of axial buckling loads of SSSS plates with various length-to-width ratios () and various thickness-to-width ratios ().

(a)

(b)

(c)

(d)

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

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, , and different combinations of edge support conditions, are examined. The plate is modelled with , and distributed particles for QH8-R element.

Figure 20: Skew plate under uniaxial in-plane loading.

Numerical results for thick plates with thickness-to-width ratios, = 0.05, 0.1, skew angles, = 0°, 15°, 30°, 45°, and various boundary conditions of , CCCC, CFCF, and , are compared with those of mesh-free method [48] (with distributed particles) and Kitipornchai et al. [46] in Table 20. , , and 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 element.

Table 20: Buckling load intensity factor of skew plates having , CCCC, CFCF, and boundary conditions subject to uniaxial loads ().

Boundary

Skew angle
Mesh-free [48]

pb-2 Ritz [46]

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

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. , , : the meshes for element QH8-R.

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 length and the thickness are employed instead. The irregularly scattered pattern of 84 elements and 308 nodes shown in Figure 21 is used. The plate is subjected to three different in-plane load cases, namely, an uniaxial compressive load in -direction, , pure shear load, and biaxial compressive loads in and directions, where , 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 Table 21. 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.8971

(a) Uniaxial SSSS

(b) Uniaxial CCCC

(c) Biaxial SSSS

(d) Biaxial CCCC

(e) Pure shear SSSS

(f) Pure shear CCCC

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

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 , 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.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The work in this paper was supported by the National Natural Sciences Foundations of China (no. 11702242 and no. 11572204), Youth Foundation of Hebei Educational Committee (QN2017151), and Foundations of Yanshan University (no. 15LGB020, no. 16LGB015, no. B922, and no. B981). Those grants are gratefully acknowledged.

References

M. Wang, “On the necessity and sufficiency of the patch test for convergence of nonconforming finite elements,” SIAM Journal on Numerical Analysis, vol. 39, no. 2, pp. 363–384, 2001.

S. Cen and Y. Shang, “Developments of Mindlin-Reissner plate elements,” Mathematical Problems in Engineering, vol. 2015, Article ID 456740, 2015.

J. Robinson and G. W. Haggenmacher, “LORA-an accurate four node stress plate bending element,” International Journal for Numerical Methods in Engineering, vol. 14, no. 2, pp. 296–306, 1979.

E. Hinton and H. C. Huang, “A family of quadrilateral Mindlin plate elements with substitute shear strain fields,” Computers & Structures, vol. 23, no. 3, pp. 409–431, 1986.

O. C. Zienkiewicz, R. L. Taylor, and J. M. Too, “Reduced integration technique in general analysis of plates and shells,” International Journal for Numerical Methods in Engineering, vol. 3, no. 2, pp. 275–290, 1971.

E. D. L. Pugh, E. Hinton, and O. C. Zienkiewicz, “A study of quadrilateral plate bending elements with ‘reduced’ integration,” International Journal for Numerical Methods in Engineering, vol. 12, no. 7, pp. 1059–1079, 1978.

T. J. R. Hughes, M. Cohen, and M. Haroun, “Reduced and selective integration techniques in the finite element analysis of plates,” Nuclear Engineering and Design, vol. 46, no. 1, pp. 203–222, 1978.

T. J. R. Hughes and T. E. Tezduyar, “Finite elements based upon Mindlin plate theory with particular reference to the four-node bilinear isoparametric element,” Transactions ASME—Journal of Applied Mechanics, vol. 48, no. 3, pp. 587–596, 1981.

T. Belytschko, C. S. Tsay, and W. K. Liu, “A stabilization matrix for the bilinear Mindlin plate element,” Computer Methods Applied Mechanics and Engineering, vol. 29, no. 3, pp. 313–327, 1981.

J.-L. Batoz and P. Lardeur, “A discrete shear triangular nine d.o.f. element for the analysis of thick to very thin plates,” International Journal for Numerical Methods in Engineering, vol. 28, no. 3, pp. 533–560, 1989.

J.-L. Batoz and I. Katili, “On a simple triangular reissner/mindlin plate element based on incompatible modes and discrete constraints,” International Journal for Numerical Methods in Engineering, vol. 35, no. 8, pp. 1603–1632, 1992.

K.-J. Bathe, A. Iosilevich, and D. Chapelle, “Evaluation of the MITC shell elements,” Computers & Structures, vol. 75, no. 1, pp. 1–30, 2000.

A. K. Soh, Z. F. Long, and S. Cen, “A new nine dof triangular element for analysis of thick and thin plates,” Computational Mechanics, vol. 24, no. 5, pp. 408–417, 1999.

A.-K. Soh, S. Cen, Y.-Q. Long, and Z.-F. Long, “A new twelve {DOF} quadrilateral element for analysis of thick and thin plates,” European Journal of Mechanics - A/Solids, vol. 20, no. 2, pp. 299–326, 2001.

C. Wanji and Y. K. Cheung, “Refined quadrilateral element based on Mindlin/Reissner plate theory,” International Journal for Numerical Methods in Engineering, vol. 47, no. 1-3, pp. 605–627, 2000.

C. Wanji and Y. K. Cheung, “Refined 9-Dof triangular Mindlin plate elements,” International Journal for Numerical Methods in Engineering, vol. 51, no. 11, pp. 1259–1282, 2001.

S. Cen, Y.-Q. Long, Z.-H. Yao, and S.-P. Chiew, “Application of the quadrilateral area co-ordinate method: a new element for Mindlin-Reissner plate,” International Journal for Numerical Methods in Engineering, vol. 66, no. 1, pp. 1–45, 2006.

C. Wang, X. Wang, X. Zhang, and P. Hu, “Assumed stress quasi-conforming technique for static and free vibration analysis of Reissner-Mindlin plates,” International Journal for Numerical Methods in Engineering, 2017.

W. Chen, “Enhanced patch test of finite element methods,” Science China Physics, Mechanics & Astronomy, vol. 49, no. 2, pp. 213–227, 2006.

W. Chen, J. Wang, and J. Zhao, “Functions for patch test in finite element analysis of the Mindlin plate and the thin cylindrical shell,” Science China Physics, Mechanics & Astronomy, vol. 52, no. 5, pp. 762–767, 2009.

G. Jelenić and E. Papa, “Exact solution of 3D Timoshenko beam problem using linked interpolation of arbitrary order,” Archive of Applied Mechanics, vol. 81, no. 2, pp. 171–183, 2011.

T. Li, X. Ma, J. Xili, and W. Chen, “Higher-order hybrid stress triangular Mindlin plate element,” Computational Mechanics, vol. 58, no. 6, pp. 911–928, 2016.

T. Li, Z. Qi, X. Ma, and W. Chen, “Higher-order assumed stress quadrilateral element for the Mindlin plate bending problem,” Structural Engineering and Mechanics, vol. 54, no. 3, pp. 393–417, 2015.

W. J. Chen and Y. K. Cheung, “The Non-Conforming Nine-Parameters Triangular Element for Thin Plate Bending and Vibration Analysis,” in Proc. COCMSGE, P. K. K. Lee, L. G. Than, and Y. K. Cheung, Eds., pp. 1580–1583, Hong Kong, December 1994.

Y. K. Cheung and W. J. Chen, “Refined nine-parameter triangular thin plate bending element by using refined direct stiffness method,” International Journal for Numerical Methods in Engineering, vol. 38, no. 2, pp. 283–298, 1995.

X. Ma and W. Chen, “Refined 18-DOF triangular hybrid stress element for couple stress theory,” Finite Elements in Analysis and Design, vol. 75, pp. 8–18, 2013.

K. M. Liew, K. C. Hung, and M. K. Lim, “Vibration of mindlin plates using boundary characteristic orthogonal polynomials,” Journal of Sound and Vibration, vol. 182, no. 1, pp. 77–90, 1995.

T. H. H. Pian, “Derivation of element stiffness matrices by assumed stress distributions,” AIAA Journal, vol. 2, no. 7, pp. 1333–1336, 1964.

T. H. H. Pian, “Element stiffness matrices for boundary compatibility and for prescribed boundary stresses,” in Proceedings of the 1st Conf. on Matrix Methods in Struct, pp. 457–477, 1965.

D. S. Malkus and T. J. R. Hughes, “Mixed finite element methods—reduced and selective integration techniques: a unification of concepts,” Computer Methods Applied Mechanics and Engineering, vol. 15, no. 1, pp. 63–81, 1978.

P. Tong, “A family of hybrid plate elements,” International Journal for Numerical Methods in Engineering, vol. 18, no. 10, pp. 1455–1468, 1982.

K. Y. Sze, Q. H. Zhang, and G. H. Liu, “Hybrid quadrilateral finite element models for axial symmetric Helmholtz problem,” Finite Elements in Analysis and Design, vol. 52, pp. 1–10, 2012.

H. A. Santos and J. P. Moitinho de Almeida, “A family of Piola-Kirchhoff hybrid stress finite elements for two-dimensional linear elasticity,” Finite Elements in Analysis and Design, vol. 85, pp. 33–49, 2014.

R. E. Erkmen, “Shear deformable hybrid finite-element formulation for buckling analysis of thin-walled members,” Finite Elements in Analysis and Design, vol. 82, pp. 32–45, 2014.

T. Nishioka and S. N. Atluri, “Assumed stress finite element analysis of through-cracks in angle-ply laminates,” AIAA Journal, vol. 18, no. 9, pp. 1125–1132, 1980.

J. P. Wolf, “Structural averaging of stresses in the hybrid stress model,” Aiaa Journal, vol. 10, no. 6, pp. 843–845, 2015.

P. Tong, “New displacement hybrid finite element models for solid continua,” International Journal for Numerical Methods in Engineering, vol. 2, no. 1, pp. 73–83, 1970.

K.-J. Bathe and E. N. Dvorkin, “A four-node plate bending element based on Mindlin/Reissner plate theory and a mixed interpolation,” International Journal for Numerical Methods in Engineering, vol. 21, no. 2, pp. 367–383, 1985.

R. Ayad, G. Dhatt, and J. L. Batoz, “A new hybrid-mixed variational approach for Reissner-Mindlin plates, The MISP model,” International Journal for Numerical Methods in Engineering, vol. 42, no. 7, pp. 1149–1179, 1998.

S. Brasile, “An isostatic assumed stress triangular element for the Reissner-Mindlin plate-bending problem,” International Journal for Numerical Methods in Engineering, vol. 74, no. 6, pp. 971–995, 2008.

A. K. Soh, L. Zhifei, and C. Song, “A Mindlin plate triangular element with improved interpolation based on Timoshenko's beam theory,” Communications in Numerical Methods in Engineering, vol. 15, no. 7, pp. 527–532, 1999.

A. W. Leissa, Vibration of Plates, National Aeronautics and Space Administration, Washington DC, Wash, USA, 1969.

T. Irie, G. Yamada, and S. Aomura, “Natural frequencies of mindlin circular plates,” Journal of Applied Mechanics, vol. 47, no. 3, pp. 652–655, 1980.

S. P. Timoshenko and J. M. Gere, Theory of Elastic Stability, McGraw-Hill, New York, NY, USA, 1970.

S. Hosseini-Hashemi, K. Khorshidi, and M. Amabili, “Exact solution for linear buckling of rectangular Mindlin plates,” Journal of Sound and Vibration, vol. 315, no. 1-2, pp. 318–342, 2008.

S. Kitipornchai, Y. Xiang, C. M. Wang, and K. M. Liew, “Buckling of thick skew plates,” International Journal for Numerical Methods in Engineering, vol. 36, no. 8, pp. 1299–1310, 1993.

H. Nguyen-Xuan, G. R. Liu, C. Thai-Hoang, and T. Nguyen-Thoi, “An edge-based smoothed finite element method ({ES}-{FEM}) with stabilized discrete shear gap technique for analysis of Reissner-Mindlin plates,” Computer Methods Applied Mechanics and Engineering, vol. 199, no. 9-12, pp. 471–489, 2010.

K. M. Liew, J. Wang, T. Y. Ng, and M. J. Tan, “Free vibration and buckling analyses of shear-deformable plates based on FSDT meshfree method,” Journal of Sound and Vibration, vol. 276, no. 3-5, pp. 997–1017, 2004.

H. Sofuoglu and H. Gedikli, “A refined 5-node plate bending element based on Reissner-Mindlin theory,” Communications in Numerical Methods in Engineering, vol. 23, no. 5, pp. 385–403, 2007.

E. Hinton, T. Rock, and O. C. Zienkiewicz, “A note on mass lumping and related processes in the finite element method,” Earthquake Engineering & Structural Dynamics, vol. 4, no. 3, pp. 245–249, 1976.

F. Wu, G. R. Liu, G. Y. Li, A. G. Cheng, and Z. C. He, “A new hybrid smoothed FEM for static and free vibration analyses of Reissner-Mindlin Plates,” Computational Mechanics, 2014.

S. J. Lee, “Free vibration analysis of plates by using a four-node finite element formulated with assumed natural transverse shear strain,” Journal of Sound and Vibration, vol. 278, no. 3, pp. 657–684, 2004.

O. G. McGee, A. W. Leissa, and C. S. Huang, “Vibrations of cantilevered skewed trapezoidal and triangular plates with corner stress singularities,” International Journal of Mechanical Sciences, vol. 34, no. 1, pp. 63–84, 1992.

T. Mlzusawa, “Vibration of rectangular mindlin plates with tapered thickness by the spline strip method,” Computers & Structures, vol. 46, no. 3, pp. 451–463, 1993.

T. Q. Bui, M. N. Nguyen, and C. Zhang, “Buckling analysis of Reissner-Mindlin plates subjected to in-plane edge loads using a shear-locking-free and meshfree method,” Engineering Analysis with Boundary Elements, vol. 35, no. 9, pp. 1038–1053, 2011.

C. M. Wang, S. Kitipornchai, Y. Xiang, and K. M. Liew, “Stability of skew mindlin plates under isotropic in-plane pressure,” Journal of Engineering Mechanics, vol. 119, no. 2, pp. 393–401, 1993.

S. F. Li and W. K. Liu, Meshfree Particle Method, Springer, Berlin, Heidelberg, NY, 2004.