SciELO - Scientific Electronic Library Online

 
vol.14 issue5Crashworthiness Analysis of S-Shaped Structures Under Axial Impact LoadingAnalytical Bending and Stress Analysis of Variable Thickness FGM Auxetic Conical/Cylindrical Shells with General Tractions author indexsubject indexarticles search
Home Pagealphabetic serial listing  

Services on Demand

Journal

Article

Indicators

Related links

Share


Latin American Journal of Solids and Structures

Print version ISSN 1679-7817On-line version ISSN 1679-7825

Lat. Am. j. solids struct. vol.14 no.5 Rio de Janeiro June 2017

https://doi.org/10.1590/1679-78253036 

Articles

A New Triangular Hybrid Displacement Function Element for Static and Free Vibration Analyses of Mindlin-Reissner Plate

Jun-Bin Huanga 

Song Cena  b  * 

Yan Shanga  c 

Chen-Feng Lid 

aDepartment of Engineering Mechanics, School of Aerospace Engineering, Tsinghua University, Beijing 100084, China

bKey Laboratory of Applied Mechanics, School of Aerospace Engineering, Tsinghua University, Beijing 100084, China

cState Key Laboratory of Mechanics and Control of Mechanical Structures, College of Aerospace Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China

dZienkiewicz Centre for Computational Engineering & Energy Safety Research Institute, College of Engineering, Swansea University, Swansea SA2 8PP, UK


Abstract

A new 3-node triangular hybrid displacement function Mindlin-Reissner plate element is developed. Firstly, the modified variational functional of complementary energy for Mindlin-Reissner plate, which is eventually expressed by a so-called displacement function F, is proposed. Secondly, the locking-free formulae of Timoshenko’s beam theory are chosen as the deflection, rotation, and shear strain along each element boundary. Thirdly, seven fundamental analytical solutions of the displacement function F are selected as the trial functions for the assumed resultant fields, so that the assumed resultant fields satisfy all governing equations in advance. Finally, the element stiffness matrix of the new element, denoted by HDF-P3-7β, is derived from the modified principle of complementary energy. Together with the diagonal inertia matrix of the 3-node triangular isoparametric element, the proposed element is also successfully generalized to the free vibration problems. Numerical results show that the proposed element exhibits overall remarkable performance in all benchmark problems, especially in the free vibration analyses.

Keywords: Finite element method; Mindlin-Reissner plate element; hybrid displacement function element method; modified principle of complementary energy; static and free vibration analyses

1 INTRODUCTION

The Mindlin-Reissner plate theory is widely used to describe the deformation and resultant fields of an elastic plate subjected to transverse loads. As the rotations ψx , ψy , and deflection w are independently defined in this theory, only C 0 continuity is required for the compatible displacement fields. But it is found that the conventional Mindlin-Reissner plate bending elements with exact integration for computing the stiffness matrix will give poor results when the plate is quite thin, which is called as ‘shear locking’. In order to overcome the shortcoming, many effective techniques have been proposed. For example, Zienkiewicz et al. (1971) proposed the reduced integration technique, and Hughes et al. (1977) proposed the selective reduced integration technique. As a result of using inadequate Gauss points, these elements often present some spurious zero energy modes, and converge pretty slowly. To eliminate these spurious modes, some stabilization methods were introduced, such as the γ-methods proposed by Belytschko et al. (1986). Besides the methods mentioned above, other approaches that can improve the performance of Mindlin-Reissner plate elements include: the ‘Assumed Natural Strain’ (ANS) method proposed by Hughes et al. (1981), also called ‘Substituted Strains Methods’, in which the transversal shear strain fields are defined independently from the approximation of kinematic variables; the mixed interpolated tensorial components (MITC) family proposed by Bathe et al. (1985, 1989); the discrete Kirchhoff-Mindlin (DKM) element method proposed by Katili (1993); the discrete shear triangle (DST) method proposed by Batoz et al. (1989, 1992); the mixed shear projected (MiSP) method proposed by Ayad et al. (1998, 2001); the linked interpolation method proposed by Taylor et al. (1993); the improved shear strain interpolation schemes derived from the formulae of the locking-free Timoshenko’s beam element proposed by Soh et al. (1999); the refined Mindlin plate elements by Chen et al. (2001); the discrete shear gap (DSG) method proposed by Bletzinger et al. (2000); the smoothed finite element method proposed by Nguyen-Xuan et al. (2008); and so on (Cen and Shang, 2015).

Recently, based on a displacement function of Mindlin-Reissner plate introduced by Hu (1984) and the principle of minimum complementary energy, Cen et al. (2014) proposed a new finite element method called hybrid displacement function (HDF) element method. In their paper, a quadrilateral hybrid displacement function Mindlin-Reissner plate element HDF-P4-11β is formulated. Numerical examples show that this HDF element is free of shear locking, presents highly accurate results for both displacement and resultants with just several elements used. And what is more interesting is that this HDF element is insensitive to severe mesh distortion.

In this paper, based on the afore mentioned HDF element method, a 3-node, 9-DOF triangular Mindlin-Reissner plate element is formulated. Compared to Reference (Cen et al., 2014), a more rigorous but complicated description of the hybrid displacement function element method from the viewpoint of the principles in mechanics is given. Firstly, the modified variational functional of complementary energy for Mindlin-Reissner plate is discussed, and it can be eventually expressed by the displacement function F. Secondly, the locking-free formulae of Timoshenko’s beam element are chosen as the deflection, rotation, and shear strain along each element boundary. Thirdly, seven fundamental analytical solutions of the displacement function F are selected as the trial functions for the assumed resultant fields, so that the assumed resultant fields satisfy all governing equations in advance. Finally, the element stiffness matrix of the new element and the related equivalent nodal load are derived from the modified principle of complementary energy. The new element is denoted by HDF-P3-7β.

By employing similar but more complicated variational principles for free vibration problems, hybrid element method can be extended into the free vibration analyses of elastic structures (Tabarrok, 1971). But it should be noted that, even following such a complicated scheme, we cannot obtain the formulae of inertia matrix, and, which is more unexpected, a nonlinear eigenvalue problem has to be solved to give the vibration modes and frequencies. In this paper, we find that the stiffness matrix of the 3-node, 9-DOF triangular HDF Mindlin plate element can be used analogously as the stiffness matrices of the displacement-based elements. That is to say, a proper inertia matrix is found (here the diagonal inertia matrix of the 3-node triangular isoparametric element is a proper choice), and only linear eigenvalue problem for vibration modes and frequencies is solved. Thus, the proposed element can be successfully applied in free vibration problems following the simplified procedure, though no rigorous mathematical proof has been given yet.

Numerical results for the benchmark problems show that the presented element strictly passes the pure bending and twisting patch tests for both thin and thick plates, avoids shear locking and exhibits excellent performance for both displacement, resultants, and free vibration frequencies, and possesses better convergence than many other similar models.

2 THEORETICAL BASIS

2.1 The Mindlin-Reissner Plate Theory

As shown in Figure 1, the positive directions of displacement components and resultants are defined. The xoy-plane represents the middle surface of Mindlin plate, and the z-axis represents the direction of thickness (the thickness is denoted by h).

Figure 1: Positive directions of displacement components and resultants. 

In the mid-surface of Mindlin plate, the transverse deflection along z-axis is denoted by w, and the rotations of normal vector in xoz-plane and yoz-plane are denoted by ψx and ψy , respectively. Obviously, the deflection w, the rotation ψx , and the rotation ψy are all functions of x and y. And the displacement components at arbitrary point in the plate are given by equation (1).

(1)

As the deflection and rotations are independently defined, the rotations ψx and ψy are no more equal to the rotations of mid-surface. The transverse shear strain vector is determined by the difference of the two aforementioned kinds of rotations.

(2)

From the displacements given in equation (1), the strains paralleled with the xoy-plane at arbitrary point in Mindlin plate can be obtained:

(3)

where κ is the curvature vector. And from equations (2) and (3), the strain compatibility equations can be derived:

(4)

The bending moments, twisting moment, and shear forces are denoted by Mx , My , Mxy , Tx , and Ty , respectively, as shown in Figure 1. And the equilibrium equations for a plate under transverse distributed load q are:

(5)

As shown in Figure 2, Mn , Ms , and Tn denote the bending moment, twisting moment, and shear force along boundary, respectively; ψn , ψs , and w denote the boundary rotations and deflection, respectively.

Figure 2: Positive directions of boundary resultants and displacement components. 

The relations between the domain and the boundary resultants, as well as the relations between the domain and the boundary displacement components, are given by

(6)

(7)

where l and m are the direction cosines of the boundary’s outer normal.

