AP04-Bittnar1.vp 1 Introduction The basic unknowns are elongations in the directions of edges �21, �32, �13 and pure rotations �iz of a small neigh- bourhood of a node (considered rigid) about the z-axis pass- ing 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 Ma- trix”) supports the rigid body motion (RBM). 2 Element displacement Derivation of the element stiffness matrix relies on Hreni- koff ’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 Li (appendix A) are used to describe the basic displacement field uij(r, s) derived from pure deformations and rotations. Here, the pure deformations are defined through elongations of the individual edges and denoted as u21, u32, …, u13. They are composed of two parts – pure extensions due to node dis- placements 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 �uij(r, s) at the point r, s by � � �u r s N r s u L uij i ij i ij( , ) ( , )� � � � . (1) Pure extensions in terms of nodal rotations. Associating the ro- tation in node 1 with �1 � 1 we write the pure extension along the edges in terms of cubic shape functions in Fig. 2. k k r s L L l r s L L l � � 21 2 1 2 21 1 13 3 1 2 13 1 ( , ) , ( , ) , � � � � (2) where lij are the lengths associated with sides ij. Rotation �i causes in node (r, s) two displacements k�ij perpendicular to lines that are parallel with sides ij (Fig. 2). Superscript k means the node about which the rotation takes place. Extensions �iuij(r, s) of three lines parallel with edges, which pass point (r, s), due to the vertex rotations �i (super- script �i) need to be expressed. These elongations �iuij(r, s) are obtained by projecting k�ij in the corresponding directions ij. Elongations from the rotation in node 1 (Fig. 3) can be written as © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 75 Acta Polytechnica Vol. 44 No. 5 – 6/2004 Membrane Triangles with Drilling Degrees of Freedom P. Fajman An accurate triangular plane element with drilling degrees of freedom is shown in this paper. This element can be successfully used for solving linear and nonlinear problems. The main advantage of this element is that the stiffness matrix is obtained from pure deformations – elongations of the edges. This aproach is very suitable for nonlinear analysis, where the unbalanced forces can be obtained directly from elongations of edges. Keywords: triangular, element, elongations, rigid, displacement. 1 2 3 r,s � u21 � u21(r,s) 1 2 3 r,s N2(r,s) =L2 s r Fig. 1: Node displacement field �u21(r, s) 1 2 3 r,s 1�13(r,s) 1�21(r,s) 1 Fig. 2: Extension from pure rotation �1 � 1 � � 1 1 21 21 1 13 13 1 21 2 1 2 21 u r s l L L ij ij ij i ( , ) � � � � � � �n t n t n t j ijl L L� �1 13 3 1 2 13n t , (3) where i � 2, 3, 1 and j � 1, 2, 3, and { } { , } {cos , sin }n ij xij yij nij nijn n� � � � is the outward unit normal vector of the element side ij, { } { , } {cos , sin }tij xij yij nij nijt t� � � � 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 � i iju r s( , ) from three node rotations are obtained by sum- mation of the separate contributions �1u r sij ( , ), �2u r sij ( , ), �3u r sij ( , ). The total displacements along the element sides in node (r, s) due to the � uij (Eq. 1) and �uij (Eq. 3) may be expressed u r s L u l L L l L Lij i ij i jk k j i j jk k i j( , ) sin sin� � � � � � � � 2 2 , (4) or more compactly u Nu* *( , )r s � , (5) where vector u*( , ) { ( , ), ( , ), ( , )}r s u r s u r s u r s T� 21 32 13 repre- sents the pure displacements of node (r, s) in the directions of sides ij, u*( , ) { , , , , , }r s u u u T� � � �21 32 13 1 2 3� � � is the vector of six element conectors – elongations of the sides and pure node rotations, and N is the (3×6) matrix composed of the linear and cubic functions obtained directly from Eq. (4) for i � 2, 3, 1 and j � 1, 2, 3. 3 Pure element stiffness matrix The individual components of the side strains are ob- tained from the pure displacement vector u*( , )r s using the geometrical equations �ij ij ij ij ij j ij ij i u r s s u r s l L u r s l L � � � d d d d d d ( , ) ( , ) ( , ) . (6) Note that �ij corresponds to a relative deformation along the side ij. Substituting Eq. (5) into Eq. (6) gives � * *� Bu , (7) where �* { , , }� � � �21 32 13 T and B is the element strain matrix with the components B N l ik ik ij � � � (k � 1, …, 6). Components of the engineering strain � � { , , }� � �x y xy T are obtained applying standard transformation as � � � � �ij ij ij ij ij� �{cos , sin , sin cos } 2 2 �, and in matrix notation � � � � * *� � � ��A A 1 . (8) Substituting Eq. (7) into (8) finally provides � � �A Bu1 *. (9) 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 E V Vi s T V T T V T T ( ) * * * ,� � � � � � � � � 1 2 1 2 1 2 1 � � d du B A D A B u K B A D A�� 1BdV V , (10) where D is the elastic material stiffness matrix. Seven Gauss integration points with �( )h5 need to be used for numerical evalutions of K*, as the stiffness matrix con- tains functions of the fourth order. 4 Element stiffness matrix The pure element stiffness matrix is extended by includ- ing the rigid body motions through the static matrix T ob- tained from the relations between the pure displacements u* { , , , , , }� � � �u u u T21 32 13 1 2 3� � � and the global node displacements u � { , , , , , , , , }u v u v u v T1 1 1 2 2 2 3 3 3� � � . Rigid body motions of the edges are decomposed into translations and rotations. The relations between elongations of edges � uij and Car- tesian displacements ui, vi at the nodes are evident from Fig. 5 � � �u r s u v u u ij ij ij ij ij i j ij ( , ) cos sin ( ) cos ( � � � � � � � � � � � v vi j ij� ) sin .� (11) 76 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 44 No. 5 – 6/2004 1 2 3 1�21 � u23 � u13 Fig. 3: Relation between node rotation 1 and elongations 1 2 3 x y n21 n13 n32 t13 t21 t32 �1 �2 �3 Fig. 4: Introducing geometrical quantity The relations among pure vertex rotations �i and node rotations �i shown in Fig. 6, are � �i i� �� , (12) where �i is the pure rotation, �i represents the node rotation, � corresponds to a rigid body rotation given by � � � � � � � � � 1 2 � � � � u y v x . (13) The contribution of the pure rotations to the nodal dis- placements 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. u L u v L v L u y r s i i i r s i i i i i i i ( , ) , ( , ) , , ,� � � � � � � � 1 3 1 3 1 2 � � � � � � � � � � � � � � � � � � � � � 1 3 1 3 1 3 1 4 , , , � � L v x A c u b v i i i i i i i i i� � � � � � � � �1 3, . (14) Combining Eqs. (12), (13), (14), provides the relation between pure displacements u* and node displacements u in the form u Tu* � , (15) where u � { , , , , , , , , }u v u v u v T1 1 1 2 2 2 3 3 3� � � is the vector of nodal degrees of freedom (ui, vi are nodal Cartesian displace- ments and �i are nodal rotations). 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 Ei s( ) * *� � � 1 2 u T K Tu K T K TT T T . 4.1 Cantilever beam under tip load A standard test example often used in the literature cor- responds 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 � 0.25. The vertical deflection at point c according to the theory of elasticity, ref. [6], is wc � 0.3558 m. Table 1 lists the numerical results obtained for the tip de- flection. 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 cal- culated here with a mesh of 128 elements is particularly note- worthy (even outperforming Allman’s element). 77 Acta Polytechnica Vol. 44 No. 5 – 6/2004 1 2 3 u1 v1 u1 �u12 v1 �v12 u3 v3 x y u2 v2 Fig. 5: Rigid body motions – translation 1 2 3 �1 x y �� �1 Fig. 6 Rigid body motions - rotation 4.8 m 1 .2 x y F = 4 0 k N c 50 4 1 .6 6 4 1 .6 6 Reserve linear loading Fig. 7: Cantilever beam The three adopted mesh patterns are shown in Fig. 8. 4.2 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, � 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 rap- idly becomes worse for aspect ratios over 2:1. The results for the present element are compared in Fig. 12 with the re- sults 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 wc at point c are also given in Table 2. 78 Acta Polytechnica Vol. 44 No. 5 – 6/2004 Number of elements Deflection wc (m) CST LST Allman Lee Present 8 0.0909 0.2696 0.2068 0.240 32 0.1988 0.3550 0.3261 0.2958 0.3151 128 0.3115 0.3556 0.3471 0.3388 0.3480 Table 1 0.240 0.315 0.348 Fig. 8: Mesh patterns 250 Exact = 0.3558 0.296 0.347 0.348 w c 0.269 0.3 50 150 300100 200 Allman 0.35 0.25 2 4 0 7 2 0.315 Present 2 4 0.326 Lee 0.240 number of equations 0.339 Fig. 9: Dependence of the deflection wc on the number of equations 32 m 2 x y M=100 c Fig. 10: Cantilever beam Number of elements Deflection wc (m) CST EFFAND Allman 3i Allman 7 i Present 2×2 1.28 100.07 0.39 0.16 1.10 4×2 4.82 99.96 5.42 2.47 3.31 8×2 15.75 99.99 38.32 24.25 25.53 16×2 36.36 99.99 76.48 69.09 69.70 Table 2: Tip deflections for a beam under end moment 1.10 25.5 Fig. 11: Mesh patterns 2×2, 8×2 5 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 matrices – the pure stiffness matrix depending on elongations of the edges and the global matrix depending on nodal displace- ments ui, vi, �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. Appendix A Area Coordinates The use of area coordinates enables the relation for one edge or one vertex to be derived. The other relations are ob- tained by cyclic substitution i from 1 to 2, from 2 to 3, from 3 to 1. L A A a b x c y A i i i i i� � � �( ) 2 , where A is the area of the triangle, a x y x y b y y c x xi j k k j i j k i j k� � � � � �, , , and xi, yi are the coordinates of the node in Fig.13. Acknowledgment Financial support from GAČR 103 02 0658 and MSM 210000001 is gratefully acknowledged. References [1] Allman D. J: “A Compatible triangular Element includ- ing vertex rigid rotations for Plane elasticity, Analysis.” Computer & structures, Vol. 19 (1984). [2] Bergan P. G., Felippa C. A.: “A triangular membrane element with rotational degrees of freedom.” Computer methods in app. Mech. and Eng., Vol. 50 (1985), p. 25–69. [3] Felippa C. A., Alexander S.: “Membrane triangles with corner drilling freedoms – III. Implementation and per- formance evalution.” Finite Elements in Analysis and Design, Vol. 12 (1992), p. 203–239. [4] Ibrahimbegovic A., Taylor R. L., Wilson E. L.: “A robust quadrilateral membrane finite element with drilling de- grees of freedom.” Int. Journal for numerical meth. In. Eng., Vol. 30 (1990), p. 445–457. [5] Lee S. C., Yoo CH. H.: “A Novel Shell element including in plane Torque effect.” Computer & structures, Vol. 28 (1988). [6] Timoshenko S., Goodier J. N.: “Theory of elasticity.” 3rd Edn., McGraw-Hill, New York, 1970. Petr Fajman phone: +420 224 354 473 e-mail: fajman@fsv.cvut.cz Department of Structural Mechanics Czech Technical University Faculty of Civil Engineering Thákurova 7 166 29 Praha 6, Czech Republic 79 Acta Polytechnica Vol. 44 No. 5 – 6/2004 number of equations Exact 100� 3 6 3.3 7 2 w c 0.0 15050 100 Present EFFAN 100 50 87.1 2 8 8 69.7 1 8 25.5 1 4 4 Allman 200 250 CST 2×2 4×2 8×2 16×2 32×2mesh Fig. 12: Dependence of deflection wc and number of equations 1 2 3 x y x2 x1 x3 y1 y2 y3 A1 A2 A3 (x,y) Fig. 13: Area coordinates