Two efficient hybrid-trefftz elements for plate bending analysis

ARTICLES

Two efficient hybrid-trefftz elements for plate bending analysis

Mohammad Rezaiee-Pajand^*; Mohammad Karkon

Civil Engineering Department, Ferdowsi University of Mashhad – Iran

ABSTRACT

This study is devoted to the analysis of the Reissner-Mindlin plate bending. In this paper, the hybrid-Trefftz strategy will be utilized. Two novel and efficient elements are formulated in details. They are a Triangular element (THT) and a quadrilateral element (QHT), which have 9 and 12 degrees of freedom, respectively. In this approach, two independent displacement fields are defined; one within the element and the other on the edges of the element. The internal field is selected in such a manner that the governing equation of thick plates could be satisfied. Boundary field is related to the nodal degree of freedoms by the boundary interpolation functions. To calculate these functions, the edges of the element are assumed to behave like a Timoshenko beam. The high accuracy and efficiency of the proposed elements and absence of the shear locking in these formulations are all proven, using various numerical tests.

Keywords: finite element, hybrid-Trefftz, thick plate, triangular element, quadrilateral element.

1 INTRODUCTION

Since the early days of the finite element method era, many elements have been proposed for the analysis of thick plates bending. Some of these elements have been formulated based on Rayleigh-Ritz method, i.e. the orthodox scheme. Routinely, polynomial functions have been used for displacement field in this technique, which do not satisfy the related governing differential equation. In this version of formulation, efficiency of the finite element is increased by enlarging the element's degrees of freedom [8, 28]. As a result, extending the degrees of freedom leads to more complex elements. It is obvious that not only the analysts would dislike this complexity, but it also adds to the costs of calculation. Between years 1965 to 1975, researchers' main concern was to achieve high accuracy by complex elements. On this base, numerous elements with various specifications have been developed to analyze thin and thick plates bending. Through these years, the ideas of the multi-field methods, such as the hybrid approach and mixed scheme, emerged and flourished considerably [11, 23].

The hybrid procedure has a vast ability to formulate efficient elements for analysis of the plate bending structures. Based on this methodology, the hybrid-Trefftz strategy (HT), were developed. Like the Ritz's process, the hybrid-Trefftz technique is suitable to solve solid mechanic problems. However, dissimilar to the Ritz's method, the assumed polynomial for the internal field establishes the governing equation of the structure. In other words, homogeneous and particular solution of the governing equation is used to set up the finite element formulation.

The hybrid-Trefftz method was first applied in examining distorted meshes by Jirousek and Leon in 1970s [15]. Because of high accuracy and efficiency, this technique soon caught the attention of other researchers as well. The mentioned approach has been successfully utilized in solving plane elasticity problems [9, 32], thin plates [14, 19], Reissner-Mindlin plates [7, 13, 22], shells [24], Helmholtz problem [33], heat conduction [10], contact problems [26], and three-dimensional structures [21]. In 1995, Jirousek and his colleagues presented a family of quadrilateral hybrid-Trefftz (HT) p-elements for moderately thick plates [16]. These researchers utilized the modified Bessel function of the first kind, and took advantage of the edge shear equilibrium for their derivation.

In hybrid-Trefftz method, weak form of compatibility between elements is achieved by applying independent boundary displacement field and using variational principles. The Main advantage of this formulation is that it allows for the integration on the boundary, which is not only easier than integration on the surface, but also makes it possible to create polygon elements. Advantages of this technique over conventional ways of the finite element procedures are as follows [14, 25]:

1. There is continuity between elements in the hybrid method, and there is no need to establish continuity as an additional condition.

2. It is possible to define numerous degrees of freedom on the boundary for the element, without involving further new and complex formulation.

3. Because of using polynomials, which establish the governing equation, it enjoys high levels of accuracy.

