《纺织复合材料》课程参考文献（Textile Composites and Inflatable Structures）06 FE Analysis of Membrane Systems Including Wrinkling and Coupling

FE Analysis of Membrane Systems Including Wrinkling and Coupling Riccardo Rossi,Vitaliani Renatol,and Eugenio Onate2 1 Universita di Padova ricrossi@caronte.dic.unipd.it 2 CIMNE,UPC,Barcelona onate cimne.upc.es 1 Abstract Current work summarizes the experience of the writer in the modeling of membrane systems.A first subsection describes an efficient membrane model, together with a reliable solution procedure.The following section addresses the simulation of the wrinkling phenomena providing details of a new solution procedure.The last one proposes an efficient technique to obtain the solution of the fluid structural interaction problem. 2 The Membrane Model A membrane is basically a 2D solid which "lives"in a 3D environment.Given the lack of flexural stiffness,membranes can react to applied load only by using their in-plane resistance "choosing"the spatial disposition that is best suited to resist to the external forces.The consequence is that membrane structures tend naturally to find the optimal shape(compatible with the ap- plied constraints)for any given load.In this "shape research",they typically undergo large displacements and rotations. From a numerical point of view,this reflects an intrinsic geometrical non- linearity that has to be taken in account in the formulation of the finite element model.In particular,an efficient Membrane Element should be able to represent correctly arbitrary rotations both of the element as a whole and internally to each element.The possibility of unrestricted rigid body motions constitutes a source of ill-conditioning or even of singularity of the tangent stiffness matrix introducing the need of carefully designed solution procedures. 2.1 Finite Element Model Current section describes a finite element model that meets all of the require- ments for the correct simulation of general membrane systems.The derivation 89 E.Onate and B.Kroplin (eds.).Textile Composites and Inflatable Structures,89-108. C 2005 Springer.Printed in the Netherlands