The constitutive equations for isotropic and linearly elastic Mindlin-Reissner plates with uniform thickness are:

(8)

where M, T, D b , and D s are the bending moment vector, shear force vector, bending elasticity matrix, and shear elasticity matrix, respectively;

(9)

(10)

(11)

in which E, G, and μ are the Young’s modulus, shear modulus, and Poisson’s ratio, respectively; and h is the thickness of the plate.

Generally, the boundary conditions (BCs) are classified into four categories: 1) Fixed boundary (S 1) condition: deflection and rotations of the boundary are given; 2) ‘Soft’ simply supported (SS 1) condition: deflection, bending moment and twisting moment of the boundary are given; 3) ‘Hard’ simply supported (SS 2) condition: deflection, bending moment and rotation ψs of the boundary are given; and 4) Free boundary (S 3) condition: all the boundary resultants are given.

2.2 The Modified Principle of Complementary Energy for Mindlin-Reissner Plate

The standard complementary energy functional for Mindlin-Reissner plate is given below (Hu, 1984):

(12)

in which , , and denote the given boundary displacement components on the displacement boundary Su ; Ω denotes the whole mid-surface region of the plate; C and R denote the elasticity matrix of compliances and the resultant vector, respectively;

(13)

(14)

Assume that the region Ω is divided into several subregions (triangles) Te (e=1 ˜ N). And the given trial functions of resultants satisfy the equilibrium equations in each subregion Te . In order to make use of the principle of minimum complementary energy, extra constrains should be imposed. According to the generalized variational principle, such extra constrains can be imposed by the so-called Lagrange multiplier method.

Extra constrains include the equilibrium conditions on the common boundaries of subregions and the force boundary conditions on the boundary of Ω, which are expressed by:

(15)

where Σ Se denotes the common boundaries of subregions; the indexes 1 and 2 denote the two subregions that share a common boundary;

(16)

where denotes the force boundary; , , and represent the given boundary resultants.

By multiplying the two constrains by corresponding Lagrange multipliers and adding them into the standard complementary energy functional (equation (12)), the modified functional can be obtained:

(17)

where λ 1 and λ 2 are Lagrange multipliers.

Let the variation of equation (17) be zero:

(18)

Here the relations (6) between domain and boundary resultants have been already considered.

Considering the following constitutive equation:

(19)

where F is an operator matrix, the area integral in equation (18) can be simplified. Integrate by parts for the area integral, we obtain

(20)

in which is the adjoint operator of F; n is the corresponding direction matrix; and ∂Te denotes the boundary of the e-th subregion;

(21)

where l and m are still the direction cosines of boundary’s outer normal.

Because the trial functions of resultants satisfy the equilibrium equations in each subregion, variation of equations (5) yields following result:

(22)

Thus, equation (18) can be simplified as

(23)

In above equation, only curvilinear integral exists. From equation (7), and using the relation

(24)

equation (23) can be further simplified to be

(23*)

Due to the arbitrariness of variation about resultants δ R, it can be easily concluded that the Lagrange multipliers are just (boundary) displacement components on corresponding boundaries. Finally, the modified functional of complementary energy can be fully determined:

(25)

2.3 The Displacement Function of Mindlin-Reissner Plate

Hu (1984) proposed the so-called displacement functions F and f for Mindlin-Reissner plates, from which the displacement components can be derived theoretically.

For a plate loaded by distributed transverse force q, the displacement components are given as

(26)

where D and C are bending stiffness and shear stiffness mentioned in equation (11), respectively. The functions F and f satisfy the following differential equations.

(27)

(28)

It should be noted that equation (27) is quite similar to the governing equation of the deflection for thin plates under Kirchhoff’s assumption. In fact, Hu pointed out that function F even satisfies the same boundary conditions as the deflection of thin plates if the plate is thin enough. Function f represents the influence of shear deformation, and describes the phenomenon of edge effect for Mindlin plates. In some cases, the function f vanishes, such as simply supported polygonal plates, circular plates of axial symmetry and other cases in which the shear forces are statically determinate (Hu, 1984). In this paper, only F is considered, because the influence region of f is near the boundary, and the contribution of f to the total energy and the displacement solutions is insignificant. Actually, in higher-order elements, the influence of displacement function f must be taken into account, otherwise the elements may present poor results for stress and strain, even for nodes away from the boundary (Bao et al., 2017).

Substitution of equation (26) into equations (2) and (3) yields

(29)

Then, substitution of equation (29) into the constitutive equations (8) yields

(30)

It is emphasized that the curvatures, shear strains, and resultants derived from displacement function F satisfy all governing equations of Mindlin-Reissner plates. Thus, it’s quite reasonable to gain proper trial functions by means of displacement functions.

Assumed that the distributed transverse force q is constant, then, the displacement function F can be solved.

(31)

in which F is expressed as the sum of the general and the particular solutions.

Considering the symmetry between x and y axes, the particular solution can be chosen as (other choices such as F * = q / 8D x 2 y 2 have been also tested, but the following one is better.)

(32)

And the general solutions of F satisfy the following governing equation,

(33)

which can be easily solved in polynomial form. The first seven general solutions and corresponding resultants are listed in Table 1. The first order completeness of bending moments and zero order completeness of shear forces are guaranteed. More analytical solutions can be found in reference (Cen et al., 2014).

Table 1: The first seven general solutions of F and corresponding resultants. 

i 1 2 3 4 5 6 7
-DF0i x 2 xy y 2 x 3 x 2 y xy 2 y 3
M 0 xi Ri 2 0 2μ 6x 2y 2μx 6μy
M 0 yi 2μ 0 2 6μx 2μy 2x 6y
M 0 xyi 0 1−μ 0 0 2(1−μ)x 2(1−μ)y 0
T 0 xi 0 0 0 6 0 2 0
T 0 yi 0 0 0 0 2 0 6

2.4 The Locking-free Formulae of Timoshenko’s Beam

As shown in Figure 3, in locking-free formulae of Timoshenko’s beam (Hu, 1984), the deflection w, rotation ψ, and shear strain γ are assumed to be cubic, quadratic, and constant, respectively, and are given by

(34)

in which

(35)

Figure 3: Positive directions of displacement components for Timoshenko’s beam. 

Here D and C represent the beam’s bending stiffness and shear stiffness, respectively, and will be replaced by corresponding stiffness of Mindlin plates in the following section.

As discussed in Section 1, the interpolation technique based on the formulae of Timoshenko’s beam has been successfully applied in many effective elements (Soh et al., 1999).

3 FORMULAE OF THE NEW ELEMENT HDF-P3-7Β

Figure 4 shows a 9-DOF triangular plate bending element, and assume that the whole plate Ω is meshed into several triangles Te (e=1˜N).

Figure 4: Nodal DOFs of the 3-node triangular hybrid displacement function element. 

Firstly, the assumed resultant fields can be derived from the aforementioned solutions of displacement function F:

(36)

in which R 0 denotes the general solution part, and is taken as a linear combination of the first seven general solutions of F listed in Table 1:

(37)

where βi (i=1˜7) are unknown resultant parameters, and R i (i=1˜7) are the resultant vectors derived from the first seven general solutions of F by employing equation (30), which have been given by Table 1; R * denotes the particular solution part, and can be derived directly from the particular solution of F and equation (30):

(38)

Thus, the assumed resultants satisfy the equilibrium equations in each subregion Te . So, the modified functional of complementary energy given in Section 2 can be employed.

Secondly, the displacement components along each element boundary are chosen as the locking-free formulae of Timoshenko’s beam theory.

Substitution of the relations

(39)

into equation (34) yields

(40)

Then, the rotation along normal direction is assumed to be linear function.

(41)

These three complicated equations can be rewritten in the form of matrix:

(42)

where q e is the nodal displacement vector of the e-th element, and is the boundary displacement interpolation function matrix or shape function matrix.

It is easy to see that the matrices on different edges are different:

(43)

All the non-zero components can be detailed directly from equations (40) and (41).

Thirdly, as the assumed resultants and boundary displacement components are all given, the modified functional of complementary energy can be fully determined.

Let

(44)

then, the modified functional of complementary energy can be simplified as:

(45)

in which β e denotes the resultant parameter vector of the e-th element. It should be emphasized that the resultant parameter vectors of different elements are independently defined, so that they can vary independently.

Using the modified principle of complementary energy, the variation of modified functional (45) is taken to be zero:

(46)

in which q is the global nodal displacement vector of the plate. And above variational equations yield

(47)

(48)

which gives the relation between β e and q e . It may be found that equation (48) is not so rigorous, and some assembling rules should be discussed.

Equation (47) yields

(49)

Substitution of equation (49) into equation (48) yields

(50)