In this study, the edges of the suggested element are designed according to Timoshenko beam, and its deflection and torsion fields are determined. In fact, the boundary interpolation function is considered as a beam component. Interpolation functions for displacement and rotation fields are cubic and quadratic functions, respectively. The boundary field is depended on the nodal displacements by shape functions. By depicting these fields in the general coordinates, x-y, the interpolation functions could be acquired. With the use of this boundary field, two new triangular and quadrangular elements are proposed. Accuracy and efficiency of these formulations are proven by utilizing several numerical tests.

2 HYBRID-TREFFTZ FORMULATION

To analyze the plate bending structures, governing finite element equation is required. First, the hybrid-Trefftz functional will be defined. In hybrid-Trefftz functional, internal displacement of (u) and boundary displacement of () are employed independently. The following equation is the hybrid functional of total energy [7, 13]:

In this relation, W, b_i , T_i , and _i are strain energy density function, body force, tractions, and prescribed traction along Γ_T , respectively. The strain-energy density function W(u) is written by W(u) = 1/2σ_ij ε_ij. In functional (1), (u) is so selected that the governing equation would be satisfied. Besides, provides the continuity and the compatibility between the elements in the weak form. In order to achieve a useful formulation for finite elements, functional (1) must be stationary in terms of the independent fields, u and . The process of functional (1) stationary will lead to the following formulas:

It should be added that for utilizing these formulas, u_i , _i , and T_i must be specified. In addition, u_i and T_i could be written as the total of two homogeneous and particular parts. Boundary field of _i could be also associated with nodal displacements by interpolation functions:

In these formulas, {u_p} is the particular solution of the governing equation, and [Φ]{c} is its homogeneous solution. Moreover, {} is the nodal displacement and [] is boundary interpolation function. Also, {T_p} represents tractions resulted from {u_p} and [Θ]{c} represents tractions resulted from homogeneous solution of the governing equation. By inserting functions (4), (5), and (6) into formulas (2) and (3), the following equations will be available:

By finding {c} from the formula (7), and putting it in the equality (8), the succeeding result can be achieved:

In this equation, {h},{g}, [H] and [G] are defined in the below form:

Eventually, the stiffness matrix [K] and the nodal forces {f} are obtained, as follows:

These formulas are applicable to plane stress, bending plates, and shell problems. In the subsequent sections, the plates bending structures, based on the Reissner-Mindlin theory, will be analyzed.

3 THICK PLATE PROBLEM

One of the most widely-used theories in analyzing thick plates is Reissner-Mindlin hypothesis. This theory utilizes first order shear deformation for analysis. Furthermore, because of using independent displacement and rotation fields, it is only needed to establish compatibility of Cº among the elements. Finite element users are aware that establishment of such a connection is by far easier than fulfilling compatibility of C¹, which is necessary in Kirchhoff's theory. Equations of displacement, stress, strain, and the relationship between them, and also tractions in thick plate could be expressed as in the next formulas:

Tractions of {T}are displayed in Figure 1. In addition, matrices [L], [D], and [A] are obtained in the following manner:

In these equations, k_s is shear correction factor, which is usually, assumed 5/6. The direction cosines of the element's side, in the x and y directions, are n_x and n_y , respectively:

The static equilibrium equation for plates is given below:

The vector {P} in the present formula, is actually the vector of external loads entering the plate, and is defined as follows:

By inserting the formula (21) into the equation (20) and simplifying it, the governing equations of thick plate are obtained in the following appearance [22]:

Present formulas are based on the Reissner-Mindlin theory. The solution of the homogeneous part of the differential equation (22), which is independent of loading, could be found by solving the bi-harmonic one, ∇⁴w_c = 0. After solving the last equation in polar coordinates, the homogeneous solution is obtained in the below shape [12]:

It is suitable to utilize the following complex variables:

Substituting these equations into (23) leads to the below solution [25]:

Therefore, the values of w_j are accessible by applying the equations below. These formulas were first introduced by Herrera [12].

The first eleven terms of the expression can be found by using k = 0,1,2. As a result, the succeeding solutions are obtained:

Particular solution of the governing equation, without dependence on support conditions, is easily calculated for known loads. For the uniformly distributed load, a particular solution of the equation (22) is obtained in the below shape [22].

