Finite Element Analysis of Membrane Structures Robert L.Taylor*1,Eugenio Onate2 and Pere-Andreu Ubach2 1 Department of Civil and Environmental Engineering University of California at Berkeley,USA rltOce.vulture.berkeley.edu 2 International Center for Numerical Methods in Engineering Edificio C1,Gran Capitan s/n 08034 Barcelona-Spain onate@cimne.upc.es Summary.This paper summarizes the development for a large displacement for- mulation of a membrance composed of three-node triangular elements.A formulation in terms of the deformation gradient is first constructed in terms of nodal variables. In particular,the use of the right Cauchy-Green deformation tensor is shown to lead to a particulary simple representation in terms of nodal quantities.This may then be used to construct general models for use in static and transient analyses. Key words:Membrane structures,finite element 1 Introduction The behavior of curved,thin bodies can be modeled by a membrane theory of shells. In such a theory only the in-plane stress resultants are included.The deformation state for a membrane may be represented by the position of points on the two- dimensional surface.General theories for shells may be specialized to those for a membrane by ignoring the resultant couples and associated changes in curvature deformations as well as any transverse shearing effects.A numerical approximation to the shell may then be constructed using a finite element approach.Examples for general shell theory and finite element solution may be found for small deforma- tions in standard books.]Theory for large deformation can proceed following the presentations of Simo et al.23.4.5]or Ramm et al.6.7.8.9.10.11,12.13 For the simplest shape finite element composed of a 3-node triangular form with displacement parameters at each vertex (a 9-degree of freedom element)it is far simpler to formulate the membrane behavior directly.This is especially true for large displacement response.Here the initially flat form of a simple triangular response is maintained throughout all deformation states.Consequently,one may proceed Visiting Professor,CIMNE,UPC,Barcelona,Spain. 47 E.Onate and B.Kroplin (eds.).Textile Composites and Inflatable Structures,47-68. 2005 Springer.Printed in the Netherlands
Finite Element Analysis of Membrane Structures Robert L. Taylor1, Eugenio O˜nate2 and Pere-Andreu Ubach2 1 Department of Civil and Environmental Engineering University of California at Berkeley, USA rlt@ce.vulture.berkeley.edu 2 International Center for Numerical Methods in Engineering Edificio C1, Gran Capit´an s/n 08034 Barcelona - Spain onate@cimne.upc.es Summary. This paper summarizes the development for a large displacement formulation of a membrance composed of three-node triangular elements. A formulation in terms of the deformation gradient is first constructed in terms of nodal variables. In particular, the use of the right Cauchy-Green deformation tensor is shown to lead to a particulary simple representation in terms of nodal quantities. This may then be used to construct general models for use in static and transient analyses. Key words: Membrane structures, finite element 1 Introduction The behavior of curved, thin bodies can be modeled by a membrane theory of shells. In such a theory only the in-plane stress resultants are included. The deformation state for a membrane may be represented by the position of points on the twodimensional surface. General theories for shells may be specialized to those for a membrane by ignoring the resultant couples and associated changes in curvature deformations as well as any transverse shearing effects. A numerical approximation to the shell may then be constructed using a finite element approach. Examples for general shell theory and finite element solution may be found for small deformations in standard books.[1] Theory for large deformation can proceed following the presentations of Simo et al. [2, 3, 4, 5] or Ramm et al. [6, 7, 8, 9, 10, 11, 12, 13] For the simplest shape finite element composed of a 3-node triangular form with displacement parameters at each vertex (a 9-degree of freedom element) it is far simpler to formulate the membrane behavior directly. This is especially true for large displacement response. Here the initially flat form of a simple triangular response is maintained throughout all deformation states. Consequently, one may proceed Visiting Professor, CIMNE, UPC, Barcelona, Spain. 47 E. Oñate and B. Kröplin (eds.), Textile Composites and Inflatable Structures, 47–68. © 2005 Springer. Printed in the Netherlands
48 Robert L.Taylor,Eugenio Onate and Pere-Andreu Ubach directly with the construction of the kinematic behavior,even in the presence of large strains.This approach is followed in the present work. The loading of a membrane is often by pressures which remain normal to the surface throughout all deformations.Such follower loading generally leads to a form which yields an unsymmetric tangent matrix.Such formulation has been presented in works by Schweizerhof and Ramm4l and by Simo et al The general approach presented in the last work is used for the special case of the flat triangular element used in this work. The formulation included in the present study includes inertial and damping terms based on second and first time derivatives of the motion.These are discretized in time using standard techniques (e.g.,the Newmark method6,171).Both explicit and implicit schemes are presented together with all linearization steps needed to implement a full Newton type solution.The inclusion of the damping term permits solution of the first order form in order to obtain a final static solution.Generally, the first order form is used until the final state is reached at which point the rate terms are deleted and the full static solution achieved using a standard Newton iterative method. The work presented is implemented in the general purpose finite element solution system FEAPlsl and used to solve example problems.The solution to some basic example problems are included to show the behavior of the element and solution strategies developed. 2 Governing Equations Reference configuration coordinates in the global Cartesian frame are indicated in upper case by X and current configuration in lower case by c.The difference between the coordinates defines a displacement u. Using standard interpolation for a linear triangle positions in the element may be specified as X=Eax (1) in the reference configuration and =ata (2) for the current configurations.If necessary,the displacement vector may be deduced as u =faua (3) In the above x,ua denote nodal values of the reference coordinates,current coordinates and displacement vector,respectively.Furthermore,the natural (area) coordinates satisfy the constraint 51+51+53=1 (4) It is convenient to introduce a surface coordinate system denoted by Yi,Y2 with normal direction N in the reference state and yi,y2 with normal direction n in the current state (see Fig.1)
48 Robert L. Taylor, Eugenio O˜nate and Pere-Andreu Ubach ˜ directly with the construction of the kinematic behavior, even in the presence of large strains. This approach is followed in the present work. The loading of a membrane is often by pressures which remain normal to the surface throughout all deformations. Such follower loading generally leads to a form which yields an unsymmetric tangent matrix. Such formulation has been presented in works by Schweizerhof and Ramm[14] and by Simo et al. [15] The general approach presented in the last work is used for the special case of the flat triangular element used in this work. The formulation included in the present study includes inertial and damping terms based on second and first time derivatives of the motion. These are discretized in time using standard techniques (e.g., the Newmark method[16, 17] ). Both explicit and implicit schemes are presented together with all linearization steps needed to implement a full Newton type solution. The inclusion of the damping term permits solution of the first order form in order to obtain a final static solution. Generally, the first order form is used until the final state is reached at which point the rate terms are deleted and the full static solution achieved using a standard Newton iterative method. The work presented is implemented in the general purpose finite element solution system FEAP[18] and used to solve example problems. The solution to some basic example problems are included to show the behavior of the element and solution strategies developed. 2 Governing Equations Reference configuration coordinates in the global Cartesian frame are indicated in upper case by X and current configuration in lower case by x. The difference between the coordinates defines a displacement u. Using standard interpolation for a linear triangle positions in the element may be specified as X = ξα X˜ α (1) in the reference configuration and x = ξα x˜α (2) for the current configurations. If necessary, the displacement vector may be deduced as u = ξα u˜ α (3) In the above X˜ α , x˜α, u˜ α denote nodal values of the reference coordinates, current coordinates and displacement vector, respectively. Furthermore, the natural (area) coordinates satisfy the constraint ξ1 + ξ1 + ξ3 = 1 (4) It is convenient to introduce a surface coordinate system denoted by Y1, Y2 with normal direction N in the reference state and y1, y2 with normal direction n in the current state (see Fig. 1)
Finite Element Analysis of Membrane Structures 49 t y 1 2 Fig.1.Description of coordinates for triangular element Placing the origin of the Yi and yi coordinates at nodal location and respectively,the unit base vectors may be constructed from the linear displacement triangle by aligning the first vector along the 1-2 side.Accordingly,we define the first unit vector as 元2-龙1 4221 v1=2--A2可 (5) where and lall -(a'a)"3. Next a vector normal to the triangle is constructed as v3=△元21×△231 (6) and normalized to a unit vector as U3 n=Tosll () The vector v3 plays a special role in later development of nodal forces for follower pressure loading as it is twice the area of the triangle times the unit normal vector n. Finally,a second orthogonal unit vector in the plane of the triangle is be com- puted as 02=n X U1 (8) The above developments have been performed based on the current configura- tion.However,reference quantities may be deduced by replacing lower case letters by upper case ones (e.g.,vi -Vi,etc.)
Finite Element Analysis of Membrane Structures 49 Fig. 1. Description of coordinates for triangular element Placing the origin of the YI and yi coordinates at nodal location X˜ 1 and x˜1, respectively, the unit base vectors may be constructed from the linear displacement triangle by aligning the first vector along the 1 − 2 side. Accordingly, we define the first unit vector as v1 = x˜2 − x˜1 x˜2 − x˜1 = ∆x˜21 ∆x˜21 (5) where ∆x˜ij = x˜i − x˜j and a = aT a 1/2 . Next a vector normal to the triangle is constructed as v3 = ∆x˜21 × ∆x˜31 (6) and normalized to a unit vector as n = v3 v3 (7) The vector v3 plays a special role in later development of nodal forces for follower pressure loading as it is twice the area of the triangle times the unit normal vector n. Finally, a second orthogonal unit vector in the plane of the triangle is be computed as v2 = n × v1 . (8) The above developments have been performed based on the current configuration. However, reference quantities may be deduced by replacing lower case letters by upper case ones (e.g., v1 → V 1, etc.). x2 x3 x1 1 2 3 y1 y2 n
50 Robert L.Taylor,Eugenio Ofate and Pere-Andreu Ubach With the above base vectors defined,positions in the plane of the triangle may be given directly as 班=(x-)·v (9) In general,an interpolation may be given as y=￡a9 (10) We note from Eq.(9)that is identically zero hence Eq.(10)reduces to y=522+33. (11) However,the constraint (4)still restricts the admissible values for 2 and 3. 2.1 Deformation Gradient From the above description of the motion of the triangle it is now possible to deduce the deformation gradient in the plane of the triangle as F= dy Ψ =1+ (12) Using the parametric representations(11)we can compute the deformation gradient from ay ay aY ay Oy OE -FOEE (13) If we define the arrays J and j as J= ay and y OE j=0 (14) then the deformation gradient is given by F=jJ-1 (15) In the above J is the Jacobian transformation for the reference frame and j that for the current frame.Expanding the relations for each Jacobian we obtain (16) and j= (△221)Tu1(421)T1 (Ai2)Tv2 ()Tv2 (17) By noting that is orthogonal to V2 and similarly for the current configu- ration that 2-is orthogonal to v2 and in addition using the definition for Vi
50 Robert L. Taylor, Eugenio O˜nate and Pere-Andreu Ubach ˜ With the above base vectors defined, positions in the plane of the triangle may be given directly as yi = x − x˜1 · vi (9) In general, an interpolation may be given as y = ξαy˜α (10) We note from Eq. (9) that y˜1 is identically zero hence Eq. (10) reduces to y = ξ2 y˜2 + ξ3 y˜3 . (11) However, the constraint (4) still restricts the admissible values for ξ2 and ξ3. 2.1 Deformation Gradient From the above description of the motion of the triangle it is now possible to deduce the deformation gradient in the plane of the triangle as F = ∂y ∂Y = 1 + ∂u ∂Y (12) Using the parametric representations (11) we can compute the deformation gradient from ∂y ∂Y ∂Y ∂ξ = F ∂Y ∂ξ = ∂y ∂ξ (13) If we define the arrays J and j as J = ∂Y ∂ξ and j = ∂y ∂ξ (14) then the deformation gradient is given by F = j J−1 (15) In the above J is the Jacobian transformation for the reference frame and j that for the current frame. Expanding the relations for each Jacobian we obtain J = ⎡ ⎣ ∆X˜ 21T V 1 ∆X˜ 31T V 1 ∆X˜ 21T V 2 ∆X˜ 31T V 2 ⎤ ⎦ (16) and j = ∆x˜21T v1 ∆x˜31T v1 ∆x˜21T v2 ∆x˜31T v2 (17) By noting that X˜ 2 − X˜ 1 is orthogonal to V 2 and similarly for the current configuration that x˜2 − x˜1 is orthogonal to v2 and in addition using the definition for V i
Finite Element Analysis of Membrane Structures 51 and v;the above simplify to - (18) and -[.a2a (19) lv3/川△金川 Using these definitions,the right Cauchy-Green deformation tensor may be ex- panded as C=FTF=J-TjTjJ-1=GTgG (20) where G is used to denote the inverse of J.In component form we have 1 2201911912]「J2-☑2 c=-加hi]9e9e0hi (21) in which the terms in the kernel array involving j may be expressed in the particu- larly simple form 911= =(4221)T4221 912-j12i1=(4221)T△231 (22) 922=+品2=(4元31)T4立31 2.2 Material Constitution-Elastic Behavior In the present work we assume that a simple St.Venant-Kirchhoff material model may be used to express the stresses from the deformations.Stresses are thus given by S=DE (23) where D are constant elastic moduli.and the Green-Lagrange strains E are given in terms of the deformation tensor as E=(C-1) (24) In each triangular element the deformation may be computed from (19)to (22), thus giving directly the stress. 3 Weak Form for Equations of Motion A weak form for the membrane may be written using a virtual work expression given by 6= δri po h:hd2+ 6xi co iid+ 6EISIs hds 2 (25) xibi dw- 6xi tidy=0
Finite Element Analysis of Membrane Structures 51 and vi the above simplify to J = ⎡ ⎣ ∆X˜ 21 , ∆X˜ 21T ∆X˜ 31 /∆X˜ 21 0 , V 3/∆X˜ 21 ⎤ ⎦ (18) and j = ∆x˜21 , ∆x˜21T ∆x˜31 /∆x˜21 0 , v3/∆x˜21 (19) Using these definitions, the right Cauchy-Green deformation tensor may be expanded as C = F T F = J −T jT jJ−1 = GT gG (20) where G is used to denote the inverse of J. In component form we have C = 1 J 2 11J 2 22 J22 0 −J12 J11 g11 g12 g12 g22 J22 −J12 0 J11 (21) in which the terms in the kernel array involving j may be expressed in the particularly simple form g11 = j2 11 = ∆x˜21T ∆x˜21 g12 = j12j11 = ∆x˜21T ∆x˜31 g22 = j2 12 + j2 22 = ∆x˜31T ∆x˜31 (22) 2.2 Material Constitution - Elastic Behavior In the present work we assume that a simple St.Venant-Kirchhoff material model may be used to express the stresses from the deformations. Stresses are thus given by S = D E (23) where D are constant elastic moduli. and the Green-Lagrange strains E are given in terms of the deformation tensor as E = 1 2 (C − I) . (24) In each triangular element the deformation may be computed from (19) to (22), thus giving directly the stress. 3 Weak Form for Equations of Motion A weak form for the membrane may be written using a virtual work expression given by δΠ = Ω δxi ρ0 h x¨i h dΩ + Ω δxi c0 x˙ i dΩ + Ω δEIJ SIJ h dΩ − Ω δxibi dω − γt δxi t ¯i dγ = 0 (25)