Above equation gives the discrete global equilibrium equation for this kind of element:

(51)

in which the summation notation means that some assembling rules of element equilibrium equations should be followed. In fact, the assembling rule is the same as what we use in traditional displacement-based element cases, which can be easily found from equation (50).

It should be noted that:

The element stiffness matrix is just H T M -1 H;

(52)

Now it can be explained that why only the first seven general solutions of resultants are employed. To avoid spurious zero energy modes, the rank of single element stiffness matrix should be at least 9−3=6, which means that at least six general solutions should be included. Completeness requirement tells us that seven solutions or eleven solutions are reasonable choices, but numerical results show that seven solutions are enough and always achieve higher accuracy.

2) The element equivalent nodal force vector is also given in equation (51);

(53)

In above equation, (V T - H T M -1 M *) denotes the equivalent nodal forces of body forces, and the last part denotes the equivalent forces of boundary forces.

Then, the element equilibrium equation can be written as P e - K e q e = 0, and equation (51) is the assembled global equilibrium equation. It is surprising but reasonable that the equivalent forces of boundary forces possess the same form as displacement-based elements.

3) The assembling rules of hybrid displacement function elements are consistent with traditional displacement-based elements. So it can be easily integrated into the standard framework of finite element programs.

The integrals in equation (44) are all calculated by numerical integration methods. One dimensional Gauss integration and Hammer integration schemes are employed. As the integrands are polynomials, all the integrals can be calculated exactly. In this paper, four integration points are employed for both line integrals and area integrals in equation (44), and it can be sure that the results are all accurate.

Following the aforementioned steps, the global nodal displacement components can be solved. Then, substitute the known element displacement vector q e into

(54)

the resultants at arbitrary point of the element can be obtained. Furthermore, the resultants at any nodes can be taken as the average value of different elements.

Once the stiffness matrix for static analysis is derived, what we need to extend the proposed element into free vibration analysis of Mindlin plate is the formulae of inertia matrix. Considering the corresponding variational principles for free vibration problems, hybrid element method can be extended into the analysis of free vibration problems theoretically. But as we mentioned in Section 1, such scheme does not give the formulae of inertia matrix and a nonlinear eigenvalue problem has to be solved.

Actually, some inertia matrices of displacement-based elements can cooperate well with the stiffness matrix of the proposed element in analyzing the free vibration of plate structures. And the diagonal inertia matrix of the 3-node triangular isoparametric element is almost the best choice, for both convergence property and computational cost. Thus, though no rigorous mathematical proof has been given yet, the generalized eigenvalue equation is still given by

(55)

where 𝒬 is the density of the plate, A is the area of the e-th element and ω is the natural frequency. Since the stiffness matrix and inertia matrix are symmetric, this eigenvalue problem can be solved by the subspace iteration method (Bathe and Wilson, 1973), and the frontal method is also employed in computer coding for saving computer memory (Liu, 1989).

4 NUMERICAL EXAMPLES

Several standard numerical examples are employed in this section to test the performance of the proposed element HDF-P3-7β. And results obtained by following models are also given for comparison.

DKT: Triangular discrete Kirchhoff element proposed by Batoz et al. (1980).

DKTM, RDKTM: Refined triangular Mindlin plate elements proposed by Chen et al. (2001).

DST-BL: A compatible triangular Mindlin plate element based on the discrete shear triangle technique proposed by Batoz et al. (1989).

DST-BK: An incompatible triangular Mindlin plate element based on the discrete shear triangle technique proposed by Batoz et al. (1992).

MITC4: A quadrilateral Mindlin-Reissner plate element based on a mixed interpolation scheme proposed by Bathe et al. (1985).

THS: A triangular hybrid stress element based on analytical solutions of thin plate equations proposed by Rezaiee-Pajand et al. (2014).

ARS-T9: A triangular Mindlin plate element based on the improved shear strain interpolation derived from the locking-free Timoshenko’s beam element proposed by Soh et al. (1999).

T3BL, T3BL(R): Triangular Mindlin plate elements based on the linked interpolation method proposed by Taylor and Auricchio (1993).

MiSP3: Hybrid-mixed Mindlin plate elements based on the mixed shear projected method proposed by Ayad et al. (1998).

MiSP3+: An improved triangular hybrid-mixed Mindlin plate element proposed by Ayad et al. (2001).

NS+ES-FEM: A ‘hybrid’ smoothed element proposed by Wu et al. (2014).

MIN3: A 3-node Mindlin plate element with improved transverse shear proposed by Tessler et al. (1985).

DSG3: Mindlin plate element based on the discrete shear gap method proposed by Bletzinger et al. (2000).

NS-DSG3, ES-DSG3: Two smoothed Mindlin plate elements of the smoothed FEM family proposed by Liu et al. (2009a, 2009b).

S4R, S3R: The quadrilateral and triangular shell elements assembled in the renowned commercial FEM software Abaqus (2009).

4.1 Numerical Examples: Static Analyses

4.1.1 Patch Tests

As shown in Figure 5, a patch is divided by four elements. Three kinds of thickness are considered. Proper constrains are imposed to eliminate rigid body motions. The size and constants are also given in the figure. Both pure bending boundary forces and pure twisting boundary forces are tested, and these two cases are shown in Figure 6.

Figure 5: Patch tests: Geometry and mesh type. 

Figure 6: Pure twisting and pure bending boundary conditions. 

Numerical results in Table 2 show that the patch tests are perfectly passed for both thin and moderately thick plates, so that the convergence should be guaranteed. In fact, as these two constant moment cases are all included in the first seven general solutions of resultants, which are the trial functions for the resultant fields, it is reasonable that exact solutions can be obtained.

Table 2: Results for patch tests. 

Boundary Conditions Pure Twisting Pure Bending
Nodal No. Mx My Mxy Mx My Mxy
h =0.04
1 1.9664E-08 2.3121E-08 1.0000E+00 1.0000E+00 1.0000E+00 1.3626E-08
2 4.9135E-07 -1.1536E-07 1.0000E+00 1.0000E+00 1.0000E+00 -3.5914E-08
3 -9.6798E-08 2.0712E-08 1.0000E+00 1.0000E+00 1.0000E+00 4.2495E-09
4 -4.5968E-08 -1.5839E-08 1.0000E+00 1.0000E+00 1.0000E+00 -7.9921E-08
5 -2.5876E-07 5.5718E-08 1.0000E+00 1.0000E+00 1.0000E+00 3.6780E-08
h =0.4
1 -1.9714E-08 2.3066E-08 1.0000E+00 1.0000E+00 1.0000E+00 1.3630E-08
2 4.9080E-07 -1.1592E-07 1.0000E+00 1.0000E+00 1.0000E+00 -3.5973E-08
3 -9.6580E-08 2.0811E-08 1.0000E+00 1.0000E+00 1.0000E+00 4.2770E-09
4 -4.5956E-08 -1.5694E-08 1.0000E+00 1.0000E+00 1.0000E+00 -7.9809E-08
5 -2.5854E-07 5.6058E-08 1.0000E+00 1.0000E+00 1.0000E+00 3.6731E-08
h =4
1 -2.3707E-08 1.9313E-08 1.0000E+00 1.0000E+00 1.0000E+00 1.4086E-08
2 4.4927E-07 -1.5838E-07 1.0000E+00 1.0000E+00 1.0000E+00 -3.8779E-08
3 -8.0601E-08 2.8371E-08 1.0000E+00 1.0000E+00 1.0000E+00 6.2446E-09
4 -4.4667E-08 -4.3771E-09 1.0000E+00 1.0000E+00 1.0000E+00 -7.2774E-08
5 -2.4020E-07 8.1505E-08 1.0000E+00 1.0000E+00 1.0000E+00 3.2917E-08

4.1.2 Square Plate Loaded by Uniform Distributed Transverse Load

As shown in Figure 7, due to the symmetry, only a quarter of plate is calculated. Clamped boundary condition and hard simply supported boundary condition are considered. The span, distributed transverse force, bending stiffness, Poisson’s ratio, and thickness of the plate are denoted by L, q, D, μ, and h, respectively, and are given by

L=1, q=1, D=1, μ=0.3, h=0.001, 0.1.

Figure 7: Typical mesh employed for a quarter of plate. 

The deflection and bending moment at the plate center

(56)

calculated by the proposed element, which are compared with the reference solutions (Taylor and Auricchio, 1993), are listed in Table 3. Comparisons with results obtained by other elements are plotted in Figure 8 and Figure 9. From these results, it can be concluded that the proposed HDF-P3-7β element shows better convergence in most of square plate cases.

Table 3: Normalized central deflection and bending moment results of the proposed element. 