To find the homogeneous solution for the equation (3), Petrolito supposed θ_x and θ_y in the following manner [22]:

The unknown functions, f₁ and f₂ , are obtained by inserting the present relationships in the governing equation (3). Therefore, the homogeneous solution of equation (3) could be expressed as in below:

It is obvious that the internal rotation functions, θ_x and θ_y , are dependent on the deflection w. Consequently, the field function of particular solution is accessible in the following appearance:

In the particular and homogeneous solution, there is no advantage over x and y, and therefore, the element will be rotational invariant.

The minimum of terms which are selected from homogeneous solution depends on the element's degrees of freedom. Number of necessary terms (and not adequate) to prevent numerical instability and deficiency of the rank of the stiffness matrix is acquired through the next condition [14, 25]:

Here, n and r are numbers of the nodal degrees of freedom of the element under consideration and of the discarded rigid body motion terms, respectively. It should be mentioned that using modes of the rigid body motion, product spurious strain energy in the element. Therefore, these modes should be removed. Number of modes of the rigid body motion in the plate is three (r = 3). This includes a transitive mode and two rotational modes in the direction of the x-axis and y-axis. Taking these points in consideration, the formula (32) for the plate changes to the following form:

It was stated that use of the minimum number (m=n-r) of the Trefftz functions in Eq. (4) does not always guarantee a stiffness matrix with full rank. Moreover, the full rank may always be achieved by suitably augmenting m. The optimal value of m for a given type of element should be found by numerical experimentation [25].

4 BOUNDARY DISPLACEMENT FIELD

The boundary displacement functions are related to the nodal degrees of freedom by shape functions. In order to find boundary shape functions, the edges of the element are assumed to behave like a Timoshenko beam. This beam element was utilized by other researchers for plate bending analysis [6, 29, 30, 35]. The shape functions of the Timoshenko beam are calculated based on Figure 2. Linear field is also used for torsion of this beam. In order to compute shape function of the beam in Figure 2, a cubic displacement polynomial and a quadratic rotational field are selected. Moreover, it is assumed that shear strain has the constant value of γ₀ . Based on these assumptions, the following equations can be written:

In these relations, β₁, β₀ , α₀ and γ₀ are unknown parameters. In order to determine their values, the equation of shear strain for Timoshenko beam is first established. By utilizing the shear strain value equal to γ₀ , the subsequent equalities will be available:

In the present formula, coefficients of the terms s and s² are equivalent to zero. Therefore, in the succeeding lines, α₀, β₁ are determined in terms of the unknown parameter γ₀ :

At this stage, there is only one unknown constant γ₀ , which can be determined from the condition of minimum strain energy. It should be added that the structural strain energy is the sum of bending and shear strain energy. Bending strain energy is calculated in the following way:

In this equation, κ represents curvature and is determined below:

Substituting these equations into (37) leads to the below bending strain energy:

Besides, the next equation determines the energy of shear strain:

By adding the bending and shear strain energy together, total strain energy is found, as follows:

Implementing ∂U/∂γ₀ = 0 ∂γ₀ = 0 will give the following results:

By substitution β₁ ,β₀ ,α₀ and γ₀ into (34), the succeeding shape functions for Timoshenko beam can be found:

By assuming a linear interpolation for torsion field, the next results will be in hand:

By transferring these fields to the general coordinates, x-y, the shape functions for the element in Figure (3-2) can be written in the following form:

Figure 4 leads to the subsequent equations:

Consequently, the next relationship can be written:

Furthermore, the boundary field's

_x and _y have the following equations:

Boundary displacement field is determined by multiplying shape functions in the nodal displacements, as follows:

All entries of the shape function matrix are presented in the ^{appendix one}.

5 NUMERICAL TESTS

The performance of the proposed elements will be assessed in this section. To evaluate the accuracy, convergence rate, and robustness of the new formulation several plate problems will be solved. In order to demonstrate the properties of the proposed elements in details, the results of THT and QHT, will be compared to the following famous elements:

