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).
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).
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.
From the displacements given in equation (1), the strains paralleled with the xoy-plane at arbitrary point in Mindlin plate can be obtained:
where κ is the curvature vector. And from equations (2) and (3), the strain compatibility equations can be derived:
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:
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.
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
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:
where M, T, D _{ b } , and D _{ s } are the bending moment vector, shear force vector, bending elasticity matrix, and shear elasticity matrix, respectively;
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}):
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;
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:
where Σ_{ Se } denotes the common boundaries of subregions; the indexes 1 and 2 denote the two subregions that share a common boundary;
where _{ Sσ } 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:
where λ _{1} and λ _{2} are Lagrange multipliers.
Let the variation of equation (17) be zero:
Here the relations (6) between domain and boundary resultants have been already considered.
Considering the following constitutive equation:
where F is an operator matrix, the area integral in equation (18) can be simplified. Integrate by parts for the area integral, we obtain
in which is the adjoint operator of F; n is the corresponding direction matrix; and _{ ∂Te } denotes the boundary of the e-th subregion;
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:
Thus, equation (18) can be simplified as
In above equation, only curvilinear integral exists. From equation (7), and using the relation
equation (23) can be further simplified to be
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:
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
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.
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
Then, substitution of equation (29) into the constitutive equations (8) yields
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.
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.)
And the general solutions of F satisfy the following governing equation,
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}).
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
in which
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).
Firstly, the assumed resultant fields can be derived from the aforementioned solutions of displacement function F:
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:
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):
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
into equation (34) yields
Then, the rotation along normal direction is assumed to be linear function.
These three complicated equations can be rewritten in the form of matrix:
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:
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
then, the modified functional of complementary energy can be simplified as:
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:
in which q is the global nodal displacement vector of the plate. And above variational equations yield
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
Substitution of equation (49) into equation (48) yields
Above equation gives the discrete global equilibrium equation for this kind of element:
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;
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);
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
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
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.
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.
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.
The deflection and bending moment at the plate center
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.
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 |
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.
For the central deflection w _{C} and bending moment _{ My } at point C, ^{Razzaque (1973}) give out the finite difference solutions:
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.
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 |
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}).
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.
Results of central deflection and two principal bending moments are listed in Table 5, and plotted in Figure 13.
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 |
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.
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
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.
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 | - |
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
_{ 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 | |
_{ Mθ } | ||||||||||
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 |
From the resultant force solutions for clamped circular plate, the following deflection equation can be obtained:
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.
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):
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.
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.
60° skew plate | ||||||||
---|---|---|---|---|---|---|---|---|
h/L=0.001 | Angle/° | 0 | 15 | 30 | 45 | 60 | 75 | 90 |
_{ wc } /10^{9} | 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 } /10^{8} | 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 } /10^{5} | 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×10^{11}, μ=0.3.
Both thin (h=0.05) and moderately thick (h=1) plates are studied and the non-dimensional frequencies are calculated by
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}).
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.
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 |
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 |
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 |
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 |
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
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.
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.
The non-dimensional frequencies are calculated by the following equation:
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.
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 |
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 |
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.
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:
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.
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 |
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 |
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.