BCs Thickness 2×2 4×4 8×8 16×16 Reference
SS 2 Deflection
h=0.001 0.9876 0.9974 0.9994 0.9999 1.0000
h=0.1 0.9914 0.9971 0.9991 0.9998
Bending Moment
h=0.001 0.9569 0.9815 0.9931 0.9975 1.0000
h=0.1 0.9568 0.9867 0.9962 0.9988
S 1 Deflection
h=0.001 0.9681 0.9915 0.9987 1.0000 1.0000
h=0.1 0.9970 1.0019 1.0017 1.0031
Bending Moment
h=0.001 1.0768 0.9937 0.9959 0.9982 1.0000
h=0.1 1.0191 1.0070 1.0046 1.0043

Figure 8: Hard simply supported boundary case: the plate central deflection and bending moment results. 

Figure 9: Clamped boundary case: the plate central deflection and bending moment results. 

4.1.3 Skew Plates Loaded by Uniform Distributed Transverse Load

Razzaque’s 60° plate

As shown in Figure 10, a 60° skew plate with two free and two soft simply supported edges was firstly used by Razzaque (1973) to test the accuracy of thin plate elements, and it has been treated as a classical benchmark for both thin and thick plate elements. The edge length, transverse force, thickness, Young’ modulus, and Poisson’s ratio are denoted by L, q, h, E, and μ, respectively, and are given by

L=100, q=1, h=0.1, E=10.92, μ=0.3.

Figure 10: A 60° skew plate: Geometry and mesh type. 

For the central deflection w C and bending moment My at point C, Razzaque (1973) give out the finite difference solutions:

(57)

Results obtained by different elements are listed in Table 4 and plotted partly in Figure 11. It can be easily concluded that the proposed element gives almost the best solutions for both deflection and bending moment.

Table 4: 60° skew plate: Normalized deflection and bending moment at point C. 

Mesh 2×2 4×4 8×8 16×16 Reference
wC/w ref HDF-P3-7 β 0.9050 0.9736 0.9899 0.9941 1.0000
DKT 0.8045 0.9474 0.9845 0.9947
ARS-T9 0.7964 0.9488 0.9838 0.9926
MiSP3 0.7218 0.9314 0.9824 0.9936
MiSP3+ 0.7823 0.9578 0.9893 0.9942
MITC4 0.4853 0.8419 0.9515 0.9814
DKTM - 0.9474 0.9919 -
RDKTM - 0.9474 0.9919 -
THS 0.9093 0.9692 0.9878 -
My /M ref HDF-P3-7 β 0.9305 0.9913 0.9999 1.0008
DKT 0.9336 1.0000 1.0007 1.0009
ARS-T9 0.9537 0.9969 1.0012 1.0012
MiSP3 0.6954 0.9313 0.9974 0.9981
MiSP3+ 0.7057 0.9490 0.9908 0.9992
MITC4 0.3940 0.8036 0.9459 0.9878

Figure 11: Results of Razzaque’s 60° plate. 

Morley’s 30° plate

As shown in Figure 12, a 30° plate with four soft simply supported edges is considered. It’s a more singular problem than the previous one. The earliest solution is calculated based on Kirchhoff theory by Morley (1963).

Figure 12: A 30° plate: Geometry and mesh type. 

The edge length, thickness, Young’s modulus, and Poisson’s ratio are denoted by L, h, E, and μ, respectively, and are given by

L=100, h=0.1, E=10.92, μ=0.3.

Then, the central deflection, maximum, and minimum bending moments at point O are evaluated. Morley’s results are employed as the reference solutions.

(58)

Results of central deflection and two principal bending moments are listed in Table 5, and plotted in Figure 13.

Table 5: 30° skew plate: Normalized deflection and principal bending moments at point O. 

Mesh 2×2 4×4 8×8 16×16 Reference
wO/(qL 4/1000D) HDF-P3-7 β 1.4473 1.1103 1.0348 1.0225 1.0000
NS+ES-FEM - 1.2049 1.0907 1.0554
ARS-T9 1.5338 1.1096 1.0397 1.0277
MIN3 - 0.6681 0.7600 0.8493
DSG3 - 0.6400 0.6507 0.7718
T3BL 1.3164 1.0316 1.0179 1.0137
T3BL(R) 1.7939 1.1978 1.0752 1.0449
DKT - 1.1103 1.0392 1.0270
DKTM - 0.8750 0.8309 0.8652
RDKTM - 1.1103 1.0392 1.0270
THS - 1.0368 1.0000 1.0000
M max/(qL 2/100) HDF-P3-7 β 1.0914 1.0727 0.9971 1.0113
NS+ES-FEM - 1.0397 1.0397 1.0185
ARS-T9 1.4763 1.1605 1.0382 1.0215
MIN3 - 0.6542 0.7970 0.8855
DSG3 - 0.6499 0.6840 0.7949
T3BL 0.7540 0.9019 0.9910 1.0012
T3BL(R) 0.7579 0.9528 1.0012 1.0171
M min/(qL 2/100) HDF-P3-7 β 0.5742 1.0953 0.9328 1.0213
ARS-T9 1.7886 1.3234 1.0331 1.0407
T3BL 0.6937 0.8846 1.0097 1.0188
T3BL(R) 0.7589 0.8931 1.0177 1.0415

Figure 13: 30° skew plate: Deflection and principal bending moments at point O. 

It can be seen that good convergence is achieved nearly by all the elements, but the proposed element provides higher accuracy for bending moments.

4.1.4 Circular Plate Under Uniform Distributed Transverse Load

Similar to the square plate problem, a quarter of circular plate is investigated and symmetric boundary conditions are imposed on corresponding edges, as shown in Figure 14.

Figure 14: Circular plate: Geometry and mesh type. 

Soft simply supported (SS 1) boundary condition is applied to the plates. For this axial symmetric problem, analytical solutions can be derived from the Mindlin-Reissner plate theory (Batoz and Dhatt, 1990). The radius, uniform distributed transverse load, thickness, Young’s modulus, and Poisson’s ratio are denoted by R, q, h, E, and μ, respectively.

The central deflection and bending moment of soft simply supported circular plates are given as

(59)

Related parameters are set to be

R=5, h=0.1, 1, E=10.92, μ=0.3.

The normalized results of the deflection and the bending moment at plate center are given in Table 6. And corresponding data are also plotted in Figure 15.

Table 6: Soft simply supported circular plate: Normalized central deflection and bending moment. 

  Number of Elements (N) 6 24 96 384 Reference
wc /w ref h=0.1 HDF-P3-7 β 1.0282 1.0074 1.0019 1.0005 1.0000
ARS-T9 0.9502 0.9894 0.9991 0.9993
T3BL(R) 1.0599 1.0158 1.0028 0.9993
T3BL 0.9900 0.9931 0.9988 0.9997
DST-BK 0.9387 0.9850 0.9962 -
DST-BL 0.9502 0.9891 0.9974 -
DKT 0.9498 0.9886 0.9971 -
DKTM 0.7877 0.9611 0.9935 -
RDKTM 0.9502 0.9890 0.9975 -
h=1 HDF-P3-7 β 1.0233 1.0058 1.0014 1.0003
ARS-T9 0.9486 0.9879 0.9991 0.9992
T3BL(R) 1.0554 1.0150 1.0038 1.0010
T3BL 0.9895 0.9959 0.9989 0.9998
DST-BK 0.9348 0.9838 0.9965 -
DST-BL 0.9900 0.9985 1.0025 -
DKTM 0.9219 0.9827 0.9961 -
RDKTM 0.9486 0.9877 0.9970 -
Mc /M ref Number of Elements (N) 6 24 96 384
h=0.1 HDF-P3-7 β 1.0294 1.0086 1.0029 1.0008
ARS-T9 1.0205 1.0199 1.0077 1.0021
TL(R) 0.9209 0.9771 0.9933 0.9976
T3BL 0.9131 0.9757 0.9938 0.9985
DST-BK 1.0531 1.0240 1.0085 -
DST-BL 1.0201 1.0143 1.0046 -
DKT 1.0199 1.0091 1.0050 -
DKTM 0.7455 0.9831 1.0025 -
RDKTM 1.0199 1.0093 1.0052 -
h=1 HDF-P3-7 β 1.0203 1.0070 1.0030 1.0010
ARS-T9 1.0343 1.0221 1.0103 1.0032
T3BL(R) 0.9204 0.9763 0.9938 0.9985
T3BL 0.9161 0.9754 0.9937 0.9985
DST-BK 1.0221 1.0124 1.0046 -
DST-BL 1.0279 1.0182 1.0124 -
DKTM 1.0465 1.0182 1.0093 -
RDKTM 1.0290 1.0160 1.0087 -