FE Analysis of Membrane Systems Including Wrinkling and Coupling Riccardo Rossi 1, Vitaliani Renato1, and Eugenio Onate2 1 Universit`a di Padova ricrossi@caronte.dic.unipd.it 2 CIMNE, UPC, Barcelona onate@cimne.upc.es 1 Abstract Current work summarizes the experience of the writer in the modeling of membrane systems. A first subsection describes an efficient membrane model, together with a reliable solution procedure. The following section addresses the simulation of the wrinkling phenomena providing details of a new solution procedure. The last one proposes an efficient technique to obtain the solution of the fluid structural interaction problem. 2 The Membrane Model A membrane is basically a 2D solid which “lives” in a 3D environment. Given the lack of flexural stiffness, membranes can react to applied load only by using their in–plane resistance “choosing” the spatial disposition that is best suited to resist to the external forces. The consequence is that membrane structures tend naturally to find the optimal shape (compatible with the ap￾plied constraints) for any given load. In this “shape research”, they typically undergo large displacements and rotations. From a numerical point of view, this reflects an intrinsic geometrical non– linearity that has to be taken in account in the formulation of the finite element model. In particular, an efficient Membrane Element should be able to represent correctly arbitrary rotations both of the element as a whole and internally to each element. The possibility of unrestricted rigid body motions constitutes a source of ill-conditioning or even of singularity of the tangent stiffness matrix introducing the need of carefully designed solution procedures. 2.1 Finite Element Model Current section describes a finite element model that meets all of the require￾ments for the correct simulation of general membrane systems. The derivation 89 E. Oñate and B. Kröplin (eds.), Textile Composites and Inflatable Structures, 89–108. © 2005 Springer. Printed in the Netherlands

90 Riccardo Rossi,Vitaliani Renato,and Eugenio Onate makes use exclusively of orthogonal bases simplifying the calculations and al- lowing to express all the terms directly in Voigt Notation,which eases the implementation. Einstein's sum on repeated indices is assumed unless specified otherwise .xI=,is the position vector of the I-th node in the cartesian space(3D space) .={describes the position of a point in the local system of coordi- nates Capital letter are used for addressing to the reference configuration Nr(g)is the value of the shape function centered on node I on the point of local coordinatesξ The use of the standard iso-parametric approach allows to express the position of any point as x()=NI()xI. In the usual assumptions of the continuum mechanics it is always possible to define the transformation between the local system of coordinates and the cartesian system as 长+，P-长，P-0x5=gG ∂灰 (1) {长，7+dr-长，nr-05l=g, (2) in which we introduced the symbols gE= aN,x灯 (3) DE gn= aN,dx灯 8n (4) the vectors ge and gn of the 3D space,can be considered linearly independent (otherwise compenetration or self contact would manifest)follows immedi- ately that they can be used in the construction of a base of the 3D space.In particular an orthogonal base can be defined as v1二ge gE (5) n= g:X gn →V2=nXV1 (6) lge×gml V3=g X gn (7) Vectors vI and v2 describe the local tangent plane to the membrane while the third base vector is always orthogonal.Follows the possibility of defining a local transformation rule that links the local coordinates and the coordinates

90 Riccardo Rossi, Vitaliani Renato, and Eugenio Onate makes use exclusively of orthogonal bases simplifying the calculations and al￾lowing to express all the terms directly in Voigt Notation, which eases the implementation. Einstein’s sum on repeated indices is assumed unless specified otherwise • xI = {xI , yI , zI } T is the position vector of the I–th node in the cartesian space (3D space) • ξ = {ξ, η} T describes the position of a point in the local system of coordi￾nates • Capital letter are used for addressing to the reference configuration • NI (ξ) is the value of the shape function centered on node I on the point of local coordinates ξ The use of the standard iso-parametric approach allows to express the position of any point as x(ξ) = NI (ξ)xI. In the usual assumptions of the continuum mechanics it is always possible to define the transformation between the local system of coordinates and the cartesian system as {ξ + dξ, η} T − {ξ, η} T → ∂x(ξ,η) ∂ξ dξ = gξdξ (1) {ξ, η + dη}T − {ξ, η} T → ∂x(ξ, η) ∂η dη = gηdη (2) in which we introduced the symbols gξ = ∂NI (ξ, η) ∂ξ xI (3) gη = ∂NI (ξ, η) ∂η xI (4) the vectors gξ and gη of the 3D space, can be considered linearly independent (otherwise compenetration or self contact would manifest) follows immedi￾ately that they can be used in the construction of a base of the 3D space. In particular an orthogonal base can be defined as v1 = gξ gξ (5) n = gξ × gη gξ × gη → v2 = n × v1 (6) v3 = gξ × gη (7) Vectors v1 and v2 describe the local tangent plane to the membrane while the third base vector is always orthogonal. Follows the possibility of defining a local transformation rule that links the local coordinates ξ and the coordinates

FE Analysis of Membrane Systems Including Wrinkling and Coupling 91 x in the local tangent plane base.This can be achieved by considering that an increment fde,0}maps to the new base as ge·V1 gm·V1 dn (8) 0 d切 gV2 this is sinthetized by the definition of the linear application d⑦ ge·V1gm●V1 d在2 (9) ge●V2gm●V2 dn →d=jd延 it should be noted that ge.v3 and gn.v3 are identically zero consequently no components are left apart in representing the membrane in the new co- ordinates system.Taking into account the definition of the base vectors the tensor j becomes (after some calculations) (10) and its determinant det)=lv3l=lgs×gl (11) As the interest is focused on the purely membranal behavior,it is not needed to take in account the deformation of the structure over the thickness,as this can be calculated "a posteriori"once known the deformation of the mid-plane. On the base of this observation,the deformation gradient which describes the deformation of the membrane as a 3D body dx F3x3=8区 (12) can be replaced by 宜2x2= 0胶0胶05 胶=c胶 =j订-1 (13) taking in account the behavior over the thickness in the definition of the (two dimensional)constitutive model to be used (for example making the assumption of plane stress).The symbol J is used here and in the following to indicate j calculated in the reference position. under this considerations the subsequent development of the finite element follows closely the standard procedure for a non-linear 2D finite element,the only difference being that the local base changes on all the domain,which makes the linearization slightly more involved. To proceed further we need therefore to write the Right Cauchy Green Strain tensor C=FTF which takes the form C=(J-TjTjJ-1)=(GTgG) (14)

FE Analysis of Membrane Systems Including Wrinkling and Coupling 91 x, in the local tangent plane base. This can be achieved by considering that an increment {dξ, 0} maps to the new base as  dξ 0  →  gξ • v1 gξ • v2  dξ ;  0 dη  →  gη • v1 gη • v2  dη (8) this is sinthetized by the definition of the linear application  dx,1 dx,2  =  gξ • v1 gη • v1 gξ • v2 gη • v2  dξ dη → dx, = jdξ (9) it should be noted that gξ • v3 and gη • v3 are identically zero consequently no components are left apart in representing the membrane in the new co￾ordinates system. Taking into account the definition of the base vectors the tensor j becomes (after some calculations) j =  gξ gξ•gη
gξ
0
v3

gξ
(10) and its determinant det(j) = v3 = gξ × gη (11) As the interest is focused on the purely membranal behavior, it is not needed to take in account the deformation of the structure over the thickness, as this can be calculated “a posteriori” once known the deformation of the mid–plane. On the base of this observation, the deformation gradient which describes the deformation of the membrane as a 3D body F3×3 = ∂x ∂X (12) can be replaced by F,2×2 = ∂x, ∂X, = ∂x, ∂ξ ∂ξ ∂X, = jJ−1 (13) taking in account the behavior over the thickness in the definition of the (two dimensional) constitutive model to be used (for example making the assumption of plane stress). The symbol J is used here and in the following to indicate j calculated in the reference position. under this considerations the subsequent development of the finite element follows closely the standard procedure for a non–linear 2D finite element, the only difference being that the local base changes on all the domain, which makes the linearization slightly more involved. To proceed further we need therefore to write the Right Cauchy Green Strain tensor C = FT F which takes the form C =  J−T j T jJ−1 =  GT gG (14)

92 Riccardo Rossi,Vitaliani Renato,and Eugenio Onate where we introduced the symbols g=jTj and G=J-1.Operator g takes, after some calculations,the simple form gE●gEgm·g g= (15) ge●gngm●gm From the definition of he Green Lagrange strain tensor E =(C-I) we obtain immediately E =6C.This allows to write the equation of vir- tual works in compact form as (taking in consideration only body forces and pressure forces) SWint 6Wert 6Wpress (16) 、 6C:S=ho x●b+ p6x●n (17) Internal Work The termC:S describes the work of internal forces during the defor- mation process.Operator G=J-1 is referred to the reference configuration and is therefore strictly constant,follows immediately that C=GTogG (18) The term 6C:S becomes in Einstein's notation 3iC:8-i0u51=06m5m 1 (19) 2 introducing the symbols 6g11 S11 29则→26g}= 6g22 SJ→{S}= S22 (20) 2ǒg12 S12 (G11)2(G12)22G11G12 GuGJ→[Q1T 0 (G22)2 0 (21) 0 G12G22G11G22 it is possible to express the (19)in Voigt form as ic:5-0g1rQrs)-gr母：o=Qrs (22) considering the definition(15),introducing the symbol x)(x)..x and taking in account the isoparametric approximation one obtains Nx·g影=6x {g 2g11= (23) {g}/

92 Riccardo Rossi, Vitaliani Renato, and Eugenio Onate where we introduced the symbols g = jT j and G = J−1. Operator g takes, after some calculations, the simple form g =  gξ • gξ gη • gξ gξ • gη gη • gη  (15) From the definition of he Green Lagrange strain tensor E = 1 2 (C − I) we obtain immediately δE = 1 2 δC. This allows to write the equation of vir￾tual works in compact form as (taking in consideration only body forces and pressure forces) δWint = δWext + δWpress (16) h0 2  Ω δC : S = h0  Ω δx • b +  ω pδx • n (17) Internal Work The term h0 2 ) Ω δC : S describes the work of internal forces during the defor￾mation process. Operator G = J−1 is referred to the reference configuration and is therefore strictly constant, follows immediately that δC = GT δgG (18) The term δC : S becomes in Einstein’s notation 1 2 δC : S = 1 2 δCIJ SIJ = 1 2 δgijGiIGjJ SIJ (19) introducing the symbols 1 2 δgij → 1 2 δ {g} = 1 2 ⎛ ⎝ δg11 δg22 2δg12 ⎞ ⎠ ; SIJ → {S} = ⎛ ⎝ S11 S22 S12 ⎞ ⎠ (20) GiIGjJ → [Q] T = ⎛ ⎝ (G11)2 (G12)2 2G11G12 0 (G22)2 0 0 G12G22 G11G22 ⎞ ⎠ (21) it is possible to express the (19) in Voigt form as 1 2 δC : S = 1 2 {δg}T [Q] T {S} = 1 2 {δg}T {s} ; {s} = [Q] T {S} (22) considering the definition(15), introducing the symbol {δx}T=  {δx1} T . . . {δxk} T  and taking in account the isoparametric approximation one obtains 1 2 δg11 = ∂NI ∂ξ δxI • gξ = {δx}T ⎛ ⎝ ∂N1 ∂ξ {gξ} ... ∂Nk ∂ξ {gξ} ⎞ ⎠ (23)

FE Analysis of Membrane Systems Including Wrinkling and Coupling 93 ga-21·8=刘 {g} (24) ON (gn) 202g12= Nx灯·g<+0c Nx·g={6x {g}+盼{ge) (25) on 股{g}+册{8<} by defining the matrix /器{g}器{g}器{g}+器{g} [b)T (26) 驶{g<}{g}{g}+{g} it is then possible to write acy7s)=EΨs=7rQrs (27) Defining the symbol B] [B]=[Q][b] (28) we finally obtain (fint}=ho [B]T (S)do (29) 6Wint={ox){fint} (30) External Work Derivation of the work of external conservative forces follows the standard procedure and can be found on any book on nonlinear finite elements.The expression of the work of follower forces (body forces)is on the other hand a little more involved.In the following the pressure is considered constant,the non linearity being introduced by the change of direction of the normal.For the derivation of the pressure contributions it is much easier to perform the integration over the actual domain then over the reference one. 6Wor p6x.nd p6x.ndet(j)dedn (31) taking in account the definition of the base vectors (6)(7),and considering (11)we obtain immediately }= Ni(,m)p(ξ，jv3(5,n)dd (32) }=（《6}T.…{.}T） (33) 6Wpr ={6x}(fpr} (34)

FE Analysis of Membrane Systems Including Wrinkling and Coupling 93 1 2 δg22 = ∂NI ∂ξ δxI • gη = {δx} T ⎛ ⎝ ∂N1 ∂η {gη} ... ∂Nk ∂η {gη} ⎞ ⎠ (24) 1 2 δ2g12 = ∂NI ∂η δxI • gξ + ∂NI ∂ξ δxI • gη = {δx}T ⎛ ⎝ ∂N1 ∂ξ {gη} + ∂N1 ∂η {gξ} ... ∂Nk ∂ξ {gη} + ∂Nk ∂η {gξ} ⎞ ⎠ (25) by defining the matrix [b] T = ⎛ ⎝ ∂N1 ∂ξ {gξ} ∂N1 ∂η {gη} ∂N1 ∂ξ {gη} + ∂N1 ∂η {gξ} ... ... ... ∂Nk ∂ξ {gξ} ∂Nk ∂η {gη} ∂Nk ∂ξ {gη} + ∂Nk ∂η {gξ} ⎞ ⎠ (26) it is then possible to write 1 2 {δC}T {S} = {δE}T {S} = {δx}T [b] T [Q] T {S} (27) Defining the symbol [B] [B] = [Q] [b] (28) we finally obtain {fint} =  Ω h0 [B] T {S} dΩ (29) δWint = {δx} T {fint} (30) External Work Derivation of the work of external conservative forces follows the standard procedure and can be found on any book on nonlinear finite elements. The expression of the work of follower forces (body forces) is on the other hand a little more involved. In the following the pressure is considered constant, the non linearity being introduced by the change of direction of the normal. For the derivation of the pressure contributions it is much easier to perform the integration over the actual domain then over the reference one. δWpr =  ω pδx • ndω =  ξ,η pδx • ndet(j)dξdη (31) taking in account the definition of the base vectors (6)(7), and considering (11) we obtain immediately {fI} =  ξ,η NI (ξ, η)p(ξ, η)v3(ξ, η)dξdη (32) {fpr} =  {f1}T . . . {fk}T T (33) δWpr = {δx} T {fpr} (34)

94 Riccardo Rossi,Vitaliani Renato,and Eugenio Onate Linearization Equation (16)is nonlinear,its practical use needs therefore its lineariza- tion.The best rate of convergence is theoretically given by Newton-Raphson technique which guarantees quadratical convergence to the solution.Defining =6Wint -oWert -oWpr each Newton-Raphson step takes the form d亚+亚=0 (35) The term can be explicitated using expression(30)(34)we therefore miss only the differential dy that can be evaluated from the linearization of the different contributions Linearization of internal work The term connected to the internal works can be linearized as follows d(Wint)= ( 6C:S 名 =空a:s+人cds6 the first terms gives,by using(22) 告46c.s=空/(g）= afog)d(x)(s) 8{x} (37) now it can be seen that a(（行og1r）=(dgcdg,gdgs+gds,)问= =(s16ge·dges226gn·dgns12(⑥gn·dge+ge·dgn) (38) substitution of the shape functions gives immediately a set of equalities in the form 2NN,d=u2086网 s1dge·dg=9110诞0c which makes possible to write On on 512 (40)

94 Riccardo Rossi, Vitaliani Renato, and Eugenio Onate Linearization Equation (16) is nonlinear, its practical use needs therefore its lineariza￾tion. The best rate of convergence is theoretically given by Newton-Raphson technique which guarantees quadratical convergence to the solution. Defining Ψ = δWint − δWext − δWpr each Newton–Raphson step takes the form dΨ + Ψ = 0 (35) The term Ψ can be explicitated using expression (30)(34) we therefore miss only the differential dΨ that can be evaluated from the linearization of the different contributions Linearization of internal work The term connected to the internal works can be linearized as follows d (Wint) = d h0 2  Ω δC : S  == h0 2  Ω d (δC) : S + h0 2  Ω δC : d (S) (36) the first terms gives, by using (22) h0 2  Ω d (δC) : S = h0 2  Ω d  {δg}T  {s} = = h0 2  Ω d  ∂ {δg}T ∂ {x} d {x} {s} (37) now it can be seen that d 1 2 {δg} T  =  δgξ • dgξ δgη • dgη δgη • dgξ + δgξ • dgη  {s} = =  s11δgξ • dgξ s22δgη • dgη s12 (δgη • dgξ + δgξ • dgη)  (38) substitution of the shape functions gives immediately a set of equalities in the form s11δgξ • dgξ = s11 ∂NI ∂ξ ∂NJ ∂ξ δijδxI • dxjJ = s11 ∂NI ∂ξ ∂NJ ∂ξ δijδxiIdxjJ (39) which makes possible to write d 1 2 {δg} T  =  s11 ∂NI ∂ξ ∂NJ ∂ξ + s22 ∂NI ∂η ∂NJ ∂η + s12 ∂NI ∂η ∂NJ ∂ξ + ∂NI ∂ξ ∂NJ ∂η  δij δxiIdxjJ (40)

FE Analysis of Membrane Systems Including Wrinkling and Coupling 95 introducing the vectors a=(.×）andb=(…跚) together with the new tensor A=(s11aa+s22bb+s12(ba+ab));aa=aa (41) we can simplify greatly the expression in the form Arj6ijoxirdxjs (42) or,in Voigt form (传er）=(cP6r(如田） =16x)T Kgeo]{dx}(43) dxk the derivation of the "material"contribution to the stiffness matrix follows the standard path.Assuming as normally as dS= E:dE→{ds}=Daml{dE吲 (44) we obtain immediately D)(dx) =6{x)Kmat]{dx} (45) [KhoD]d (46) Linearization of external forces The linearization of the work Wert is not needed as it describes the work of constant forces.The only term missing is therefore the one relative to the work of the follower forces a([nix nd）-iPa面=iPa(ngs×g}do）an differentiating pNI [gsx gn}we obtain d(pN{ge×g)=pN,{dge×gn}+pNi{ge×dgn}= =pNi{ge×dgn}-pNi{gn×dge} (48)

FE Analysis of Membrane Systems Including Wrinkling and Coupling 95 introducing the vectors a =  ∂N1 ∂ξ . . . ∂Nk ∂ξ  and b =  ∂N1 ∂η ... ∂Nk ∂η  together with the new tensor A = (s11aa + s22bb + s12 (ba + ab)) ; aa = a ⊗ a (41) we can simplify greatly the expression in the form d 1 2 {δg}T  = AIJ δijδxiIdxjJ (42) or, in Voigt form d 1 2 {δg}T  =  {δx1}T . . . {δxk}T   A11 [I] ... A1k [I] . .. ... ... Ak1 [I] ... Akk [I]  ⎛ ⎝ {dx1} ... {dxk} ⎞ ⎠ = {δx}T [Kgeo] {dx}(43) the derivation of the “material” contribution to the stiffness matrix follows the standard path. Assuming as normally dS = ∂S ∂E : dE → {dS} = [Dtan] {dE} (44) we obtain immediately  Ω h0 2 δC : d (S) =  Ω h0δ {x}T [B] T [Dtan] [B] dΩ {dx} = δ {x}T [Kmat] {dx} (45) [Kmat] =  Ω h0 [B] T [Dtan] [B] dΩ (46) Linearization of external forces The linearization of the work Wext is not needed as it describes the work of constant forces. The only term missing is therefore the one relative to the work of the follower forces d  ω pδx • ndω → δ {xI} T dfI = δ {xI}T d  ω pNI {gξ × gη} dω (47) differentiating pNI {gξ × gη}we obtain d (pNI {gξ × gη}) = pNI {dgξ × gη} + pNI {gξ × dgη} = = pNI {gξ × dgη} − pNI {gη × dgξ} (48)

96 Riccardo Rossi,Vitaliani Renato,and Eugenio Onate Considering that it is possible to write the cross product of two vectors in Voigt format as C 0 -a3 a2 c=ax b- C2 a3 -a1 b →{c}=[ax]{b} (49) C3 -a2a1 0 and taking in account (3)and(4)we obtain d(pN,{ge×gn}) (50) K11… [Kpr] Kk1... Kkk Kdl=(wk-n (51) d(oWpr)={ox)Kpr]{dx} (52) Linearized formulation The only step missing is to merge all the terms in (35)to find the final expression.The result of this operation is {ox)([Kgeo]+[Kmat]-[Kprl){dx}={ox)({fext}-{fint}) (53) invoking the arbitrariety of fox and introducing the definitions Ktan]Kgeo]+Kmat]-Kpr] (54) {R}={fext}-{fint} (55) the principle of virtual works gives for each element [Ktan]{dx}=(R) (56) 2.2 Solution Procedure As briefly outlined at the beginning of the section,membrane systems are possibly subjected to large rigid body motions which reflects in singular or ill-conditioned "static"stiffness matrices.In addition,convergence of the Newton-Raphson algorithm is often difficult as the final solution can be very "far"from the initial guess even for little variations of the applied loads. Dynamic solution techniques on the other hand are not affected by such problems.Mass and damping contributions remove the singularities from the system and generally provide a better conditioning to the problem.The intro- duction of dynamic terms provides as well an excellent source of stabilization for the solution (physically the solution can't change much in a small time), ending up with better convergence properties inside each solution step

96 Riccardo Rossi, Vitaliani Renato, and Eugenio Onate Considering that it is possible to write the cross product of two vectors in Voigt format as c = a × b → ⎛ ⎝ c1 c2 c3 ⎞ ⎠ = ⎛ ⎝ 0 −a3 a2 a3 0 −a1 −a2 a1 0 ⎞ ⎠ ⎛ ⎝ b1 b2 b3 ⎞ ⎠ → {c} = [a×] {b} (49) and taking in account (3) and (4) we obtain d (pNI {gξ × gη}) =  pNI ∂NJ ∂η [gξ×] − pNI ∂NJ ∂ξ [gξ×]  {dxJ} (50) [Kpr] = ⎛ ⎝ K11 ... K1k ... ... ... Kk1 . . . Kkk ⎞ ⎠ ; [KIJ] =  pNI ∂NJ ∂η [gξ×] − pNI ∂NJ ∂ξ [gξ×]  (51) d(δWpr) = {δx} T [Kpr] {dx} (52) Linearized formulation The only step missing is to merge all the terms in (35) to find the final expression. The result of this operation is {δx} T ([Kgeo] + [Kmat] − [Kpr]) {dx} = {δx}T ({fext} − {fint}) (53) invoking the arbitrariety of {δx} and introducing the definitions [Ktan] = [Kgeo] + [Kmat] − [Kpr] (54) {R} = {fext} − {fint} (55) the principle of virtual works gives for each element [Ktan] {dx} = {R} (56) 2.2 Solution Procedure As briefly outlined at the beginning of the section, membrane systems are possibly subjected to large rigid body motions which reflects in singular or ill-conditioned “static” stiffness matrices. In addition, convergence of the Newton–Raphson algorithm is often difficult as the final solution can be very “far” from the initial guess even for little variations of the applied loads. Dynamic solution techniques on the other hand are not affected by such problems. Mass and damping contributions remove the singularities from the system and generally provide a better conditioning to the problem. The intro￾duction of dynamic terms provides as well an excellent source of stabilization for the solution (physically the solution can’t change much in a small time), ending up with better convergence properties inside each solution step

FE Analysis of Membrane Systems Including Wrinkling and Coupling 97 Any standard (non-linear)time integration technique can be theoretically used in conjunction with the proposed FE model for the study of dynamic response of the systems of interest.Some care should be however taken in the choice because the high geometric non-linearities tend to challenge the stability of the time-integration scheme chosen. Generally speaking,"statics"can be seen as the limit to which a dynamic process tends (under a given constant load).Dynamic systems show a "tran- sient"phase that vanishes in time to reach the so called "steady state";the presence of damping in the system reduces gradually the oscillations making the system tend to a constant configuration that is the "static"solution.The time needed for the system to reach this final configuration is controlled by the amount of damping.For values of system's damping exceeding a critical value,the transient phase disappears and the systems reaches directly the final solution without any oscillation. In many situations the main engineering interest is focused on "static" solutions rather than on the complete dynamic analysis of the system.The previews considerations suggests immediately that "statics"could be obtained efficiently by studying the dynamics of over damped systems.This could be obtained by simply adding a fictitious damping source to the "standard"dy- namic problem.The "only"problem is therefore the choice of an idoneous form for such damping.Unfortunately this choice is not trivial,however it possible to observe [10],[8 that the "steady state"solution of the system M+D+Kx=f (x) (57) is (statically)equivalent to that of the system D+Kx=f(x) (58) which can be seen as the previews for the case of zero density.The advantage of this equivalent system is that the inertia terms are always zero,consequently the system converges smoothly in time to its solution.This final solution is not affected by the particular choice of the damping,however in the author's experience,an effective choice is D=BM as proposed by [10]. Table (1)gives the details of the proposed solution procedure,making use of Newmark's integration scheme.The procedure described differs from a"real"dynamics simulation only on the choice of the damping and of the mass matrix.Any other choice is possible for the time integration scheme to be used.It is of interest to observe that the system described is highly dissipative,energy stability of the time integration scheme is therefore not crucial. 3 Wrinkling Simulation Given the lack of flexural stiffness,membrane systems are easily subjected to buckling in presence of any compressive load.The idea is that when a

FE Analysis of Membrane Systems Including Wrinkling and Coupling 97 Any standard (non–linear) time integration technique can be theoretically used in conjunction with the proposed FE model for the study of dynamic response of the systems of interest. Some care should be however taken in the choice because the high geometric non–linearities tend to challenge the stability of the time–integration scheme chosen. Generally speaking, “statics” can be seen as the limit to which a dynamic process tends (under a given constant load). Dynamic systems show a “tran￾sient” phase that vanishes in time to reach the so called “steady state”; the presence of damping in the system reduces gradually the oscillations making the system tend to a constant configuration that is the “static” solution. The time needed for the system to reach this final configuration is controlled by the amount of damping. For values of system’s damping exceeding a critical value, the transient phase disappears and the systems reaches directly the final solution without any oscillation. In many situations the main engineering interest is focused on “static” solutions rather than on the complete dynamic analysis of the system. The previews considerations suggests immediately that “statics” could be obtained efficiently by studying the dynamics of over damped systems. This could be obtained by simply adding a fictitious damping source to the “standard” dy￾namic problem. The “only” problem is therefore the choice of an idoneous form for such damping. Unfortunately this choice is not trivial, however it possible to observe [10],[8] that the “steady state” solution of the system Mx¨ + Dx˙ + Kx = f (x) (57) is (statically) equivalent to that of the system Dx˙ + Kx = f (x) (58) which can be seen as the previews for the case of zero density. The advantage of this equivalent system is that the inertia terms are always zero, consequently the system converges smoothly in time to its solution. This final solution is not affected by the particular choice of the damping, however in the author’s experience, an effective choice is D = βM as proposed by [10]. Table (1) gives the details of the proposed solution procedure, making use of Newmark’s integration scheme. The procedure described differs from a “real” dynamics simulation only on the choice of the damping and of the mass matrix. Any other choice is possible for the time integration scheme to be used. It is of interest to observe that the system described is highly dissipative, energy stability of the time integration scheme is therefore not crucial. 3 Wrinkling Simulation Given the lack of flexural stiffness, membrane systems are easily subjected to buckling in presence of any compressive load. The idea is that when a

98 Riccardo Rossi,Vitaliani Renato,and Eugenio Onate Table 1.Pseudo-Static solution procedure for pseudo-static strategy:calculate the constant matrices D=BM set M =0 after initializing the damping matrix.(if M is not set to 0,"real" dynamic simulation can be performed) choose Newmark constants:a classical choice is 1 1 6=2；a=4 evaluate the constants a41 i a1=6 1 1 a0= -at i 02=adt s=六-1：a=日-1；s=兰（很-2） predict the solution at time t+At using for example x8+a1=XL+文△t x+△t= :=0 iterate until convergence -calculate the system's contributions Key=[Ktan]+do[M]+1 [D] Ray")=(R)-[M][D] solve the system for the correction dx update the results as x=xitar+dx Ax=x牛4-x t+At=a1△x-a4xt-a5戈t 戈4+At=a0△x-a2:-a3戈： go to next time step compressive stress tends to appear on a part of a structure,it is immediately removed by local instability phenomena,that manifest with the formation of little"waves"of direction perpendicular to the direction of stresses.Predic- tion of the size of those "waves"commonly called "wrinkles"is not generally possible as their disposition is somehow random and connected to initial im- perfections.However their average size is strictly connected to the bending stiffness meaning in particular that for the problems of interest the wrinkle

98 Riccardo Rossi, Vitaliani Renato, and Eugenio Onate Table 1. Pseudo–Static solution procedure • for pseudo–static strategy: calculate the constant matrices D = βM set M = 0 after initializing the damping matrix. (if M is not set to 0, “real” dynamic simulation can be performed) • choose Newmark constants: a classical choice is δ = 1 2 ; α = 1 4 • evaluate the constants a0 = 1 α∆t2 ; a1 = δ α∆t ; a2 = 1 α∆t a3 = 1 2α − 1 ; a4 = δ α − 1 ; a5 = ∆t 2  δ α − 2  • predict the solution at time t + ∆t using for example x0 t+∆t = xt + x˙ t∆t x˙ t+∆t = x˙ t x¨t+∆t = 0 • iterate until convergence – calculate the system’s contributions  Kdyn tan  = [Ktan] + a0 [M] + a1 [D] 1 Rdyn 2 = {R} − [M] 1 x¨i t+∆t 2 − [D] 1 x˙ i t+∆t 2 – solve the system for the correction dx – update the results as xi+1 t+∆t = xi t+∆t + dx ∆x = xi+1 t+∆t − xt x˙ t+∆t = a1∆x − a4x˙ t − a5x¨t x¨t+∆t = a0∆x − a2x˙ t − a3x¨t • go to next time step compressive stress tends to appear on a part of a structure, it is immediately removed by local instability phenomena, that manifest with the formation of little ”waves” of direction perpendicular to the direction of stresses. Predic￾tion of the size of those ”waves” commonly called ”wrinkles” is not generally possible as their disposition is somehow random and connected to initial im￾perfections. However their average size is strictly connected to the bending stiffness meaning in particular that for the problems of interest the wrinkle