Two efficient hybrid-Trefftz elements for plate bending analysis

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.


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.

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 Latin American Journal of Solids and Structures 9(2012) 43 -67 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 Ti 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, { d} 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: Latin American Journal of Solids and Structures 9(2012) 43 -67 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.

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 1 , 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: Latin American Journal of Solids and Structures 9(2012) 43 -67 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, ∇ 4 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 1 and f 2 , 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].

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 γ 0 .Based on these assumptions, the following equations can be written: In these relations, β 1 , β 0 , α 0 and γ 0 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 γ 0 , the subsequent equalities will be available: Latin American Journal of Solids and Structures 9(2012) 43 -67 In the present formula, coefficients of the terms s and s 2 are equivalent to zero.Therefore, in the succeeding lines, α 0 , β 1 are determined in terms of the unknown parameter γ 0 : At this stage, there is only one unknown constant γ 0 , 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: Latin American Journal of Solids and Structures 9(2012) 43 -67 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 β 1 , β 0 , α 0 and γ 0 into (34), the succeeding shape functions for Timoshenko beam can be found: Latin American Journal of Solids and Structures 9(2012) 43 -67 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: 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.

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: 8. The quadrilateral element of Mindlin's discrete plate (DKMQ) [17].

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.[34] 0.4062

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.[34] 0.1265  [34] 0.2310

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.[34] 0.6656   [34] 0.1848  [34] 0.8125

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.[27] 0.9589

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.Babuka 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).[20] 1.080

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.

APPENDIX: SHAPE FUNCTIONS FOR BOUNDARY DISPLACEMENT
In the following, the entries of the interpolation matrix for the boundary displacement field are introduced:

Figure 1
Figure 1 Boundary forces on an arbitrary direction.

Figure 4
Figure 4 Position of axes on the side.

Figure 6
Figure 6 Three mesh types for quadrant of a circular plate.

( a )
Displacement at centre (b) Bending moment Mmax at centre (c) Bending moment Mmin at centre