Figure 15: Normalized central deflection and bending moment of the circular plate. 

Furthermore, the shear force and bending moments along the radius of both thin (R/h=50) and thick (R/h=5) clamped plate are also compared with the analytical solutions (Batoz and Dhatt, 1990). These numerical results are calculated by 96 proposed elements. And all data are listed in Table 7 and plotted in Figure 16. The analytical solutions are given as

(60)

Table 7: Resultants along the radius of clamped circular plate obtained by the proposed element. 

Mr
r 0 0.625 1.25 1.875 2.5 3.125 3.75 4.375 5
R/h=50 Result 2.0265 1.9389 1.6944 1.2932 0.7199 −0.0066 −0.8991 −1.9727 −3.1388
R/h=5 Result 2.0266 1.9363 1.6935 1.2951 0.7237 −0.0028 −0.8908 −1.9499 −3.1398
Exact 2.0313 1.9507 1.7090 1.3062 0.7422 0.0171 −0.8691 −1.9165 −3.125
r 0 0.625 1.25 1.875 2.5 3.125 3.75 4.375 5
R/h=50 Result 2.0265 1.9806 1.8438 1.6143 1.2906 0.8760 0.3635 −0.2636 −0.9611
R/h=5 Result 2.0266 1.9818 1.8436 1.6098 1.2854 0.8731 0.3613 −0.2476 −0.8619
Exact 2.0313 1.9849 1.8457 1.6138 1.2891 0.8716 0.3613 −0.2417 −0.9375
Tr
r 0 0.625 1.25 1.875 2.5 3.125 3.75 4.375 5
R/h=50 Result −0.0105 −0.3180 −0.6202 −0.9235 −1.2091 −1.5137 −1.8237 −2.1107 −2.5265
R/h=5 Result −0.0123 −0.3104 −0.6110 −0.9096 −1.2131 −1.5172 −1.8242 −2.1357 −2.4301
  Exact 0 −0.3125 −0.625 −0.9375 −1.25 −1.5625 −1.875 −2.1875 −2.5

Figure 16: Resultants along the radius of clamped circular plate obtained by the proposed element. 

From the resultant force solutions for clamped circular plate, the following deflection equation can be obtained:

(61)

Then, the strain energy of a quarter of plate can be calculated by employing the virtual work principle, that is to say, the strain energy is equal to a half of the work of the distributed load.

(62)

In present model, the strain energy (the same as complementary energy for linear elastic material) is evaluated via the resultant force results given by equation (54):

(63)

where the index h denotes numerical solutions.

Thus, the energy error will be calculated by (U Uh )/U, and the energy norm error is the square root of (U Uh )/U. Following these notations, the energy norm error results can be obtained and are listed in Table 8. It can be simply concluded that the convergence rate is about 1 for both thin and moderately thick plates.

Table 8: The results of energy norm error for clamped circular plates. 

Number of Elements 6 24 96 384 Analytical Solution
h=0.1
Strain Energy 23113.9 29560.3 31402.9 31882.3 32046.0
Energy norm error 0.5279 0.2785 0.1417 0.0715 0.0000
h =1
Strain Energy 30.0292 37.5967 39.8493 40.4912 40.7235
Energy norm error 0.5125 0.2771 0.1465 0.0755 0.0000

High accuracy is achieved by the proposed element HDF-P3-7β for both central deflection and bending moment of the circular plate. And it is of particular interest that the resultants along the radius coincide perfectly with the analytical solutions.

4.1.5 Rotational Frame Dependence Test

In order to show that there is no rotational frame dependence, 30° and 60° skew plates are tested with the global frame rotated. Central deflection of the plate is calculated using 2×8×8 elements with the frame rotated counterclockwise in steps of 15°. Parameters of these plates are almost the same as the aforementioned ones. Taking the distributed load to be 1, deflection results are listed in Table 9.

Table 9: Central deflection of skew plates with frame rotated. 

60° skew plate
h/L=0.001 Angle/° 0 15 30 45 60 75 90
wc /109 0.786461 0.786461 0.786461 0.786461 0.786462 0.786462 0.786461
30° skew plate
h/L=0.001 Angle/° 0 15 30 45 60 75 90
wc /108 0.422175 0.422176 0.422175 0.422174 0.422173 0.422174 0.422175
h/L=0.01 Angle/° 0 15 30 45 60 75 90
wc /105 0.423212 0.423212 0.423212 0.423211 0.423210 0.423211 0.423212

With the numerical error considered, it can be concluded from the above table that the solutions calculated by HDF-P3-7β are independent of the global frame rotation, which is because that corresponding completeness of resultants is guaranteed.

4.2 Numerical Examples: Free Vibration Analyses

4.2.1 Free Vibration of Square Plates

As shown in Figure 17, two plates with four hard simply supported (SSSS) edges and four fixed (CCCC) edges are calculated to evaluate the performance of the proposed element. The span, thickness, density, Young’s modulus, and Poisson’s ratio of the plate are denoted by L, h, ρ, E, and μ, respectively, and are given by

L=10, h=0.05, 1, 𝒬=8000, E=2×1011, μ=0.3.

Figure 17: SSSS square plate and CCCC square plate. 

Both thin (h=0.05) and moderately thick (h=1) plates are studied and the non-dimensional frequencies are calculated by

(64)

where D is the aforementioned bending stiffness of Mindlin plate given in equation (11).

For the SSSS thin square plate, the free vibration frequencies can be derived theoretically and given by the following equation (Leissa, 1969).

(65)

where m and n are natural numbers.

The reference solutions for CCCC thin plates are proposed by Young (1950). And in other cases, the reference solutions can be found in the book of Abassian et al. (1987).

Normalization of the five lowest non-dimensional frequencies are listed in Table 10, Table 11, Table 12 and Table 13, where those results of some other elements are also given for comparison. The results of various elements using a mesh of 8×8×2 elements are plotted in Figure 18 and Figure 19.

Table 10: Normalized non-dimensional frequency parameters for square plates with SSSS BCs (L/h=200). 

Mesh Mode 1 2 3 4 5
4×4 HDF-P3-7 β 1.0030 1.0091 1.0240 1.0147 1.0398
NS+ES-FEM 1.0310 0.9996 1.0287 0.9327 0.8426
S4R 1.0265 1.1323 1.1323 1.1066 1.4858
S3R 1.5936 1.2401 2.5619 2.3850 2.7674
MIN3 1.1346 1.2358 1.4331 1.4608 1.7325
DSG3 1.2520 1.2548 1.6837 1.5094 1.8314
8×8 HDF-P3-7 β 1.0006 1.0016 1.0049 1.0026 1.0099
NS+ES-FEM 1.0105 1.0049 1.0121 1.0172 1.0013
S4R 1.0065 1.0289 1.0289 1.0263 1.0809
S3R 1.1121 1.0721 1.2797 1.1649 1.1865
MIN3 1.0309 1.0541 1.0888 1.1148 1.1351
DSG3 1.0652 1.0666 1.1706 1.1564 1.1773
16×16 HDF-P3-7 β 1.0001 1.0002 1.0010 1.0004 1.0019
NS+ES-FEM 1.0045 1.0023 1.0039 1.0063 1.0010
S4R 1.0016 1.0069 1.0069 1.0064 1.0185
S3R 1.0093 1.0147 1.0240 1.0400 1.0313
MIN3 1.0074 1.0132 1.0207 1.0292 1.0314
DSG3 1.0158 1.0178 1.0416 1.0537 1.0445
32×32 HDF-P3-7 β 1.0000 0.9999 1.0001 0.9999 1.0002
S4R 1.0003 1.0016 1.0016 1.0014 1.0043
S3R 1.0014 1.0027 1.0038 1.0058 1.0058
Exact Frequencies 4.4430 7.0250 7.0250 8.8860 9.9350

Table 11: Normalized non-dimensional frequency parameters for square plates with SSSS BCs (L/h=10). 