1. The triangular element of Kirchhoff's discrete plate (DKT) [3].

2. The triangular element of Mindlin's discrete plate (DKMT) [18].

3. The triangular element according to Reissner-Mindlin's theory (ARS-T9) [31].

4. The triangular 3-node element with 9 DOF (T3BL) [34].

5. The triangular element with reduced integration (T3BL(R)) [34].

6. The triangular element according to Mindlin's theory (RDKTM) [6].

7. The quadrilateral element of Kirchhoff's discrete plate (DKQ) [4].

8. The quadrilateral element of Mindlin's discrete plate (DKMQ) [17].

9. The quadrilateral element with reduced integration (Q4BL) [36].

10. Bathe-Dvorkin quadrilateral element (MITC4) [2].

11. The quadrilateral element for Reissner-Mindlin's plates (ARS-Q12) [30].

12. The quadrilateral 4-node element (AC-MQ4) [5].

5.1 Study of Trefftz function effect

In this part, the number of Trefftz functions will be examined carefully through various numerical examples in view of accuracy and efficiency of computation. For this purpose, a square plate with simply support is solved under a uniformly distributed load q. The obtained solutions confirm a typical pattern of accuracy when the number of function terms will increase. To demonstrate the effect of Trefftz function terms, the integer n represents the number of the used terms, in THT-n and QHT-n. Table 1 reveals the minimum number of Trefftz function terms, which are sufficient for present formulation, and compatible with the necessary condition of Eq. (33). The findings of this study show that the triangular element THT with 7 Trefftz functions and the quadrangular element QHT with 11 Trefftz functions have the highest level of accuracy for plate analysis.

5.2 Analysis of square plate

In this section, the solution for the square plate bending is conducted under a uniformly distributed load q, with simply and clamped supports. The thickness and side length of the plate are denoted by h and l, respectively; and the Poisson's ratio ν of the material is assumed to be 0.3. To reach a general conclusion, the results of the central displacements and moments from the very thin plate (h/l = 0.001) to the thick plate (h/l = 0.2) are presented in Tables 2 to 9. According to the numerical results, it is obvious that the proposed elements, THT and QHT, have a high precision and fast convergence behavior, whether it is used for thin or thick plate analysis. It is important to note that no shear locking happens in the thin plate limit. Furthermore, the suggested elements are insensitive to mesh distortion. The comparisons of the present results with those obtained by other authors are also presented in Tables 10 to 13.

5.3 Circular plate test

A circular plate subjected to a uniformly distributed load q, with two different boundary conditions, simply supported and clamped conditions, was analyzed with mesh subdivisions shown in Figure 6. The radius, modulus of elasticity and Poisson's ratio are, a, E, and ν, respectively. The results obtained from the analyses are presented in Tables 14 through 19. The accuracy and reliability of the proposed element answers are clearly illustrated by the findings. As it is seen, these results obtained using meshes with 3, 12 and 48 elements. Similar to the square plate case, the solutions obtained for THT and QHT elements, are close to exact answers when they are compared to other well-known elements.

5.4 Razzaque skew plate

Figure 7 shows a 60º skew plate which was originally studied by Razzaque [27]. This plate has simply supported on two opposite edges and free on the other two edges, and subjected to a uniformly distributed load. The results obtained using the THT and QHT elements, are shown in Tables 20 and 21. The referenced value was obtained by Razzaque with 16×16 finite difference mesh.

5.5 Morley skew plate

The simply supported skew plate illustrated in Figure (8) has been solved by Morley. This researcher used polar co-ordinates and a least-squares solution procedure to analyze an acute skew plate subject to a uniformly distributed load. It is worth emphasizing that this problem poses severe difficulties for numerical methods, since there is a singularity in the bending moments at the obtuse corner. Despite these difficulties, in the case of the thin-plate theory, Morley's solution is assumed to be near exact. Babuška and Scapolla also solved this problem as a 3-D elastic structure. The results were obtained using suggested elements, and compared to some other famous formulations, as shown in Tables (22) through (24) and Figure (9).

