Implementation of a simple wrinkling model into argyris’ membrane finite element

142 VI International Conference on Textile Composites and Inflatable Structures STRUCTURAL MEMBRANES 2013 K.-U.Bletzinger, B. Kröplin and E. Oñate (Eds) IMPLEMENTATION OF A SIMPLE WRINKLING MODEL INTO ARGYRIS’ MEMBRANE FINITE ELEMENT RUY M.O. PAULETTI*, DANIEL M. GUIRARDI† *Polytechnic School of the University of São Paulo P.O. Box 61548, 05424-970 São Paulo, Brazil Email: [email protected]r, web page: http://www.lmc.ep.usp.br/people/pauletti †IPT - Institute for Technological Research P.O. Box 05508-901 São Paulo, Brazil Email: [email protected], web page: www.ipt.br Key words: nonlinear analysis, membrane structures, Argyris’ membrane element, membrane wrinkling and slackening. Summary. This paper presents the implementation of a simple wrinkling/slackening model into the classical Argyris membrane element, comparing the solution performance by Newton’s iterations, using either the tangent stiffness matrix (numerically evaluated through a finite-difference approximation), or a secant stiffness matrix (obtained through the modification of the elasticity matrix, according to a projection technique which decompose deformations into elastic and wrinkle components). 1 INTRODUCTION There may be cases in which wrinkles on a membrane must be represented with accuracy. In these cases, a shell-like formulation is generally required, and numerical solutions require heavy computation. In the case of architectonical membranes, however, wrinkling is to be avoided, at least on the initial configuration. To guarantee a taut surface, shape finding processes seek configurations where the minimum membrane’s principal stresses are positive everywhere. In past times, this requirement was also extended to the membrane under design loads, and the onset of wrinkling was considered as a type of structural failure. Thus, in order to avoid wrinkling, high initial stresses were usually prescribed at the initial membrane configuration. Nevertheless, in practical architectonical applications, there is not an unconditional reason for the membrane to be free of wrinkles, or even slack zones, in extreme load cases, nor it is required to know the precise pattern of wrinkling, when it occurs. What does is necessary is a method capable to determine the correct load transfer mechanisms, which are distorted if the adopted finite elements do not avoid spurious compression states. Besides, if wrinkling and slackening are allowed in some conditions, the resulting anchorage loads are reduced, because smaller initial membrane stresses are then required, and also because larger membrane displacements contribute to a more favorable load distribution. Implementation of a simple wrinkling model into Argyris’ membrane finite element 143 Ruy M.O. Pauletti, Daniel M. Guirardi. 2 In this paper, we present the implementation of a simple wrinkling/slackening model for the classical Argyris membrane element, already available in the SATS program1,2, and compare the performance of the tangent stiffness matrix (consistent with the proposed wrinkling/slackening model, but numerically evaluated through a finite-difference approximation), to the performance of a secant stiffness matrix (obtained when the classical linear elasticity matrix is replaced by the modified elasticity matrix derived by Akita et al.3, using a projection technique to decompose deformations into elastic and wrinkle components). 2 NONLINEAR EQUILIBRIUM. VECTOR NOTATION Upon discretization, the problem of equilibrium of a membrane structure can be expressed as finding a displacement vector * u such that       * ** gu pu fu 0 (1) where r  uxx is the global displacement vector with respect to a reference configuration r x ,   gu is a residual force vector,   pu is the internal load vector and   fu is the external load vector, all of order dof n , the number of degrees of freedom of the system. These functions can be computed in sub-domains e  , or elements (Figure 1(a)), as a function of the element displacement vectors ee    uu , 1, , e n n   , 1, , e en , where e n is the number of elements and e n n is the number of nodes defining the eth element. Defined onto every element, there exist an element residual force vector ee    gg , where   e ee  g gu is the contribution of the eth element to the residual vector, evaluated at its th  node. The equilibrium problem (1) can be solved –within a vicinity of a solution * u – iterating Newton’s recurrence formula,       1 -1 1 k k k kk k k          u g u u gu u K gu u (2) where we define the tangent stiffness matrix , , 1, , i dof j gij n u           g Ku (3) Recurrence (2) may converge to a solution even when the consistent linearization of (1) is replaced by some approximation (i.e., a secant stiffness matrix), usually at the price of reduced convergence rate. 144 Ruy M.O. Pauletti, Daniel M. Guirardi. 3 The element displacement vectors e u can be extracted from the global one according to ee u Au , where e A is a Boolean incidence matrix for that element. Likewise, the global residual force vector and the global tangent stiffness matrix can be assembled according to 1 e n eT e e  g Ag and , , 1, , e e ee i dof ee j gij n u           g ku (4), (5) where e dof n denotes the number of degrees of freedom of the eth element. 3 INTERNAL LOAD VECTOR FOR A MEMBRANE ELEMENT Figure 1 shows the Argyris’ natural triangular membrane finite element, defined in an initial configuration 0  , in which it is already under an initial stress field. A reference configuration r  usually considers stress-free conditions. For small strains, we assume 0r   . The element’s current configuration is denoted by c  . Element nodes and edges are numbered anticlockwise, with edges facing nodes of same number. Nodal coordinates are referred to a global Cartesian system, and a local coordinate system, indicated by an upper hat, is adapted to every element configuration, such that the ˆ x axis is always aligned with edge 3, oriented from node 1 to node 2, whilst the ˆ z axis is normal to the element plane. e  o xx ˆ , o yy ˆ , o zz ˆ , r  0 ˆ x 0 ˆ z y ˆ x ˆ z ˆ c  0 ˆ y 0  1 l 3 l 2 l 0 1 l 0 3 l 0 2 l r l 1 r l 3 r l2 Configuração Corrente - c Configuração Inicial - 0 Configuração de 2 1 2 1 3 1 3 2 3 12 P 3 r  2 0 ˆ x 0 ˆ z 3 y ˆ x ˆ z ˆ 2  0 P x 3 1  P x t  P  P u 1 0 P 0 ˆ y 0  r yy ˆ , r xx ˆ , r zz ˆ , Initial configuration Reference configuration Current configuration 0 P x P x P u o xx ˆ , o yy ˆ , o zz ˆ , r  0 ˆ x 0 ˆ z y ˆ x ˆ z ˆ c  0 ˆ y0  1 l 3 l 2 l 0 1 l 0 3 l 0 2 l r l 1 r l 3 r l 2 Configuração Corrente - c Configuração Inicial - 0 Configuração de 2 1 2 1 3 1 3 2 3 12 P 3 r  2 0 ˆ x 0 ˆ z 3 y ˆ x ˆ z ˆ 2  0 P x 3 1  P x t  P  P u 1 0 P 0 ˆ y 0  r yy ˆ , r xx ˆ , r zz ˆ , Initial configuration Reference configuration Current configuration 0 P x P x P u Fig. 1: (a) a domain  , discretized into elements e  ; (b) a triangular element in three different configurations; (c) position vector P x and displacement vector P u ; (d) internal angles ,,  and unit vectors i v along the edges; (e) internal nodal forces i p decomposed into natural forces ii Nv . 145 Ruy M.O. Pauletti, Daniel M. Guirardi. 4 Fig. 1(c) displays the current global coordinates of e P . Nodal coordinates are given by 01, 2, 3, iii i xxu , where i u are the element’s nodal displacements. Fig. 1(d) displays the lengths of element edges, given by ii kj   l xx , with indexes , , 1, 2, 3i jk in cyclic permutation. Unit vectors parallel to the element edges are denoted by i ii v ll . With these definitions, the vector of internal nodal forces can be decomposed into forces parallel to the element edges, according to 1 22 33 2 3 1 2 33 11 1 3 2 3 11 2 2 1 2 3 e NN N NN N NN N                         p v v 0vv p p v v v 0 v CN p v v vv0 (6) where C is a geometric operator, which collects the unit vectors parallel to the element edges and   123 T NNNN is the vector of natural forces (see Fig. 1(d)). We assume that the taut behavior of the element is linear-elastic, thus a linear relationship exists, such that 0 r n  N ka N , (7) where   123 T a , is the vector of natural displacements (with 0 1, 2, 3 , i ii i  ) and the element natural stiffness is a constant matrix given by -1 1 -1 ˆ rT n rr r r r V  k TTL LD , (8) where r V is the element volume,   diag r ri L , ˆ D collects the coefficients of Hooke’s law for plane stresses, such that ˆ ˆˆ σε=D and, finally, r T is a transformation matrix, relating the linear Green strains ˆ ε to the natural strains , i.e., ˆ rn Tε ε , highlighting the fact that Argyris’ natural membrane element is akin to a strain gauge rosette. Explicitly we have 22 22 cos sin sin cos cos sin sin cos 10 0 r r rr r r r rr             T , (9) where the internal angles r  and r  are depicted in Figure 1(d). The upper or lower index 'r' indicates that computations are performed in the reference configuration, which, for small strains kinematics, can be superimposed to the initial one. A more detailed deduction of r n k is provided in references1,2. Since r n k has only six independent components, its storage is usually economic, reducing the number of operations required to calculate the internal loads and tangent stiffness, and thus the overall computing time. Inserting (7) into (6), the vector of internal forces at each configuration is given by   0 er n  p Cka N . (10) 146 Ruy M.O. Pauletti, Daniel M. Guirardi. 5 It is interesting to define also an external wind load vector, according to   333 3 T ee w pA  f IIIn , (11) where p is a normal wind pressure acting on the element, A is its area and e n its normal unit vector, in the current configuration. Thus, the element error vector becomes e ee w  gpf (12) Proceeding with derivation of (5), the consistent tangent stiffness matrix of the membrane element is obtained:       33 22 23 3 2 33 11 33 2 1 12 22 23 3 2 123 3 1 3 1 123 11 123 2 1 12 11 6 , NN NN NN NN e rT n N N NN e c g ext p                       MM M M ΛΛΛ k Ck C M M M M ΛΛΛ ΛΛΛ M M MM kkkk (13) where 3 T i ii  M I vv , 1, 2, 3i , and skew( ) ii Λl are skew-symmetric matrices, whose axial vectors are given by i ii lv , and where each component define, respectively, the constitutive, the geometric and the external components of the element tangent stiffness matrix. The contributions of e g and e k to (4) and (5) are given according to 1 2 33 ee e i jk   AA A I and 123 ee e mmm   AAAO ,     1, 2, , \ , , n m n i jk . 4 A SIMPLE WRINKLING MODEL Matrix r n k is constant, and provides a fast way to compute the element’s internal loads and tangent stiffness, when the membrane is fully under tension, avoiding explicitly calculation of strains or stresses during solution, as they can be post-processed after equilibrium is achieved. However, before developing compressive stresses, membranes become wrinkled or slack. Criteria for identification of the status of a membrane element (taut, wrinkled or slack) require consideration of stresses, strains, or both. Therefore, when wrinkling or slackening are possible, equations (10) and (13) have to be replaced for lengthier calculations. For isotropic materials, the principal stress and strain directions are parallel, and stress, strain or mixed wrinkling criteria are equivalent [7]. We choose a pure stress criterion and calculate element stresses directly from natural displacements. In order to further speed up calculations, we first decompose the element natural stiffness r n k according to     -1 1 -1 ˆ r T rr n rr r r r a V    k T T kkL LD , (14) 147 Ruy M.O. Pauletti, Daniel M. Guirardi. 6 where we observe that r  k and   1ˆT rr ar V   kkD are two symmetric constant matrices that can be conveniently stored during the pre-processing phase. Then stresses in any configuration can be evaluated according to 0 ˆˆ r a σ ka σ= , (15) and, after calculating the principal stresses 1,2  and the “principal angle” 1  (between axis ˆ x and the 1  direction), we determine the element status and eventually modify stresses according to the following criterion:       2 1 12 1 1 1 1 ˆˆ 0 TAUT ˆ 0 0 WRINKLED 1 cos 1 cos sin 2 2 ˆ 0 SLACK T                          σσ σ σ0 (16) Thereafter, we replace (10) by (17): ˆ er  p Ck σ . (17) A mixed stress-strain criterion could avoid the intermediate calculation of negative stresses, which have no physical meaning for membranes, but by using equations (15) and (17) we determine stresses and internal loads without explicit calculation of strains, thus a stress criterion becomes more effective in our framework. 5 A FINITE-DIFFERENCE ESTIMATIVE OF THE CONSISTENT STIFFNESS MATRIX Equation (17), inserted into (12) is all what is needed to solve the equilibrium problem through dynamic relaxation, as done for instance in reference4. In order to use Newton’s iterations, however, also the system’s tangent stiffness matrix is required. Nevertheless, instead of performing the consistent linearization of (17), we opt to evaluate it numerically, through the finite-difference approximation developed in references5,6, which was shown to yield the same convergence rate and precision as the consistent tangent stiffness matrix, for highly nonlinear problems of cables, membrane and shell structures, with fairly acceptable extra computational cost. Such type of finite-difference procedures is commonly used in nonlinear mechanics, to compute approximate consistent tangent moduli, for complicated material laws   σ σε . In our approach, however, instead of approximating the tangent modulus, we directly approximate the tangent stiffness matrix, including all nonlinear effects that might affect the global error vector g . First, we first partition the element tangent stiffness matrix (13) according to 12 e dof e ee e n    k kk k (18) 148 Ruy M.O. Pauletti, Daniel M. Guirardi. 7 where column j j x   g k can be interpreted as the directional derivative of g with respect to the jth component of x . Thereafter, we approximate these e dof n directional derivatives e j k by a central-difference scheme, according to     1, 1, , 2 e ee e e e e j j j dof h hj n h     k gx δ gx δ , (19) where h is a finite scalar parameter and e j hδ are backward and forward perturbations of the jth element degree of freedom, such that , 1, , ee j ij dof in      δ , and where ij  is the Kronecker delta. Finally, inserting approximations (19) into (18), we obtain numerical estimates for the tangent stiffness matrices e k , which are then assembled into a global a stiffness matrix K . In references5,6 we have shown that this method provides excellent approximations for the consistent tangent stiffness matrix, as long as h is small enough. Moreover, since no extra cost is associated to reducing the size of h , it can be taken as a function of the machine precision  , as small as possible, but without introducing numerical noise. In MATLAB environment, where SATS was implemented, we found that 4 h   is a good compromise estimative. 6 A SECANT STIFFNESS MATRIX We have also investigated the use of a secant approximation to the stiffness matrix, by which instead of performing the consistent linearization of equation (17), we simply modify the elements’ natural stiffness, equation (14), according to a modified elasticity matrix, borrowed from the paper by Akita et al.3. After decomposing the total strains on a membrane element into elastic and wrinkled fractions, according to ˆˆ ˆ ew  εε ε , these authors finally arrive to 22 22 ˆ ˆˆ ˆ T wT   ssD ε Qε s Ds , (20) where Q is a projection matrix, that extracts the wrinkled portion of the deformation from the total one. Vectors 1 s and 2 s are such that 11 2 2   εs s , with 1  and 2  being the element’s principal strains (never actually calculated, in our method). Explicitly,         1 1 11 2 1 11 11 cos2 1 cos2 2sin 2 2 11 cos2 1 cos 2 2sin 2 2 T T                 s s (21) Now since ˆw ε does not rise stresses,     ˆˆ ˆ ˆ ˆˆ ˆ ew     σ Dε D ε ε D I Q ε Dε , where   ˆˆ  D I QD (22) is the modified elasticity matrix we seek. Inserting (22) into (14) and that into (13), we finally obtain a modified, secant stiffness matrix 149 Ruy M.O. Pauletti, Daniel M. Guirardi. 8 e c g ext kkkk , (23) in which only the constitutive part c k is approximate. In principle, this non-consistent stiffness matrix may slow down convergence rates, but its use in reference7 allowed easy solution, as shown in the following benchmark. 7 A FIRST BENCHMARK – THE ‘MEMORIAL DOS POVOS’ In references4,7 we presented results obtained by the SATS program2 for the membrane roof of the “Memorial dos Povos de Belém do Pará”, shown Figure 2, under a uniform upward wind load acting over the whole membrane surface, considering both fully-adherent and frictionless sliding border cables. A discretization much coarser than the one used in actual design was adopted, to ease the visualization of results. In the present paper, we consider only results obtained for the fully-adherent model, which can be compared also with results given by the Ansys FEM code (a sliding cable is not directly available in Ansys, thus the sliding condition was not analyzed by that program). Figure 2: The membrane roof of the “Memorial dos Povos de Belém do Pará” We also compare results obtained by SATS, through Newton’s iterations, using both the tangent or the secant stiffness matrices ( e k or e k ), with those obtained through dynamic relaxation4, a method which only requires the definition of the modified error vector g , yielding an independent checking for results, besides Ansys model. Table 1: Comparison of some selected results, for   5 0 / 10  gg Program SATS Ansys Method Newton ( tangent K ) Newton ( secant K ) Dynamic Relaxation Newton + linesearch Maximum displacement [m] 0.40403 0.40403 0.40343 0.40285 Maximum 1  [MPa] 13.60979 13.60974 13.59338 13.6 Minimum 1  [MPa] 5.47444 5.47442 5.47499 5.47 Maximum 2  [MPa] 5.37598 5.37598 5.37532 5.37 Minimum 2  [MPa] 0.00000 0.00000 0.00000 0.00014 Number of iterations 5 7 -- 12 Time to solution [s] 4.586 4.353 -- 1.747 Time per iteration [s] 0.9172 0.6218 -- 0.1456 150 Ruy M.O. Pauletti, Daniel M. Guirardi. 9 1 MN MX 0 .044756 .089512 .134268 .179024 .223781 .268537 .313293 .358049 .402805 APR 1 2009 09:50:53 NODAL SOLUTION STEP=3 SUB =6 TIME=3 USUM (AVG) RSYS=0 DMX =.402805 SMX =.402805 1 MNMX .547E+07 .637E+07 .727E+07 .817E+07 .907E+07 .997E+07 .109E+08 .118E+08 .127E+08 .136E+08 APR 1 2009 09:51:54 ELEMENT SOLUTION STEP=3 SUB =6 TIME=3 S1 (NOAVG) DMX =.402805 SMN =.547E+07 SMX =.136E+08 1 MN MX 800 597378 .119E+07 .179E+07 .239E+07 .298E+07 .358E+07 .418E+07 .477E+07 .537E+07 APR 1 2009 09:57:32 ELEMENT SOLUTION STEP=3 SUB =6 TIME=3 S2 (NOAVG) DMX =.402805 SMN =141.967 SMX =.537E+07 Figure 3: Results by SATS (both Newton’s iterations dynamic relaxation) compared to results by Ansys; Top to bottom: displacement norms; 1  on elements; 2  on elements, wrinkled elements shown in grey. Converged results obtained considering these four different methods were all in good agreement, as can be seen in Table 1, which compares some selected results obtained with SATS and Ansys models. In the Ansys model, the wrinkled elements detected in SATS presented very low but still positive 2nd principal stresses, possibly due to the post-processing extrapolations adopted by that program. Figure 3 compares displacement’s norms and principal stresses on elements, obtained using SATS (differences of results through different methods are visually imperceptible, so figures are not repeated) and Ansys. In the plotting generated by SATS, principal directions are shown with lines of length proportional to the