Mesh Mode 1 2 3 4 5
4×4 HDF-P3-7 β 1.0012 0.9991 1.0132 1.0003 0.9890
NS+ES-FEM 1.0143 0.9670 0.9714 0.8750 0.8055
S4R 1.0242 1.1148 1.1148 1.0838 1.3126
S3R 1.0972 1.1669 1.2767 1.2845 1.3543
MIN3 1.1242 1.2467 1.3931 1.4544 1.6928
DSG3 1.1435 1.2091 1.3993 1.3598 1.5324
8×8 HDF-P3-7 β 0.9999 1.0012 1.0040 1.0009 1.0019
NS+ES-FEM 1.0024 0.9947 0.9953 0.9883 0.9718
S4R 1.0054 1.0264 1.0264 1.0223 1.0655
S3R 1.0202 1.0393 1.0564 1.0728 1.0802
MIN3 1.0324 1.0707 1.0914 1.1248 1.1560
DSG3 1.0273 1.0489 1.0761 1.0930 1.1084
16×16 HDF-P3-7 β 0.9994 1.0010 1.0016 1.0009 1.0010
NS+ES-FEM 0.9994 0.9981 0.9987 0.9965 0.9930
S4R 1.0007 1.0069 1.0069 1.0060 1.0157
S3R 1.0043 1.0101 1.0142 1.0189 1.0195
MIN3 1.0137 1.0328 1.0370 1.0529 1.0637
DSG3 1.0056 1.0123 1.0176 1.0233 1.0256
32×32 HDF-P3-7 β 0.9992 1.0007 1.0009 1.0006 1.0004
S4R 0.9995 1.0022 1.0022 1.0019 1.0040
S3R 1.0005 1.0030 1.0040 1.0052 1.0050
Reference Frequencies 4.3700 6.7400 6.7400 8.3500 9.2200

Table 12: Normalized non-dimensional frequency parameters for square plates with CCCC BCs (L/h=200). 

Mesh Mode 1 2 3 4 5
4×4 HDF-P3-7 β 1.0291 1.0545 1.0625 1.0565 1.0346
NS+ES-FEM 1.0292 0.9473 0.9710 0.8244 0.7602
S4R 1.0759 1.3165 1.3165 1.2463 5.0766
S3R 3.2192 3.2547 4.6649 4.4407 4.8971
MIN3 1.2351 1.3729 1.5641 1.5718 1.8285
DSG3 1.4035 1.4907 1.7466 1.6586 1.8642
8×8 HDF-P3-7 β 1.0063 1.0122 1.0129 1.0132 1.0243
NS+ES-FEM 1.0080 1.0032 1.0114 1.0137 0.9970
S4R 1.0167 1.0519 1.0519 1.0446 1.1268
S3R 1.3460 1.3442 1.6246 1.4788 1.6046
MIN3 1.0578 1.0892 1.1309 1.1635 1.1881
DSG3 1.1195 1.1422 1.2333 1.2492 1.2665
16×16 HDF-P3-7 β 1.0013 1.0024 1.0026 1.0025 1.0046
NS+ES-FEM 1.0023 1.0007 1.0029 1.0042 1.0001
S4R 1.0039 1.0118 1.0118 1.0102 1.0270
S3R 1.0286 1.0394 1.0534 1.0875 1.0644
MIN3 1.0138 1.0213 1.0306 1.0421 1.0431
DSG3 1.0299 1.0359 1.0584 1.0807 1.0649
32×32 HDF-P3-7 β 1.0001 1.0003 1.0003 1.0001 1.0006
S4R 1.0008 1.0026 1.0026 1.0021 1.0060
S3R 1.0032 1.0049 1.0067 1.0098 1.0093
Reference Frequencies 5.9992 8.5680 8.5680 10.4053 11.4734

Table 13: Normalized non-dimensional frequency parameters for square plates with CCCC BCs (L/h=10). 

Mesh Mode 1 2 3 4 5
4×4 HDF-P3-7 β 1.0195 1.0193 1.0254 1.0164 0.9736
NS+ES-FEM 0.9932 0.9130 0.9322 0.8066 0.7540
S4R 1.0627 1.1951 1.1951 1.1460 1.3704
S3R 1.1812 1.2414 1.3680 1.3245 1.3783
MIN3 1.2246 1.3692 1.5329 1.5583 1.7870
DSG3 1.2040 1.2556 1.4067 1.3543 1.4909
8×8 HDF-P3-7 β 1.0040 1.0061 1.0071 1.0050 1.0043
NS+ES-FEM 0.9968 0.9857 0.9872 0.9748 0.9597
S4R 1.0133 1.0389 1.0389 1.0310 1.0803
S3R 1.0364 1.0534 1.0784 1.0877 1.0976
MIN3 1.0708 1.1180 1.1479 1.1802 1.2196
DSG3 1.0429 1.0611 1.0951 1.1039 1.1196
16×16 HDF-P3-7 β 1.0004 1.0017 1.0019 1.0013 1.0012
NS+ES-FEM 0.9977 0.9954 0.9961 0.9929 0.9899
S4R 1.0023 1.0089 1.0089 1.0072 1.0181
S3R 1.0080 1.0128 1.0187 1.0224 1.0235
MIN3 1.0382 1.0632 1.0696 1.0870 1.1007
DSG3 1.0090 1.0144 1.0219 1.0265 1.0281
32×32 HDF-P3-7 β 0.9992 1.0002 1.0002 1.0000 0.9997
S4R 0.9997 1.0019 1.0019 1.0014 1.0038
S3R 1.0011 1.0029 1.0043 1.0053 1.0052
Reference Frequencies 5.7100 7.8800 7.8800 9.3300 10.1300

Figure 18: Non-dimensional frequency parameters of SSSS square plates. 

Figure 19: Non-dimensional frequency parameters of CCCC square plates. 

Thin square plates (h=0.05) with various boundary conditions (as shown in Figure 20) are also analyzed by the proposed element using a 16×16 mesh. These problems are more challenging and two sets of reference solutions are considered, first of which is calculated by 500×500 S4R elements in the commercial FEM software Abaqus (2009), while the second is suggested by Leissa (1969). The non-dimensional frequencies are given by

(66)

Figure 20: Square plates with various boundary conditions. 

Normalized results of the lowest four non-dimensional frequencies (normalized by the first set of reference frequencies) are listed in Table 14 and compared with the results obtained by some other elements.

Table 14: Normalized non-dimensional frequency parameters for square plates with various BCs (L/h=200). 

BCs Mode 1 2 3 4
SSSF HDF-P3-7 β 1.0004 0.9971 1.0021 0.9994
S4R 1.0039 1.0044 1.0181 1.0119
S3R 1.0063 1.0204 1.0199 1.0421
DSG3 1.0079 1.0232 1.0190 1.0422
ES-DSG3 1.0003 1.0038 1.0061 1.0111
NS-DSG3 0.9973 1.0006 0.9964 1.0011
Reference 11.6801 27.7330 41.1799 59.0191
Reference 11.68 27.76 41.20 59.07
SFSF HDF-P3-7 β 1.0007 0.9987 0.9914 1.0025
S4R 1.0050 1.0032 1.0052 1.0198
S3R 1.0042 1.0071 1.0212 1.0171
DSG3 1.0039 1.0146 1.0274 1.0147
ES-DSG3 1.0013 1.0004 1.0063 1.0073
NS-DSG3 0.9992 0.9967 1.0036 0.9995
Reference 9.6297 16.1169 36.6738 38.9315
Reference 9.631 16.13 36.72 38.94
CCCF HDF-P3-7 β 1.0031 0.9994 1.0066 0.9957
S4R 1.0104 1.0069 1.0318 1.0202
S3R 1.0240 1.0644 1.0432 1.0831
DSG3 1.0158 1.0454 1.0290 1.0563
ES-DSG3 0.9995 1.0061 1.0054 1.0162
NS-DSG3 0.9873 0.9987 0.9791 1.0065
Reference 23.9064 39.9572 63.1736 76.6350
CFCF HDF-P3-7 β 1.0034 1.0022 0.9951 1.0069
S4R 1.0125 1.0082 1.0053 1.0340
S3R 1.0125 1.0383 1.0659 1.0314
DSG3 1.0084 1.0304 1.0541 1.0227
ES-DSG3 1.0006 1.0018 1.0092 1.0296
NS-DSG3 0.9873 0.9908 1.0027 1.0290
Reference 22.1576 26.3782 43.5274 61.1340
Reference 22.17 26.40 43.6 61.2

As shown in above numerical examples, the free frequencies of square plates obtained by the proposed element exhibit higher precision over other elements. And ideal solutions can be obtained by element HDF-P3-7β for relatively higher order frequencies with only coarse meshes.

4.2.2 Free Vibration of Parallelogram Plates

As shown in Figure 21, a cantilevered 60° parallelogram plate with a fixed edge and three free edges is investigated. The density, Young’s modulus, and Poisson’s ratio are the same as those in previous square plate example, and two kinds of thickness-span ratio, thin and moderately thick cases, are considered:

L=100, h=0.1, 20.

Figure 21: Cantilevered parallelogram plate and the typical mesh. 

The non-dimensional frequencies are calculated by the following equation:

(67)

Compared with the reference solutions given by Karunasena et al. (1996) using pb-2 Ritz method, as presented in Table 15 and Table 16, better performance can be produced by the proposed element than other elements like S4R, DSG3 and MIN3. And the proposed element HDF-P3-7β possesses nearly the same (even better) accuracy as NS+ES-FEM. A direct comparison is plotted in Figure 22.

