Membrane Triangles with Drilling Degrees of Freedom

The basic unknowns are elongations in the directions of edges 21, 32, 13 and pure rotations iz of a small neighbourhood of a node (considered rigid) about the z-axis passing through vertex i, where i 1, 2, 3. The effect of pure swing is transferred to the elongations of the edges and is taken into account in formulating of the element stiffness matrix. The static matrix (see section “Element Stiffness Matrix”) supports the rigid body motion (RBM).


Element displacement
Derivation of the element stiffness matrix relies on Hrenikoff 's idea that a triangular element can be represented by a system of lines parallel to the edges. Owing to this idea we further assume that displacements along individual element sides are mutually independent. Triangular coordinates L i (appendix A) are used to describe the basic displacement field u ij (r, s) derived from pure deformations and rotations. Here, the pure deformations are defined through elongations of the individual edges and denoted as u 21 , u 32 , …, u 13 . They are composed of two parts -pure extensions due to node displacements and pure extensions caused by node rotations.
Pure extensions in terms of nodal elongations. With the help of these quantities we may express the elongations of edges Pure extensions in terms of nodal rotations. Associating the rotation in node 1 with j 1 = 1 we write the pure extension along the edges in terms of cubic shape functions in Fig. 2.
where l ij are the lengths associated with sides ij.
Rotation j i causes in node (r, s) two displacements k D ij perpendicular to lines that are parallel with sides ij (Fig. 2). Superscript k means the node about which the rotation takes place. Extensions ji u ij (r, s) of three lines parallel with edges, which pass point (r, s), due to the vertex rotations j i (superscript j i ) need to be expressed. These elongations ji u ij (r, s) are obtained by projecting k D ij in the corresponding directions ij.
Elongations from the rotation in node 1 (Fig. 3) can be written as where i = 2, 3, 1 and j = is the unit vector of the element side ij (Fig. 4).
The other contributions from the rotations in nodes 2 and 3 are obtained by cyclic index replacement. Elongation The total displacements along the element sides in node (r, s) due to the D u ij (Eq. 1) and j u ij (Eq. 3) may be expressed or more compactly

Pure element stiffness matrix
The individual components of the side strains are obtained from the pure displacement vector u * ( , ) r s using the geometrical equations Note that e ij corresponds to a relative deformation along the side ij. Substituting Eq. (5) into Eq. (6) gives where e * { , , } = e e e The finite element formulation used in this paper is based on the principle of minimum potential energy. The element stiffness matrix is calculated by considering the strain energy in the form where D is the elastic material stiffness matrix.
Seven Gauss integration points with Â( ) h 5 need to be used for numerical evalutions of K * , as the stiffness matrix contains functions of the fourth order.

Element stiffness matrix
The pure element stiffness matrix is extended by including the rigid body motions through the static matrix T obtained from the relations between the pure displacements The relations among pure vertex rotations j i and node rotations Q i shown in Fig. 6, are where j i is the pure rotation, Q i represents the node rotation, Q corresponds to a rigid body rotation given by The contribution of the pure rotations to the nodal displacements u, v in Eq. (13) is neglected, since they cause no translations of the apex of a triangle. When using linear functions, the displacement field assumes the form.
Combining Eqs. (12), (13), (14), provides the relation between pure displacements u * and node displacements u in the form Note that the static matrix T is constant and does not cause any artifical strain or stress during the pure rigid body motions.
Relationship (15) is used in the expression for the strain energy (10), from which the element stiffness matrix K is obtained

Cantilever beam under tip load
A standard test example often used in the literature corresponds to a shear-loaded cantilever beam. In particular a cantilever beam of rectangular cross section is subjected to parabolic traction at its tip with the resultant F equal to 40 kN. The geometry and boundary conditions of the beam are shown in Fig. 7. The beam thickness is 1 m. The following material properties are used: Young's modulus E = 30000 kPa and Poisson's ratio n = 0.25. The vertical deflection at point c according to the theory of elasticity, ref. [6], is w c = 0.3558 m. Table 1 lists the numerical results obtained for the tip deflection. These are compared with the results obtained using the constant strain triangle (CST), the linear strain triangle (LST), the element of Allman [1] and Lee [5]. The values of the tip deflection obtained with the present element are closer to the exact solution than those obtained using the CST and the element due to Lee. The high accuracy of the results calculated here with a mesh of 128 elements is particularly noteworthy (even outperforming Allman's element).  The three adopted mesh patterns are shown in Fig. 8.

Cantilever beam under end moment
The next example is a cantilever beam subjected to the end moment M = 100. The specimen dimensions together with the support and loading conditions are illustrated in Fig. 10. The modulus of elasticity and the Poisson ratio are set equal to E = 768, n = 0.25 so that the exact tip deflection is 100. The adopted mesh patterns ranging from 32×2 to 2×2 (Fig. 11) are used to examine the element aspect ratio ranging from 1:1 to 16:1.
This element performs well for unit aspect ratios, but rapidly becomes worse for aspect ratios over 2:1. The results for the present element are compared in Fig. 12 with the results obtained using the constant strain triangle (CST), the element of Allman [1] and the EFFAND element due to Felippa and Alexander [2,3]. Numerical results for the tip deflexion w c at point c are also given in Table 2.

Conclusions
This paper presents the derivation of a simple triangular membrane finite element using the principle of minimum potential energy. The formulation of two stiffness matricesthe pure stiffness matrix depending on elongations of the edges and the global matrix depending on nodal displacements u i , v i , j i -enables nonlinear problems to be solved more effectively and quickly due to the simple calculation of the unbalanced forces. It is advantage in comparison with quadrilateral elements [4].
The element can represent exactly all constant strain states. The numerical results from the selected test example show that quite acceptable accuracy is achieved for practical applications.

Area Coordinates
The use of area coordinates enables the relation for one edge or one vertex to be derived. The other relations are obtained by cyclic substitution i from 1 to 2, from 2 to 3, from 3 to 1.