Applications of a Rotation-Free Triangular Element for Finite Strain Analysis of Thin Shells and Membranes Fernando G.Flores!and Eugenio Ofate2 1 Structures Department National University of Cordoba Casilla de Correo 916 5000 Cordoba-Argentina fflores@efn.unc.edu.ar Web page:http://www.efn.unc.edu.ar 2 International Centre for Numerical Methods in Engineering (CIMNE) Universidad Politecnica de Cataluna Campus Norte UPC,08034 Barcelona,Spain onate@cimne.upc.es Web page:http://www.cimne.upc.es Summary.This paper shows applications of a recently developed shell element to the analysis of thin shell and membrane structures.The element is a three node triangle with only translational DOFs (rotation free)that uses the configuration of the three adjacent elements to evaluate the strains.This allows to compute (constant) bending strains and (linear)membrane strains.A total Lagrangian formulation is used.Strains are defined in terms of the principal stretches.This allows to consider rubber materials and other type of materials using the Hencky stress-strain pair. An explicit central difference scheme is used to integrate the momentum equations. Several eramples,including inflation and deflation of membranes show the excellent convergence properties and robustness of the element for large strain analysis of thin shells and membranes. Key words:airbag inflation,deflation,shell triangular elements,rotation free, membranes 1 Introduction The simulation of the inflation of membrane structures is normally performed with membrane finite elements,i.e.no bending stiffness included.The formulation of such elements is simple as they only require Co continuity [1],in contrast with elements based on thin shell theory where Cl continuity implies important obstacles [2]in 69 E.Onate and B.Kroplin (eds.).Textile Composites and Inflatable Structures,69-88. 2005 Springer.Printed in the Netherlands
Applications of a Rotation-Free Triangular Element for Finite Strain Analysis of Thin Shells and Membranes Fernando G. Flores1 and Eugenio O˜nate2 1 Structures Department National University of C´ordoba Casilla de Correo 916 5000 Cordoba - Argentina ´ fflores@efn.unc.edu.ar Web page: http://www.efn.unc.edu.ar 2 International Centre for Numerical Methods in Engineering (CIMNE) Universidad Polit´ecnica de Catalu˜na Campus Norte UPC, 08034 Barcelona, Spain onate@cimne.upc.es Web page: http://www.cimne.upc.es Summary. This paper shows applications of a recently developed shell element to the analysis of thin shell and membrane structures. The element is a three node triangle with only translational DOFs (rotation free) that uses the configuration of the three adjacent elements to evaluate the strains. This allows to compute (constant) bending strains and (linear) membrane strains. A total Lagrangian formulation is used. Strains are defined in terms of the principal stretches. This allows to consider rubber materials and other type of materials using the Hencky stress-strain pair. An explicit central difference scheme is used to integrate the momentum equations. Several examples, including inflation and deflation of membranes show the excellent convergence properties and robustness of the element for large strain analysis of thin shells and membranes. Key words: airbag inflation, deflation, shell triangular elements, rotation free, membranes 1 Introduction The simulation of the inflation of membrane structures is normally performed with membrane finite elements, i.e. no bending stiffness included. The formulation of such elements is simple as they only require C0 continuity [1], in contrast with elements based on thin shell theory where C1 continuity implies important obstacles [2] in 69 E. Oñate and B. Kröplin (eds.), Textile Composites and Inflatable Structures, 69–88. © 2005 Springer. Printed in the Netherlands
70 Fernando Flores and Eugenio Onate the development of conforming elements.Triangular elements are naturally preferred as they can easily adapt to arbitrary geometries and due to the robustness of the associated mesh generators. When only the final configuration of the membrane is of interest implicit pro- grams are normally used,including special algorithms due to the lack of stiffness of the membrane when no tensile stresses are yet present.When the inflation/deflation process is of interest,the use of programs with explicit integration of the momentum equations are largely preferred.In the latter case linear strain triangles are normally not effective,specially when contact between surfaces is present.This implies a fine discretization of constant strain triangles to capture the details.what makes simula- tion quite expensive due to the time increment limitations.In this paper a triangular finite element with similar convergence properties to the linear strain triangle,but without its drawbacks,is used. Membrane structures components have some,although small,bending stiffness that in most of the cases is sensibly disregarded.However in many cases it may be convenient to include bending energy in the models due to the important regulariza- tion effect it supposes.Shell elements are of course more expensive due the increase in degrees of freedom(rotations)and integration points (through the thickness). In the last few years shell elements without rotation degrees of freedom have been developed (see 3]-10]among others),which make shell elements more efficient for both implicit and explicit integrators. The outline of this papers is as follows.Next section summarizes the rotation-free shell triangle used [10].Sec.3 shows convergence properties of the element in 2-d plane stress problems and 3-d linear bending/membrane problems.Sec.4 presents examples of inflation/deflation of membranes with and without bending stiffness. Finally Sec.5 summarizes some conclusions. 2 Formulation of the Rotation Free Shell Triangle The rotation-free EBST(for Enhanced Basic Shell Triangle)element has three nodes with three displacement degrees of freedom at each node.An element patch is defined by the central triangle and the three adjacent elements (Fig.1).This patch helps to define the membrane strains and curvature field within the central triangle (the EBST element)in terms of the displacement of the six patch nodes The node-ordering in the patch is the following (see Fig.1.a) The nodes in the main element (M)are numbered locally as 1,2 and 3.They are defined counter-clockwise around the positive normal The sides in the main element are numbered locally as 1,2,and 3.They are defined by the local node opposite to the side The adiacent elements(which are part of the cell)are numbered with the number associated to the common side The extra nodes of the cell are numbered locally as 4,5 and 6,corresponding to nodes on adjacent elements opposite to sides 1,2 and 3 respectively The connectivities in the adjacent elements are defined beginning with the extra node. Any convenient in plane (ti,t2)local cartesian coordinate system can be defined for the patch,with t3 the unit normal to the plane.The main features of the element formulation are the following:
70 Fernando Flores and Eugenio O˜nate the development of conforming elements. Triangular elements are naturally preferred as they can easily adapt to arbitrary geometries and due to the robustness of the associated mesh generators. When only the final configuration of the membrane is of interest implicit programs are normally used, including special algorithms due to the lack of stiffness of the membrane when no tensile stresses are yet present. When the inflation/deflation process is of interest, the use of programs with explicit integration of the momentum equations are largely preferred. In the latter case linear strain triangles are normally not effective, specially when contact between surfaces is present. This implies a fine discretization of constant strain triangles to capture the details, what makes simulation quite expensive due to the time increment limitations. In this paper a triangular finite element with similar convergence properties to the linear strain triangle, but without its drawbacks, is used. Membrane structures components have some, although small, bending stiffness that in most of the cases is sensibly disregarded. However in many cases it may be convenient to include bending energy in the models due to the important regularization effect it supposes. Shell elements are of course more expensive due the increase in degrees of freedom (rotations) and integration points (through the thickness). In the last few years shell elements without rotation degrees of freedom have been developed (see [3]–[10] among others), which make shell elements more efficient for both implicit and explicit integrators. The outline of this papers is as follows. Next section summarizes the rotation-free shell triangle used [10]. Sec. 3 shows convergence properties of the element in 2-d plane stress problems and 3-d linear bending/membrane problems. Sec. 4 presents examples of inflation/deflation of membranes with and without bending stiffness. Finally Sec. 5 summarizes some conclusions. 2 Formulation of the Rotation Free Shell Triangle The rotation-free EBST (for Enhanced Basic Shell Triangle) element has three nodes with three displacement degrees of freedom at each node. An element patch is defined by the central triangle and the three adjacent elements (Fig. 1). This patch helps to define the membrane strains and curvature field within the central triangle (the EBST element) in terms of the displacement of the six patch nodes. The node-ordering in the patch is the following (see Fig. 1.a) • The nodes in the main element (M) are numbered locally as 1, 2 and 3. They are defined counter-clockwise around the positive normal • The sides in the main element are numbered locally as 1, 2, and 3. They are defined by the local node opposite to the side • The adjacent elements (which are part of the cell) are numbered with the number associated to the common side • The extra nodes of the cell are numbered locally as 4, 5 and 6, corresponding to nodes on adjacent elements opposite to sides 1, 2 and 3 respectively • The connectivities in the adjacent elements are defined beginning with the extra node. Any convenient in plane (t1, t2) local cartesian coordinate system can be defined for the patch, with t3 the unit normal to the plane.The main features of the element formulation are the following:
Applications of a Rotation-Free Triangular Shell Element 71 4 1 5: 3 2 G? 2 M 3 3 5 6 6 (a) (b) Fig.1.Patch of elements for strain computation.(a)in spatial coordinates (b)in natural coordinates 1.The geometry of the patch formed by the central element and the three adjacent elements is guadratically interpolated from the position of the six nodes in the patch 2.The membrane strains are assumed to vary linearly within the central triangle and are expressed in terms of the(continuous)values of the deformation gradient at the mid side points of the triangle 3.An assumed constant curvature field within the central triangle is obtained using the values of the(continuous)deformation gradient at the mid side points. Details of the derivation of the EBST element are given below. 2.1 Definition of the Element Geometry and Computation of Membrane Strains As mentioned above a quadratic approximation of the geometry of the patch of four elements is chosen using the position of the six nodes.It is useful to define the patch in the isoparametric space using the nodal positions given in the Table 1 (see also Fig.1.b). Table 1.Isoparametric coordinates of the six nodes in the patch of Fig.1.b 123456 ξ0101-11 700111-1 The quadratic interpolation for the geometry is defined by (1)
Applications of a Rotation-Free Triangular Shell Element 71 Fig. 1. Patch of elements for strain computation. (a)in spatial coordinates (b)in natural coordinates 1. The geometry of the patch formed by the central element and the three adjacent elements is quadratically interpolated from the position of the six nodes in the patch 2. The membrane strains are assumed to vary linearly within the central triangle and are expressed in terms of the (continuous) values of the deformation gradient at the mid side points of the triangle 3. An assumed constant curvature field within the central triangle is obtained using the values of the (continuous) deformation gradient at the mid side points. Details of the derivation of the EBST element are given below. 2.1 Definition of the Element Geometry and Computation of Membrane Strains As mentioned above a quadratic approximation of the geometry of the patch of four elements is chosen using the position of the six nodes. It is useful to define the patch in the isoparametric space using the nodal positions given in the Table 1 (see also Fig. 1.b). Table 1. Isoparametric coordinates of the six nodes in the patch of Fig. 1.b 1 2 34 5 6 ξ 0101 -1 1 η 0011 1 -1 The quadratic interpolation for the geometry is defined by ϕ = *6 i=1 Niϕi (1) 1 2 5 3 4 6 η ξ . .. G G G 2 1 3 (b) 3 2 4 6 M 1 2 3 2 3 1 (a)
72 Fernando Flores and Eugenio Onate with i being the position vector of node i,s=1-g-n and N1=(+nN4=多(G-1) N2=ξ+n%N5=(E-1) (2) N3=I+SξN6=是(7-1) this interpolation allows to compute the displacement gradients at selected points in order to use an assumed strain approach.The computation of the gradients is performed at the mid side points of the central element (M)denoted by G1,G2 and G3 in Fig.1.b.This choice has the advantage that gradients at the three mid side points depend only on the nodes belonging to the two elements adjacent to each side.When gradients are computed at the common mid-side point of two ad- jacent elements,the same values are obtained,as the coordinates of the same four points are used.This in practice means that the gradients at the mid-side points are independent of the element where they are computed. The deformation gradient at the mid-side points of the element are obtained from the quadratic interpolations (1)as (a)G= Ni+3.ap+3,a=1,2,i=1,2,3(3) In Eg.(3)()'denotes values computed at the ith mid-side point. The cartesian derivatives of the shape functions are computed at the original configuration by the standard expression (4) where the Jacobian matrix at the original configuration is J= %tit taot (5) Once the deformation gradient is obtained,any convenient strain measure can be coupled.The membrane strains within the central triangle are now obtained using a linear assumed membrane strain field em,i.e. Em =Em (6) with 3 m=(1-2)eh+(1-2)e+(1-2n)ea.=∑N:e (7) where emm are the membrane strains computed at the three mid side points Gi (i=1,2,3 see Fig.2).In(7) N1=(1-2c),N2=(1-2ξ)andN3=(1-2m) (8) If,for example,Green-Lagrange strains are used, Emy =(pnpj-8yj) (9)
72 Fernando Flores and Eugenio O˜nate with ϕi being the position vector of node i, ζ = 1 − ξ − η and N1 = ζ + ξη N4 = ζ 2 (ζ − 1) N2 = ξ + ηζ N5 = ξ 2 (ξ − 1) N3 = η + ζξ N6 = η 2 (η − 1) (2) this interpolation allows to compute the displacement gradients at selected points in order to use an assumed strain approach. The computation of the gradients is performed at the mid side points of the central element (M) denoted by G1, G2 and G3 in Fig. 1.b. This choice has the advantage that gradients at the three mid side points depend only on the nodes belonging to the two elements adjacent to each side. When gradients are computed at the common mid-side point of two adjacent elements, the same values are obtained, as the coordinates of the same four points are used. This in practice means that the gradients at the mid-side points are independent of the element where they are computed. The deformation gradient at the mid-side points of the element are obtained from the quadratic interpolations (1) as (ϕα)Gi = ϕi α = *3 j=1 Ni j,αϕj + Ni i+3,αϕi+3 , α = 1, 2 , i = 1, 2, 3 (3) In Eq.(3) (·) i denotes values computed at the ith mid-side point. The cartesian derivatives of the shape functions are computed at the original configuration by the standard expression Ni,1 Ni,2 = J−1 Ni,ξ Ni,η (4) where the Jacobian matrix at the original configuration is J = ϕ0 ξ · t1 ϕ0 η · t1 ϕ0 ξ ξ · t2 ϕ0 η η · t2 (5) Once the deformation gradient is obtained, any convenient strain measure can be coupled. The membrane strains within the central triangle are now obtained using a linear assumed membrane strain field εˆm, i.e. εm = εˆm (6) with εˆm = (1 − 2ζ)ε1 m + (1 − 2ξ)ε2 m + (1 − 2η)ε3 m = *3 i=1 N¯iεi m (7) where εi m are the membrane strains computed at the three mid side points Gi (i = 1, 2, 3 see Fig. 2). In (7) N¯1 = (1 − 2ζ) , N¯2 = (1 − 2ξ) and N¯3 = (1 − 2η) (8) If, for example, Green-Lagrange strains are used, εmij = 1 2 (ϕi · ϕj − δij ) (9)
Applications of a Rotation-Free Triangular Shell Element 73 substituting(3)into(9)and using the usual membrane strain vector [11] Em emu,Emi2,Emia] (10) equation (7)gives pi1*p1-1 Em pi2pi2-1 (11) 2 i=1 2pi1·p2 The virtual membrane strains are expressed by pi1·p1 p2pi2 (12) 6pi1·pi2+pi1·p2, 2.2 Computation of Curvatures The curvatures(second fundamental form)of the middle surface are defined by[ Kag=5(pa·tgwg+pgt3a)=-tg·pag,a,B=1,2 (13) We will assume the following constant curvature field within each element Kia8 KaB (14) where RaB is the assumed constant curvature field obtained as 1 Ra8=- ta·PBo dA0 (15) and A is the area(in the original configuration)of the central element in the patch. Substituting (15)into(14)and integrating by parts the area integral gives the curvature vector [11]within the element in terms of the following closed line integral K11 K22 3·p1 dro (16) 2K12 t3·P2 -2-1 where ni are the components (in the local system)of the normals to the element sides in the initial configuration T. For the definition of the normal vector t3,the linear interpolation of the position vector over the central element is used. 3 =LM (17) i=1 where L are the standard linear shape functions of the central triangle(area co- ordinates)[11].In this case the tangent plane components are e-∑p:，a=1,2 (18) i=1
Applications of a Rotation-Free Triangular Shell Element 73 substituting (3) into (9) and using the usual membrane strain vector [11] εm = [εm11 , εm12 , εm12 ] T (10) equation (7) gives εm = *3 i=1 1 2 N¯i ' ϕi 1 · ϕi 1 − 1 ϕi 2 · ϕi 2 − 1 2ϕi 1 · ϕi 2 ( (11) The virtual membrane strains are expressed by δεm = *3 i=1 N¯i ' ϕi 1 · δϕi 1 ϕi 2 · δϕi 2 δϕi 1 · ϕi 2 + ϕi 1 · δϕi 2 ( . (12) 2.2 Computation of Curvatures The curvatures (second fundamental form) of the middle surface are defined by [] καβ = 1 2 (ϕα · t3β + ϕβ · t3α) = −t3 · ϕαβ , α, β = 1, 2 (13) We will assume the following constant curvature field within each element καβ = ˆκαβ (14) where ˆκαβ is the assumed constant curvature field obtained as κˆαβ = − 1 A0 M A0 M t3 · ϕβα dA0 (15) and A0 M is the area (in the original configuration) of the central element in the patch. Substituting (15) into (14) and integrating by parts the area integral gives the curvature vector [11] within the element in terms of the following closed line integral κ = ' κ11 κ22 2κ12 ( = 1 A0 M Γ0 M −n1 0 0 −n2 −n2 −n1 t3 · ϕ1 t3 · ϕ2 dΓ0 (16) where ni are the components (in the local system) of the normals to the element sides in the initial configuration Γ0 M. For the definition of the normal vector t3, the linear interpolation of the position vector over the central element is used. ϕM = *3 i=1 LM i ϕi (17) where LM i are the standard linear shape functions of the central triangle (area coordinates) [11]. In this case the tangent plane components are ϕMα = *3 i=1 LM i,αϕi , α = 1, 2 (18)
74 Fernando Flores and Eugenio Onate N x oN t3= loxx (19) From these expressions it is also possible to compute in the original configuration the element area A,the outer normals (n,n2)'at each side and the side lengths I.Equation (19)also allows to evaluate the thickness ratio A in the deformed configuration and the actual normal t3. Direction t3 can be seen as a reference direction.If a different direction than that given by (19)is chosen,at an angle 0 with the former,this has an influence of order 02 in the computation of the curvatures (see (23)below).This justifies (19) for the definition of ts as a function exclusively of the three nodes of the central triangle,instead of using the 6-node isoparametric interpolation. The numerical evaluation of the line integral in(16)results in a sum over the integration points at the element boundary which are,in fact,the same points used for evaluating the gradients when computing the membrane strains.As one integration point is used over each side,it is not necessary to distinguish between sides(i)and integration points(Gi). The explicit form of the gradient evaluated at each side Gi(3)from the quadratic interpolation is P2 Ni.2N2N.2N+3.2 (20) Pi+3 We note again that the gradient at each mid side point Gi depends only on the coordinates of the three nodes of the central triangle and on those of an additional node in the patch,associated to the side i where the gradient is computed. In this way the curvatures can be computed by LM 0 0 L (21) P2 An alternative form to express the curvatures,which is useful when their varia- tions are needed.is to define the vectors 3 hij (22) k=1 This gives K=h订·t3 (23) The variation of the curvatures can be obtained as 0 =空 0 LM LM +[a] (L出h+L) (LM22+L22) t3·dui) (24) (LMi2+L20i2)
74 Fernando Flores and Eugenio O˜nate t3 = ϕM 1 × ϕM 2 |ϕM 1 × ϕM 2 | = λ ϕM 1 × ϕM 2 . (19) From these expressions it is also possible to compute in the original configuration the element area A0 M , the outer normals (n1, n2) i at each side and the side lengths l M i . Equation (19) also allows to evaluate the thickness ratio λ in the deformed configuration and the actual normal t3. Direction t3 can be seen as a reference direction. If a different direction than that given by (19) is chosen, at an angle θ with the former, this has an influence of order θ2 in the computation of the curvatures (see (23) below). This justifies (19) for the definition of t3 as a function exclusively of the three nodes of the central triangle, instead of using the 6-node isoparametric interpolation. The numerical evaluation of the line integral in (16) results in a sum over the integration points at the element boundary which are, in fact, the same points used for evaluating the gradients when computing the membrane strains. As one integration point is used over each side, it is not necessary to distinguish between sides (i) and integration points (Gi). The explicit form of the gradient evaluated at each side Gi (3) from the quadratic interpolation is ϕi 1 ϕi 2 = Ni 1,1 Ni 2,1 Ni 3,1 Ni i+3,1 Ni 1,2 Ni 2,2 Ni 3,2 Ni i+3,2 ⎡ ⎢ ⎣ ϕ1 ϕ2 ϕ3 ϕi+3 ⎤ ⎥ ⎦ . (20) We note again that the gradient at each mid side point Gi depends only on the coordinates of the three nodes of the central triangle and on those of an additional node in the patch, associated to the side i where the gradient is computed. In this way the curvatures can be computed by κ = 2 *3 i=1 ⎡ ⎣ LM i,1 0 0 LM i,2 LM i,2 LM i,1 ⎤ ⎦ t3 · ϕi 1 t3 · ϕi 2 (21) An alternative form to express the curvatures, which is useful when their variations are needed, is to define the vectors hij = *3 k=1 LM k,iϕk j + LM k,jϕk i (22) This gives κij = hij · t3 (23) The variation of the curvatures can be obtained as δκ = 2 *3 i=1 ⎡ ⎣ LM i,1 0 0 LM i,2 LM i,2 LM i,1 ⎤ ⎦ '*3 i=1 Ni j,1(t3 · δuj ) Ni j,2(t3 · δuj ) + Ni i+3,1(t3 · δui+3) Ni i+3,2(t3 · δui+3) ( − *3 i=1 ⎡ ⎣ (LM i,11 11 + LM i,22 11) (LM i,11 22 + LM i,22 22) (LM i,11 12 + LM i,22 12) ⎤ ⎦(t3 · δui) (24)
Applications of a Rotation-Free Triangular Shell Element 75 where the projections of the vectors hi;over the contravariant base vectors have been included e%=hp0,a,i,j=1,2 (25) with p州=入p当xts (26) p9=-入pH×tg (27) In above expressions superindexes in L and u refer to element numbers whereas subscripts denote node numbers.As before the superindex M denotes values in the central triangle (Fig.1.a).Note that as expected the curvatures (and their variations)in the central element are a function of the nodal displacements of the six nodes in the four elements patch. Details of the derivation of (12)and(24)can be found in [10].The explicit expressions of the membrane and curvature matrices can be found in [12].The derivation of the element stiffnes matrix is described in [10,12].Also in [10,12] details of the quasi-static formulation and the fully explicit dynamic formulation are given. It must be noted that while the membrane strains are linear the curvature strains are constant.A full numerical integration of the stiffness matrix terms requires three points for the membrane part and one point for the bending part.Numerical experiments show that: when using one or three integration points the element is free of spurious energy modes and passes the patch test for initial curved surfaces the element with full (three point)integration leads to some membrane locking.This defect dissapears if one integration point is used for the membrane stiffness term. It can also be observed that: for large strain elastic or elastic-plastic problems membrane and bending parts can not be integrated separately,and a numerical integration trought the thick- ness must be performed for explicit integrators (hydro codes)is much more effective to use only one integration point for both the membrane and bending parts. Above arguments lead to reccomended the use of one integration point for both membrane and bending parts.This element is termed EBSTI to distinguish from the fully integrated one. 2.3 Boundary Conditions Elements at the domain boundary,where an adjacent element does not exist,de- serve a special attention.The treatment of essential boundary conditions associated to translational constraints is straightforward,as they are the degrees of freedom of the element.The conditions associated to the normal vector are crucial in this for- mulation for bending.For clamped sides or symmetry planes,the normal vector t3 must be kept fixed(clamped case),or constrained to move in the plane of symmetry (symmetry case).The former case can be seen as a special case of the latter,so we
Applications of a Rotation-Free Triangular Shell Element 75 where the projections of the vectors hij over the contravariant base vectors ϕ˜ Mα have been included α ij = hij · ϕ˜Mα , α, i, j = 1, 2 (25) with ϕ˜M1 = λ ϕM2 × t3 (26) ϕ˜M2 = −λ ϕM1 × t3 (27) In above expressions superindexes in LM j and δuk j refer to element numbers whereas subscripts denote node numbers. As before the superindex M denotes values in the central triangle (Fig. 1.a). Note that as expected the curvatures (and their variations) in the central element are a function of the nodal displacements of the six nodes in the four elements patch. Details of the derivation of (12) and (24) can be found in [10]. The explicit expressions of the membrane and curvature matrices can be found in [12]. The derivation of the element stiffnes matrix is described in [10, 12]. Also in [10, 12] details of the quasi-static formulation and the fully explicit dynamic formulation are given. It must be noted that while the membrane strains are linear the curvature strains are constant. A full numerical integration of the stiffness matrix terms requires three points for the membrane part and one point for the bending part. Numerical experiments show that: • when using one or three integration points the element is free of spurious energy modes and passes the patch test • for initial curved surfaces the element with full (three point) integration leads to some membrane locking. This defect dissapears if one integration point is used for the membrane stiffness term. It can also be observed that: • for large strain elastic or elastic-plastic problems membrane and bending parts can not be integrated separately, and a numerical integration trought the thickness must be performed • for explicit integrators (hydro codes) is much more effective to use only one integration point for both the membrane and bending parts. Above arguments lead to reccomended the use of one integration point for both membrane and bending parts. This element is termed EBST1 to distinguish from the fully integrated one. 2.3 Boundary Conditions Elements at the domain boundary, where an adjacent element does not exist, deserve a special attention. The treatment of essential boundary conditions associated to translational constraints is straightforward, as they are the degrees of freedom of the element. The conditions associated to the normal vector are crucial in this formulation for bending. For clamped sides or symmetry planes, the normal vector t3 must be kept fixed (clamped case), or constrained to move in the plane of symmetry (symmetry case). The former case can be seen as a special case of the latter, so we
76 Fernando Flores and Eugenio Onate symmetry 40 plane n= S= 0.s M original deformed Fig.2.Local cartesian system for the treatment of symmetry boundary conditions will consider symmetry planes only.This restriction can be imposed through the definition of the tangent plane at the boundary,including the normal to the plane of symmetrythat does not change during the process. The tangent plane at the boundary (mid-side point)is expressed in terms of two orthogonal unit vectors referred to a local-to-the-boundary Cartesian system (see Fig.2)defined as [p9,p] (28) where vector is fixed during the process while directionemerges from the intersection of the symmetry plane with the plane defined by the central element (M).The plane (gradient)defined by the central element in the selected original convective Cartesian system (t1,t2)is [p州，p] (29) the intersection line(side i)of this plane with the plane of symmetry can be written in terms of the position of the nodes that define the side (j and k)and the original length of the side i.e. 1 p,=应(pk-p) (30) That together with the outer normal to the side n=[m,n=[nt,nta (resolved in the selected original convective Cartesian system)leads to (31) where,noting that A is the determinant of the gradient,the normal component of the gradient can be approximated by in= (32)
76 Fernando Flores and Eugenio O˜nate Fig. 2. Local cartesian system for the treatment of symmetry boundary conditions will consider symmetry planes only. This restriction can be imposed through the definition of the tangent plane at the boundary, including the normal to the plane of symmetry ϕ0 n that does not change during the process. The tangent plane at the boundary (mid-side point) is expressed in terms of two orthogonal unit vectors referred to a local-to-the-boundary Cartesian system (see Fig. 2) defined as ϕ0 n, ϕ¯s (28) where vector ϕ0 n is fixed during the process while direction ϕ¯ s emerges from the intersection of the symmetry plane with the plane defined by the central element (M). The plane (gradient) defined by the central element in the selected original convective Cartesian system (t1, t2) is ϕM1 , ϕM2 (29) the intersection line (side i) of this plane with the plane of symmetry can be written in terms of the position of the nodes that define the side (j and k) and the original length of the side l M i , i.e. ϕi s = 1 l M i (ϕk − ϕj ) (30) That together with the outer normal to the side ni = [n1, n2] T = [n · t1, n · t2] T (resolved in the selected original convective Cartesian system) leads to ϕiT 1 ϕiT 2 = n1 −n2 n2 n1 ϕiT n ϕiT s (31) where, noting that λ is the determinant of the gradient, the normal component of the gradient ϕi n can be approximated by ϕi n = ϕ0 n λ|ϕi s| (32) M j k i j n=ϕ0 ’n t 1 t 2 symmetry plane original t 3 0i Y Z X i j ϕ’s ϕι =ϕ0 ’n i ϕΜ deformed ’1 ϕΜ ’2 t 3 i ’n
Applications of a Rotation-Free Triangular Shell Element 77 In this way the contribution of the gradient at side i to vectors has(22)results h L 0 i.1 0 21 =2 (33) 2h12 M 2 i,2 i.1 For a simple supported(hinged)side,the problem is not completely defined.The simplest choice is to neglect the contribution to the side rotations from the adjacent element missing in the patch in the evaluation of the curvatures via(16)[6,8.This is equivalent to assume that the gradient at the side is equal to the gradient in the central element.More precise changes can be however introduced to account for the different natural boundary conditions.One may assume that the curvature normal to the side is zero,and consider a contribution of the missing side to introduce this constraint.As the change of curvature parallel to the side is zero along the hinged side,both things lead to zero curvatures in both directions.For the case of a triangle with two sides associated to hinged sides.the normal curvatures to both sides must be set to zero. For a free edge the same approximation can be used but due to Poisson's effect this will lead to some error. For the membrane formulation of the EBST element,the gradient at the mid-side point of the boundary is assumed to be equal to the gradient of the main triangle. 2.4 Constitutive Models In the numerical experiments presented below two constitutive models have been used.A standard linear elastic orthotropic material and a hyper-elastic material for rubbers. For the case of rubbers,the Ogden [13 model extended to the compressible range is considered.The material behaviour is characterized by the strain energy density per unit undeformed volume defined as =P+∑ (34) 0=1 where K is the bulk modulus of the material,Ai are the principal strain ratios,J is the determinant of the deformation gradient F (J=A1A2A3),N,ui and a;are material parameters,ui,ai are real numbers such that ui;>0 (Vi =1,N)and N is a positive integer. 3 Convergence Studies In this section three examples are presented to show the convergence properties and the performance of present element.Examples are solved with a implicit program capable of dealing static/dynamic problems with moderate non-linearities
Applications of a Rotation-Free Triangular Shell Element 77 In this way the contribution of the gradient at side i to vectors hαβ (22) results in ⎡ ⎣ hT 11 hT 22 2hT 12 ⎤ ⎦ i = 2 ⎡ ⎣ LM i,1 0 0 LM i,2 LM i,2 LM i,1 ⎤ ⎦ ϕiT 1 ϕiT 2 = 2 ⎡ ⎣ LM i,1 0 0 LM i,2 LM i,2 LM i,1 ⎤ ⎦ n1 −n2 n2 n1 ϕiT n ϕiT s (33) For a simple supported (hinged) side, the problem is not completely defined. The simplest choice is to neglect the contribution to the side rotations from the adjacent element missing in the patch in the evaluation of the curvatures via (16) [6, 8]. This is equivalent to assume that the gradient at the side is equal to the gradient in the central element. More precise changes can be however introduced to account for the different natural boundary conditions. One may assume that the curvature normal to the side is zero, and consider a contribution of the missing side to introduce this constraint. As the change of curvature parallel to the side is zero along the hinged side, both things lead to zero curvatures in both directions. For the case of a triangle with two sides associated to hinged sides, the normal curvatures to both sides must be set to zero. For a free edge the same approximation can be used but due to Poisson’s effect this will lead to some error. For the membrane formulation of the EBST element, the gradient at the mid-side point of the boundary is assumed to be equal to the gradient of the main triangle. 2.4 Constitutive Models In the numerical experiments presented below two constitutive models have been used. A standard linear elastic orthotropic material and a hyper-elastic material for rubbers. For the case of rubbers, the Ogden [13] model extended to the compressible range is considered. The material behaviour is characterized by the strain energy density per unit undeformed volume defined as ψ = K 2 (ln J) 2 + *N p=1 µp αp J− αp 3 *3 i=1 λαp−1 i − 3 (34) where K is the bulk modulus of the material, λi are the principal strain ratios, J is the determinant of the deformation gradient F (J = λ1λ2λ3), N, µi and αi are material parameters, µi , αi are real numbers such that µiαi > 0 (∀i = 1, N) and N is a positive integer. 3 Convergence Studies In this section three examples are presented to show the convergence properties and the performance of present element. Examples are solved with a implicit program capable of dealing static/dynamic problems with moderate non-linearities
78 Fernando Flores and Eugenio Onate 3.1 Cook's Membrane Problem One of the main targets of the rotation-free triangular element is to obtain a mem- brane element with a behaviour similar to the linear strain triangle (LST).Such capacity is studied in this example [14]corresponding to a problem with an impor- tant amount of shear energy involved.This problem is also intended also to assess the ability of the element to distort.Figure 3.a shows the geometry of a tapered panel clamped on one side and with a uniformly distributed shear load in the oppo- site side.Figure 3.b presents the vertical displacement of point C(mid point of the loaded side)for the uniformly refined meshes considered as a function of the total number of nodes in the mesh. 25 -- 9 20 O 1 15 CST LST E=1 -D--EBST v=0.33 --4-=EBST1 10 (a) (b) 48 ● 20 40 60 80 100 Number of nodes Fig.3.Cook's membrane problem.(a)Geometry and load (b)Vertical displacement of point C for different meshes For the EBST element with three integration points,it can be seen that for the coarsest mesh (two linear elements),the measured displacement is slightly superior than the constant strain triangle(CST);but when the mesh is refined,values rapidly catch up with those obtained with the linear strain triangle.The element with only one integration point (EBST1)shows excellent predictions for coarse meshes and fast convergence properties for the reported displacement. 3.2 Cylindrical Roof In this example an effective membrane interpolation is of primary importance.Hence this is good test to assess the rotation-free element.The geometry is a cylindrical roof supported by a rigid diaphragm at both ends and it is loaded by a uniform dead weight (see Fig.4.a).Only one quarter of the structure is meshed due to symmetry conditions.Unstructured and structured meshes are considered.In the latter case two orientations are possible (Fig.4.a shows orientation B). Figure 4.b shows the normalized vertical displacement of the midpoint of the free side(point-C)over both (A and B)structured and (U)unstructured meshes as
78 Fernando Flores and Eugenio O˜nate 3.1 Cook’s Membrane Problem One of the main targets of the rotation-free triangular element is to obtain a membrane element with a behaviour similar to the linear strain triangle (LST). Such capacity is studied in this example [14] corresponding to a problem with an important amount of shear energy involved. This problem is also intended also to assess the ability of the element to distort. Figure 3.a shows the geometry of a tapered panel clamped on one side and with a uniformly distributed shear load in the opposite side. Figure 3.b presents the vertical displacement of point C (mid point of the loaded side) for the uniformly refined meshes considered as a function of the total number of nodes in the mesh. Fig. 3. Cook’s membrane problem. (a) Geometry and load (b) Vertical displacement of point C for different meshes For the EBST element with three integration points, it can be seen that for the coarsest mesh (two linear elements), the measured displacement is slightly superior than the constant strain triangle (CST); but when the mesh is refined, values rapidly catch up with those obtained with the linear strain triangle. The element with only one integration point (EBST1) shows excellent predictions for coarse meshes and fast convergence properties for the reported displacement. 3.2 Cylindrical Roof In this example an effective membrane interpolation is of primary importance. Hence this is good test to assess the rotation-free element. The geometry is a cylindrical roof supported by a rigid diaphragm at both ends and it is loaded by a uniform dead weight (see Fig. 4.a). Only one quarter of the structure is meshed due to symmetry conditions. Unstructured and structured meshes are considered. In the latter case two orientations are possible (Fig. 4.a shows orientation B). Figure 4.b shows the normalized vertical displacement of the midpoint of the free side (point-C) over both (A and B) structured and (U) unstructured meshes as F 48 44 16 E = 1 ν = 0.33 ( ) a C (b) Number of nodes p C Dislacement C Di l 0 20 40 60 80 100 5 10 15 20 25 CST LST EBST EBST1