Table 15: Normalized non-dimensional frequency parameters for cantilevered parallelogram plates (L/h=1000). 

Mesh Mode 1 2 3 4 5
4×4 HDF-P3-7 β 0.9942 0.9590 0.8982 0.9485 0.8633
NS+ES-FEM 1.0105 1.0599 0.8825 0.9405 0.8898
S4R 1.0013 1.0144 1.1157 1.1181 1.1607
S3R 1.0697 1.1515 1.1704 1.3600 1.3789
MIN3 1.0449 1.1970 1.2493 1.5553 1.3544
DSG3 1.0699 1.2987 1.1442 1.6714 1.2836
8×8 HDF-P3-7 β 0.9980 0.9916 0.9756 0.9876 0.9623
NS+ES-FEM 1.0008 1.0194 0.9666 0.9825 0.9449
S4R 0.9987 1.0036 1.0338 1.0253 1.0410
S3R 1.0188 1.0676 1.0676 1.1053 1.1200
MIN3 1.0158 1.0682 1.0816 1.1528 1.1102
DSG3 1.0361 1.1572 1.0851 1.3968 1.1608
16×16 HDF-P3-7 β 0.9980 0.9983 0.9948 0.9965 0.9904
NS+ES-FEM 0.9965 1.0028 0.9936 0.9907 0.9736
S4R 0.9977 1.0004 1.0080 1.0061 1.0099
S3R 1.0028 1.0106 1.0170 1.0162 1.0159
MIN3 1.0033 1.0210 1.0252 1.0417 1.0312
DSG3 1.0133 1.0670 1.0427 1.1457 1.0635
32×32 HDF-P3-7 β 0.9977 0.9994 0.9985 0.9989 0.9973
S4R 0.9975 0.9997 1.0014 1.0013 1.0020
S3R 0.9987 1.0013 1.0028 1.0022 1.0021
Reference 0.399 0.954 2.564 2.628 4.190

Table 16: Normalized non-dimensional frequency parameters for cantilevered parallelogram plates (L/h=5). 

Mesh Mode 1 2 3 4 5
4×4 HDF-P3-7 β 1.0016 0.9946 0.9496 0.9467 0.9148
NS+ES-FEM 0.9992 0.9868 0.8976 0.8383 0.8217
S4R 1.0050 1.0200 1.1055 1.0798 1.1200
S3R 1.0379 1.0364 1.1012 1.0882 1.0734
MIN3 1.0735 1.2464 1.3598 1.5329 1.4603
DSG3 1.0178 1.0950 1.0678 1.3536 1.1525
8×8 HDF-P3-7 β 1.0037 1.0045 0.9920 0.9889 0.9821
NS+ES-FEM 1.0000 0.9863 0.9690 0.9359 0.9458
S4R 1.0019 1.0056 1.0273 1.0212 1.0314
S3R 1.0127 1.0129 1.0310 1.0226 1.0272
MIN3 1.0459 1.1257 1.1894 1.1957 1.2242
DSG3 1.0111 1.0315 1.0194 1.1054 1.0705
16×16 HDF-P3-7 β 1.0019 1.0018 0.9988 0.9977 0.9960
NS+ES-FEM 1.0008 0.9929 0.9899 0.9801 0.9816
S4R 1.0003 1.0013 1.0068 1.0055 1.0079
S3R 1.0042 1.0038 1.0086 1.0061 1.0080
MIN3 1.0347 1.0862 1.1360 1.1095 1.1482
DSG3 1.0053 1.0062 0.9989 1.0282 1.0209
32×32 HDF-P3-7 β 1.0005 1.0005 0.9997 0.9996 0.9990
S4R 0.9995 1.0002 1.0015 1.0015 1.0019
S3R 1.0011 1.0010 1.0022 1.0017 1.0020
Reference 0.377 0.817 1.981 2.165 3.103

Figure 22: Non-dimensional frequency parameters of cantilevered parallelogram plates. 

4.2.3 Free Vibration of Circular Plates

The lowest eight frequencies of clamped circular plates with the mesh shown in Figure 23 are investigated. 393 proposed elements and 221 nodes are used to discretize the whole plate, while the compared solutions given by other elements employ even finer meshes. The frequencies given by S3R use the same mesh as the proposed element, but the solutions of MIN3 and NS+ES-FEM employ a mesh of 394 elements and 222 nodes, and the solutions of DSG3, NS-DSG3 and ES-DSG3 employ a much finer mesh of 848 elements and 460 nodes.

Figure 23: Clamped circular plate and the mesh. 

The material property parameters are remained unchanged, and the radius and thickness of the circular plates are:

R=5, h=0.1, 1.

Non-dimensional parameters of the frequencies are computed by the following equation:

(68)

Reference solutions of this benchmark are suggested by Leissa (1969) (for the thin plate) and Irie et al. (1980) (for the thick case).

It is noticeable that the proposed element produces the most accurate solutions while employing the coarsest mesh, which owes to the analytical trial function method used to formulate the element. And the solutions given by different elements are listed in Table 17 and Table 18 and plotted in Figure 24, from which a visualized comparison may be obtained.

Table 17: Normalized non-dimensional frequency parameters for circular plates with clamped BCs (2R/h=100). 

Mode 1 2 3 4 5 6 7 8
HDF-P3-7 β 1.0023 1.0021 1.0031 1.0013 1.0031 1.0045 1.0017 1.0028
MIN3 1.0188 1.0451 1.0463 1.0822 1.0832 1.0821 1.1342 1.1380
NS+ES-FEM 1.0053 0.9950 0.9954 0.9875 0.9879 0.9843 0.9734 0.9741
S3R 0.9861 1.0167 1.0268 1.0513 1.0633 1.0421 1.0930 1.1009
DSG3 1.0077 1.0184 1.0188 1.0318 1.0321 1.0356 1.0470 1.0485
ES-DSG3 1.0024 1.0066 1.0070 1.0121 1.0128 1.0150 1.0191 1.0208
NS-DSG3 1.0041 1.0095 1.0103 1.0167 1.0178 1.0208 1.0255 1.0275
Reference 10.2158 21.2600 21.2600 34.8800 34.8800 39.7710 51.0400 51.0400

Table 18: Normalized non-dimensional frequency parameters for circular plates with clamped BCs (2R/h=10). 

Mode 1 2 3 4 5 6 7 8
HDF-P3-7 β 1.0040 0.9996 1.0000 0.9934 0.9944 1.0023 0.9888 0.9896
MIN3 1.0788 1.1386 1.1397 1.2045 1.2049 1.2300 1.2812 1.2827
NS+ES-FEM 1.0006 0.9893 0.9895 0.9763 0.9765 0.9852 0.9582 0.9593
S3R 0.9847 1.0015 1.0090 1.0148 1.0236 1.0178 1.0326 1.0345
DSG3 1.0066 1.0095 1.0099 1.0142 1.0145 1.0257 1.0226 1.0235
ES-DSG3 1.0014 1.0002 1.0005 1.0007 1.0009 1.0101 1.0047 1.0055
NS-DSG3 1.0042 1.0048 1.0058 1.0079 1.0088 1.0192 1.0152 1.0165
Reference 9.2400 17.834 17.834 27.214 27.214 30.211 37.109 37.109

Figure 24: Non-dimensional frequency parameters of clamped circular plates. 

5 CONCLUSION

A new triangular hybrid displacement function element is formulated and evaluated in this paper. This element is based on the modified principle of complementary energy and the so-called displacement function. With the locking-free formulae of Timoshenko’s beam employed and the assumed resultants satisfying all governing equations, the proposed element HDF-P3-7β achieves high accuracy for both displacement components and resultants, and it is totally locking-free, which is because that it’s a stress-based element and the thin plate solutions are included in the selected resultant fields. Cooperating with the inertia matrix of 3-node triangular isoparametric element, the proposed element can be successfully applied in free vibration problems and give precise solutions for vibration frequency. This element can be also integrated into the standard framework of finite element program conveniently.

This hybrid displacement function element method was firstly introduced by Cen et al. (2014). In this paper, a more rigorous but complicated description of the hybrid displacement function element method from the viewpoint of the principles in mechanics is given, which gives the formulae of the equivalent nodal forces directly.

Numerical examples show that the proposed element performs well, especially for the resultants and free vibration frequencies. For the circular plate problem, resultants along the radius coincide perfectly with the corresponding analytical solutions, which is hardly satisfied by many classical displacement-based Mindlin plate elements. Last but not least, this element is independent of the frame rotation, which is because that corresponding completeness of resultants is guaranteed.