Table 25

6 CONCLUSION

In this study, the triangular element THT and the quadrangular element QHT are proposed for the analysis of the thick bending plates. Formulations of these elements are based on the hybrid-Trefftz method. The two displacement fields, one within the element and the other on the borders of the element, are supposed to be independent. The independent internal field is defined by utilizing the Trefftz scheme. By assuming that the element's edges are Timoshenko beam, the required interpolation functions and the boundary displacement field for this beam component are calculated. In order to examine the accuracy and efficiency of the proposed elements, several numerical tests were conducted. As it was proven by the numerical results, the answers of the suggested elements quickly converged in all the tests. The present formulations lead to good answers, even when the coarse meshes are utilized. It should be added that the accuracy and convergence rate of the element QHT is higher than the element THT.

Received 22 May 2011;

In revised form 23 Jul 2011

In the following, the entries of the interpolation matrix for the boundary displacement field are introduced:

• [1] I. Babuka and T. Scapolla. Benchmark computation and performance evaluation for a rhombic plate bending problem. Int. J. Num. Meth. Eng., 28:155179, 1989.
• [2] K.J. Bathe and E.H. Dvorkin. A four-node plate bending element based on Mindlin/Reissner plate theory and mixed interpolation. Int. J. Num. Meth. Eng., 21:367383, 1985.
• [3] J.L. Batoz, K.J. Bathe, and L.W. Ho. A study of three-node triangular plate bending element. Int. J. Num. Meth. Eng., 15:17711812, 1980.
• [4] J.L. Batoz and M.B. Tahar. Evaluation of a new quadrilateral thin plate bending element. Int. J. Numer. Meth. Eng., 18:16551677, 1982.
• [5] S. Cen, Y.Q. Long, Z.H. Yao, and S.P. Chiew. Application of the quadrilateral area co-ordinate method: A new element for MindlinReissner plate. Int. J. Numer. Meth. Eng., 66:145, 2006.
• [6] W. Chen and Y.K. Cheung. Refined 9-Dof triangular mindlin plate elements. Int. J.Num. Meth. Eng., 51:12591281, 2001.
• [7] Y.S. Choo, N. Choi, and B.C. Lee. A new hybrid-Trefftz triangular and quadrilateral plate elements. Applied Mathematical Modeling, 34:1423, 2010.
• [8] H.R. Dhananjaya, J. Nagabhushanam, P.C. Pandey, and Zamin Jumaat Mohd. New twelve node serendipity quadrilateral plate bending element based on mindlin-reissner theory using integrated force method. Structural Engineering and Mechanics, 36(5):625642, 2010.
• [9] M. Dhanasekar, J.J. Han, and Q.H. Qin. A hybrid-Trefftz element containing an elliptic hole. Finite Elements in Analysis and Design, 42:13141323, 2006.
• [10] Z.J. Fu, Q.H. Qin, and W. Chen. Hybrid-Trefftz finite element method for heat conduction in nonlinear functionally graded materials. Engineering Computations, 28(5), 2011.
• [11] A. Ghali and J. Chieslar. Hybrid finite elements. Journal of Structural Engineering, ASCE, 112:24782493, 1986.
• [12] I. Herrera. Boundary methods: an algebraic theory Pitman, London, 1984.
• [13] J. Jirousek. Hybid Trefftz plate bending elements with p-method capabilities. Int. J. Numer. Meth. Eng., 24:13671393, 1987.
• [14] J. Jirousek and L. Guex. The hybrid-Trefftz finite element model and its application to plate bending. Int. J. Numer. Meth. Eng., 23:651693, 1986.
• [15] J. Jirousek and J.N. Leon. A powerful finite element for plate bending. Comput. Struct., 12:7796, 1977.
• [16] J. Jirousek, A. Wrcjblewski, Q.H. Qin, and X.Q. He. A family of quadrilateral hybrid Trefftz p-element for thick plate analysis. Comput. Methods Appl. Mech. Engrg., 127:315344, 1995.
• [17] I. Katili. A new discrete kirchhoffmindlin element based on mindlinreissner plate theory and assumed shear strain fields part ii: An extended dkq element for thick-plate bending analysis. Int. J. Num. Meth. Eng., 36:18851908, 1993.
• [18] I. Katili. A new discrete kirchhoff-mindlin element based on mindlin-reissner plate theory and assumed shear strain fields - part i: An extended dkt element for thick-plate bending analysis. Int. J. Num. Meth. Eng., 36:18591883, 1993.
• [19] N. Leconte, B. Langrand, and M. Markiewicz. On some features of a plate hybrid-Trefftz displacement element containing a hole. Finite Elements in Analysis and Design, 46:819828, 2010.
• [20] Morley L.S.D. Skew plates and structures. Pergamon Press, Oxford, 1963.
• [21] K. Peters and W. Wagner. New boundary-type finite element for 2-D and 3-D elastic structures. Int. J. Numer. Methods Eng., 37:10091025, 1994.
• [22] J. Petrolito. Hybrid-Trefftz quadrilateral elements for thick plate analysis. Comput. Methods Appl. Mech. Eng., 78:331351, 1990.
• [23] T.H.H. Pian. State of the art development of hybrid/mixed finite element method. Finite Element in Analysis and Design, 21:520, 1995.
• [24] R. Piltner. On the systematic construction of trial functions for hybrid Trefftz shell elements. Computer Assisted Mechanics and Engineering Sciences, 4:633644, 1997.
• [25] Q.H. Qin. Trefftz finite element method and its applications. ASME, 58:316337, 2005.
• [26] Q.H. Qin and K.Y. Wang. Application of hybrid-Trefftz finite element method to frictional contact problems. Computer Assisted Mechanics and Engineering Sciences, 15:319336, 2008.
• [27] A. Razzaque. Program for triangular bending elements with derivative smoothing. Int. J. Numer. Meth. Eng., 6:333345, 1973.
• [28] M. Rezaiee-Pajand and M. Akhtary. A family of thirteennode plate bending triangular elements. Com. in Num. Meth. in Eng., 14:529537, 1998.
• [29] H. Sofuoglu and H. Gedikli. A refined 5-node plate bending element based on ReissnerMindlin theory. Commun. Numer. Meth. Engng., 23:385403, 2007.
• [30] A.K. Soh, S. Cen, Y.Q. Long, and Z.F. Long. A new twelve DOF quadrilateral element for analysis of thick and thin plates. Eur. J. Mech. A/Solids, 20:299326, 2001.
• [31] A.K. Soh, Z.F. Long, and S. Cen. A new nine DOF triangular element for analysis of thick and thin plates. Computational Mechanics, 24:408417, 1999.
• [32] C.O. Souza and S.P.B. Proen¸ca. A hybrid-Trefftz formulation for plane elasticity with selective enrichment of the approximations. Commun. Numer. Meth. Eng., 27(5):785804, 2011.
• [33] K. Sze and G. Liu. Hybrid-Trefftz six-node triangular finite element models for helmholtz problem. Computational Mechanics, 46:455470, 2010.
• [34] R.L. Taylor and F. Auricchio. Linked interpolation for reissner-mindlin plate element: part II a simple triangle. Int. J.Num. Meth. Eng., 36:30573066, 1993.
• [35] H. Zhang and J.S. Kuang. Eight-node ReissnerMindlin plate element based on boundary interpolation using Timoshenko beam function. Int. J.Num. Meth. Eng., 69:13451373, 2007.
• [36] O.C. Zienkiewicz, Z. Xu, F.Z. Ling, A. Samuelsson, and N.E. Wiberg. Linked interpolation for ReissnerMindlin plate element: part I a simple quadrilateral. Int. J. Num. Meth. Eng., 36:30433056, 1993.

appendix one

Publication Dates

• Publication in this collection
  20 Apr 2012
• Date of issue
  2012

History

• Received
  22 May 2011
• Accepted
  23 July 2011