Although the patch tests are passed, the convergence may not be guaranteed using existed theory, because the proposed element is not displacement-based. A following work is the rigorous proof for the convergence of the hybrid displacement function element method, and is still being studied. Since the employment of the inertia matrices of displacement-based element has not been studied before, it’s also of remarkable importance.

Acknowledgments

The authors would like to acknowledge the financial supports of the National Natural Science Foundation of China (Project No. 11272181), the Specialized Research Fund for the Doctoral Program of Higher Education of China (Project No. 20120002110080), and the Tsinghua University Initiative Scientific Research Program (Project No. 2014z09099).

References

Abassian, F., Hawswell, D.J., Knowles, N.C. (1987). Free vibration benchmarks, Glasgow: Department of Trade and Industry, National Engineering Laboratory. [ Links ]

Ayad, R., Dhatt, G., Batoz, J.L. (1998). A new hybrid-mixed variational approach for Reissner-Mindlin plates. The MiSP model. International Journal for Numerical Methods in Engineering 42(7):1149-1179. [ Links ]

Ayad, R., Rigolot, A., Talbi, N. (2001). An improved three-node hybrid-mixed element for Mindlin/Reissner plates. International Journal for Numerical Methods in Engineering 51(8):919-942. [ Links ]

Bao, Y., Cen, S., Li, C.F. (2017). Distortion-resistant and locking-free eight-node elements effectively capturing the edge effects of Mindlin-Reissner plates. Engineering Computations 34(2): in press. [ Links ]

Bathe, K.J. and Dvorkin, E.N. (1985). A four-node plate bending element based on Mindlin-Reissner plate theory and a mixed interpolation. International Journal for Numerical Methods in Engineering 21(2):367-383. [ Links ]

Bathe, K.J. and Wilson, E.L. (1973). Solution methods for eigenvalues problems in Structural Mechanics. International Journal for Numerical Methods in Engineering 6(2):213-226. [ Links ]

Bathe, K.J., Cho, S.W., Buchalem, M.L., Brezzi, F. On our MITC plate bending/shell elements, in: Analytical and Computational Models for shells, Noor A.K. et al., Eds., (1989) CED 3:261-278, ASME Special Publication. [ Links ]

Batoz, J.L. and Dhatt, G. (1990). Modèlisation des Structures par Eléments Finis, Vol. 2: Poutres et Plaques, Editions Hermès, Paris. [ Links ]

Batoz, J.L. and Kaliti, I. (1992). On a simple triangular Reissner/Mindlin plate element based on incompatible modes and discrete constrains. International Journal for Numerical Methods in Engineering 35(8):1603-1632. [ Links ]

Batoz, J.L. and Lardeur, P. (1989). A disrete shear triangular nine dof element for the analysis of thick to very thin plates. International Journal for Numerical Methods in Engineering 28(3):533-560. [ Links ]

Batoz, J.L., Bathe, K.J., Ho, L.W. (1980). A study of three-node triangular plate bending elements. International Journal for Numerical Methods in Engineering 15(12):1771-1812. [ Links ]

Belytschko, T. and Bachrach, W.E. (1986). Efficient implementation of quadrilaterals with high coarse-mesh accuracy. Computer Methods in Applied Mechanics and Engineering 54(3):279-301. [ Links ]

Bletzinger, K.U., Bischoff, M., Ramm, E. (2000). A unified approach for shear-locking-free triangular and rectangular shell finite elements. Computer & Structure 75(3):321-334. [ Links ]

Cen, S. and Shang, Y. (2015). Developments of Mindlin-Reissner plate elements. Mathematical Problems in Engineering 2015:1-12. [ Links ]

Cen, S., Shang, Y., Li, C.F., Li, H.G. (2014). Hybrid displacement function element method: a simple hybrid-Trefftz stress element method for analysis of Mindlin-Reissner plate. International Journal for Numerical Methods in engineering 98(3):203-234. [ Links ]

Chen, W.J. and Cheung, Y.K. (2001). Refined 9-Dof triangular Mindlin plate elements. International Journal for Numerical Methods in Engineering 51(11):1259-1281. [ Links ]

Dassault Systèmes Simulia Corp., (2009). Abaqus 6.9 HTML Documentation, Providence, RI. [ Links ]

Hu, H.C. (1984). Variational Principle of Theory of Elasticity with Applications, Science press, Beijing. [ Links ]

Hughes, T.J.R. and Tezduyar, T.E. (1981). Finite elements based upon Mindlin plate theory with particular reference to the four-node bilinear isoparametric element. Journal of Applied Mechanics 48(3):587-596. [ Links ]

Hughes, T.J.R., Taylor, R.L., Kanoknukulchai, W. (1977). A simple and efficient finite element for plate bending. International Journal for Numerical Methods in Engineering 11(10):1529-1543. [ Links ]

Irie, T., Yamada, G., Aomura, S. (1980). Natural frequencies of Mindlin circular plates. Journal of Applied Mechanics 47(3):652-655. [ Links ]

Karunasena, W., Liew, K.M., Al-Bermani, F.G.A. (1996). Natural frequencies of thick arbitrary quadrilateral plates using the pb-2 Ritz method. Journal of Sound and Vibration 196(4):371-385. [ Links ]

Katili, I. (1993). A new discrete Kirchhoff-Mindlin element based on Mindlin-Reissner plate theory and assumed strain fields-Part I: an extended DKT element for thick-plate bending analysis. International Journal for Numerical Methods in Engineering 36(11):1859-1883. [ Links ]

Leissa, A.W. (1969). Vibration of plates, OHIO STATE UNIV COLUMBUS. [ Links ]

Liu, G. (1989). A method for large scale finite element dynamic analysis. Chinese Journal of Computational Mechanics 6(1):247-251 (in Chinese). [ Links ]

Liu, G.R., Nguyen-Thoi, T., Lam, K.Y. (2009b). An edge-based smoothed finite element method (ES-FEM) for static, free and forced vibration analyses of solids. Journal of Sound and Vibration 320(4-5):1100-1130. [ Links ]

Liu, G.R., Nguyen-Thoi, T., Nguyen-Xuan, H., Lam, K.Y. (2009a). A node-based smoothed finite element (NS-FEM) for upper bound solutions to solid mechanics problems. Computer & Structures 87(1-2):14-26. [ Links ]

Morley, L.S.D. (1963). Skew Plates and Structures, Macmillan, New York. [ Links ]

Nguyen-Xuan, H., Rabczuk, T., Bordas, S., Debongnie, J.F. (2008). A smoothed finite element method for plate analysis. Computer Methods in Applied Mechanics and Engineering 197(13-16):1184-1203. [ Links ]

Razzaque, A. (1973). Program for triangular bending elements with derivative smoothing. International Journal for Numerical Methods in Engineering 6(3):333-343. [ Links ]

Rezaiee-Pajand, M. and Karbon, M. (2014). Hybrid stress and analytical functions for analysis of thin plate bending. Latin American Journal of Solids and Structures 11(4):556-579. [ Links ]

Soh, A.K., Long, Z.F., Cen, S. (1999). A new nine DOF triangular element for analysis of thick and thin plates. Computational Mechanics 24(5):408-417. [ Links ]

Tabarrok, B. (1971). A variational principle for the dynamic analysis of continua by hybrid element method. International Journal of Solids and Structures 7(3):251-268. [ Links ]

Taylor, R.L. and Auricchio, F. (1993). Linked interpolation for Reissner-Mindlin plate elements: Part II-a simple triangle. International Journal for Numerical Methods in Engineering 36(18):3057-3066. [ Links ]

Tessler, A. and Hughes, T.J.R. (1985). A three-node Mindlin plate element with improved transverse shear. Computer Methods in Applied Mechanics and Engineering 50(1):71-101. [ Links ]

Wu, F., Liu, G.R., Cheng, A.G., He, Z.C. (2014). A new hybrid smoothed FEM for static and free vibration analyses of Reissner-Mindlin plates. Computational Mechanics 54(3):865-890. [ Links ]

Young, D. (1950). Vibration of rectangular plates by the Ritz method. Journal of Applied Mechanics-Transactions of the ASME 17(4):448-453. [ Links ]

Zienkiewicz, O.C., Taylor, R.L., Too, J.M. (1971). Reduced integration technique in general analysis of plates and shells. International Journal for Numerical Methods in Engineering 3(2):275-290. [ Links ]

Received: April 30, 2016; Revised: February 05, 2017; Accepted: March 08, 2017

*Corresponding Author, censong@tsinghua.edu.cn

Creative Commons License This is an open-access article distributed under the terms of the Creative Commons Attribution License