A finite element interior-point implementation of tension field theory

Computers and Structures 151 (2015) 30–41

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

A finite element interior-point implementation of tension field theory R. de Rooij ⇑, M.M. Abdalla Aerospace Structures & Computational Mechanics, Delft University of Technology, Kluyverweg 1, 2629 HS, The Netherlands

a r t i c l e

i n f o

Article history: Received 11 December 2013 Accepted 7 January 2015

Keywords: Membranes Relaxed strain energy function Interior-point method Wrinkling Variational formulation Finite element method

a b s t r a c t This paper presents an interior-point implementation method of the relaxed strain energy function to model the stress distribution in membranes in the presence of wrinkles. The relaxed strain energy function is reformulated using the interior-point method. This formulation is proven to be quasi-convex with respect to the deformation gradient. A set of governing equations for the membrane is developed in the standard finite element format. The solution to this set of equations is obtained using the interior-point method, and several numerical examples are presented for validation. The proposed method is shown to converge and to be robust and efficient. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction Membrane structures have a rich history and are widely used in aerospace and structural engineering applications. A few examples can be found in solar sails, airbags, atmospheric balloons and parachutes. These membranes have many advantages, including their ability to take complex shapes and their low mass to surface ratio, which is especially important in aerospace engineering. Although membranes can carry tensile loads very well, they tend to wrinkle under the slightest compressive load. It is important to be able to predict the stress distribution in the membrane, taking into account the possible presence of wrinkles, as this affects the load carrying capabilities of a membrane. Wrinkles in a structure are often of a much smaller scale than the structure itself. In order to model these wrinkles the mesh size in e.g. a finite element model needs to be smaller than the wrinkle size, leading to high computational costs. A possible solution is to model the wrinkles as a continuous in-plane contraction of the membrane. The exact geometry of the wrinkles is lost by this procedure, however the stress field in the membrane is modeled accurately, and the computational costs are reduced significantly. This solution, known as the tension field theory, is widely used and was first introduced by Wagner [1]. A few geometrically exact models for membrane wrinkling exist in which the tension field theory method is not applied, yet the ⇑ Corresponding author. Tel.: +31 15 27 87308; fax: +31 15 27 85337. E-mail addresses: [email protected] (R. de Rooij), [email protected] (M.M. Abdalla). URL: http://ae.tudelft.nl (R. de Rooij). http://dx.doi.org/10.1016/j.compstruc.2015.01.007 0045-7949/Ó 2015 Elsevier Ltd. All rights reserved.

information about the detailed wrinkling configuration is maintained. Puntel [2] gives an analytical solution for stretching a flat membrane and an analysis of a point load on a pressurized sphere is given in Vella [3]. Geometrically exact numerical models are presented in Flores et al. [4] and Weinberg et al. [5]. Capturing wrinkling geometry requires the mesh size in a finite element model to be smaller than the wrinkling in the membrane which results in high computational costs for very thin membranes [6]. Additionally the detailed wrinkling geometry of the membrane often is of small importance, therefore the tension field theory method is used in most models for membranes. The governing equations for the load carrying capabilities of a membrane are presented in the tension field theory method [1]. It is assumed that the membrane has zero bending stiffness and all compressive stresses are eliminated. This elimination is enabled by the introduction of an in-plane contraction which carries zero strain energy. Analytical solutions to these equations exist for a few specific cases. Mansfield [7,8] presents solutions for the tension field of a rectangular strip under shear and an annular ring under torsion loading. These solutions have been validated with experiments and serve well for verification purposes. Pipkin [9,10] in his work incorporates the tension field theory into the theory of finite deformations as a special case of non-linear elasticity. A modified stored energy function is developed by relaxing the strain energy function of the material to obtain a nondecreasing and convex relaxed strain energy function. Based on this work Steigmann [11] presents the tension field theory in terms of two potential functions. Different membrane states are derived and it is demonstrated that application of the relaxed strain energy leads to stable solutions for the stress

31

R. de Rooij, M.M. Abdalla / Computers and Structures 151 (2015) 30–41

distribution. The first numerical results based on the relaxed strain energy function are also presented by Steigmann [12]. Several numerical methods in which the tension field theory method is used have been developed over the years. In one class of methods the projection technique of Akita et al. [13] is implemented. Here the deformation of the membrane is split into a strain carrying and zero-strain carrying part. This is realized by projecting the zero-strain energy components out of the stiffness matrix, which modifies the constitutive law. Jarasjarungkiat et al. [14] have extended this method for use with orthotropic materials. In a second class of numerical methods for modeling the tension field in membranes the constitutive law is maintained, but the deformation gradient is modified. Roddeman et al. [15] present the theory and analysis for this modification in isotropic materials while Schoop et al. [16] implement it using finite elements. Raible et al. [17] and Epstein [18] extend the method to enable the modeling of orthotropic and general anisotropic materials respectively. Pagitz et al. [6] propose a new finite element that includes a rotational degree of freedom to control the direction of the tension in the membrane. The performance of this element is verified even for the membrane regions close to slacking. Dynamics methods to model tension fields are developed in Kang [19] and Miyazaka [20]. An extension into non-linear material laws to account for plasticity is presented by Mosler et al. [21] using a variational formulation of the internal strain energy. The tension field theory method is formulated using the relaxed strain energy function by Pipkin [10], and an explicit parameterization of the wrinkling strain tensor is applied. The convexity properties of the energy functional are not preserved by this parameterization. In this paper a new implementation of the relaxed strain energy function by Pipkin [10] will be presented. Convexity properties of the relaxation are preserved and exploited, and the final set of governing equations is presented in the standard finite element format. For this purpose, a reformulation of the relaxed strain energy function will be derived for which the convexity properties yield a computationally efficient method. The resulting reformulation will be a constrained, convex optimization problem, which will be solved using the interior-point method. With this implementation no if-else statements are required in the proposed method, as they are in many of available methods to distinguish between the three possible states of a membrane: taut, wrinkled or slack. The outline of the paper is as follows: the kinematics of membranes is presented in Section 2, which is followed by the general equilibrium conditions in Section 3. These general equilibrium conditions lead to an unstable energy minimization, as shown in Pipkin [10]. Therefore, a variational principle based on the relaxed strain energy function is applied in Section 4 to obtain a stable set of governing equations. The discretization of the membrane structure and the linearization of the equations is performed in Section 5. The numerical implementation of the discretized equations is discussed in Section 6. Several numerical results are presented in Section 7 in which the proposed method is verified and validated. The conclusions are presented in Section 8. In addition, two appendices have been attached for completion. The mathematical conditions for quasi-convexity are presented in Appendix A. These conditions are related to the relaxation of a function in Appendix B, and a closed form expression of this relaxation is derived.

wrinkles are modeled as an in-plane contraction of the membrane, meaning that all strains are membranal. As a result only in-plane strains need to be considered in this section. Consider the membrane in the reference/undeformed configuration occupies a bounded plane with covariant basis vectors Aa ða ¼ 1; 2Þ. Note: Greek characters are used to denote the reference configuration. The reference position of any material point can then be given as r 0 ¼ Aa na , with na being the contravariant components of the material point. In the current/deformed configuration the membrane can occupy the space normal to the reference configuration which is described by the covariant basis vectors ai ði ¼ 1; 2; 3Þ. The current position vector then becomes r ¼ ai xi . Let Gab and g ij be the metric tensors in the reference and current configurations respectively:

Gab ¼ Aa  Ab

ð1Þ

g ij ¼ ai  aj

ð2Þ

The deformation from the reference to the current configuration is described by the deformation gradient F ¼ F ia Aa ai :

F ia ¼

@xi @na

ð3Þ

This shows that the components of F form a 3  2 matrix, which means F has a null space along the normal of the membrane. The null space is spanned by the vectors u for which the pull back by the deformation gradient equals zero: F T u ¼ 0. The Green–Lagrange strain is a symmetric 2  2 tensor defined in the reference configuration:

E¼

 1 T F F I 2

ð4Þ

In convected coordinates a material point in the current configuration is described by the same set of coordinates as in the reference configuration. The Green–Lagrange strain tensor in convected coordinates is given as:

E¼

1 ðg  G Þ 2

ð5Þ

Let u be the displacement vector of this material point, then:

ai ¼ Ai þ

@u ¼ Ai þ u;i ; @ni

i ¼ 1; 2

ð6Þ

The engineering strain in Voigt form then becomes:

0

1

A1  u;1

e¼B @

A2  u;2 A1  u;2 þ A2  u;1

0

C 1B Aþ @ 2

u;1  u;1 u;2  u;2

1 C A

ð7Þ

u;1  u;2 þ u;2  u;1

3. Equilibrium equations Let /ðEÞ ¼ /ðEðFÞÞ ¼ wðFÞ be the strain energy function of the material, which is assumed to be convex in E as is the case when the stress–stain relation is one-to-one [22,10]. The first and second Piola–Kirchhoff stress tensors are given in (8) and (9) respectively: i

P ¼ Pai aa b

Pai ¼

S ¼ Sab aa ab

Sa b

@w

@F ia @/ ¼ @Eab

ð8Þ ð9Þ

2. Kinematics

Because of symmetry in E, also S is symmetric. The relation between P and S is given as:

Membranes are very thin shells and thus they have a low bending stiffness, which is neglected in the analysis of the tension field in membranes; this causes the membrane to wrinkle under the slightest compressive deformation. In tension field theory these

Pai ¼ Gij F jc Sca ) P ¼ FS

ð10Þ

The strong form of the equilibrium equations for any material point in the membrane is given as:

32

R. de Rooij, M.M. Abdalla / Computers and Structures 151 (2015) 30–41

Pai;a þ bi ¼ 0 ) r  P þ b ¼ 0

ð11Þ

with b being the external force. Alternatively a weak form of this equilibrium equation can be formulated, which in this case is the minimization of the total potential energy P in the membrane:

PðuÞ ¼

Z

/ðEðuÞÞdX 

Z

X

X

b  udX 

Z

t  ud@ X

ð12Þ

@X

in which u is the displacement field and b and t represent body forces and boundary tractions respectively. The displacement field minimizing the potential energy is a stable state of the membrane:

Þ u ¼ arg inf Pðu  u

A unique solution for u can only exist when boundary conditions are prescribed:

P  n ¼ tp u ¼ up

at @ Xr at @ Xu

ð13Þ ð14Þ

in which @ Xr and @ Xu represent the parts of the boundary where tractions t p or displacements up are prescribed and n is the normal vector to the boundary. The existence and uniqueness of an energy minimizing displacement field requires quasi-convexity of wðFÞ, which is weaker than convexity and stronger than rank-one convexity [23,24]. Pipkin [10] shows that when compressive stresses, S < 0, develop in a membrane this requirement is violated and the energy minimization becomes unstable. This results from the zero bending stiffness of the membranes causing immediate wrinkling upon compression. The physical reason for this to happen is the negligible bending stiffness of the membrane, meaning that any compressive stresses immediately lead to buckling/wrinkling of the membrane. Consequently the only load carrying capabilities of the membrane are in tension. Three distinct membrane states can be distinguished: taut, wrinkled and slack. In the taut state the membrane is under pure tension and behaves as a regular shell, the slack state represents a pure compression in which the membrane has no load carrying capabilities. In the wrinkled state the membrane is tense in one direction and slack in the perpendicular direction and can thus only carry loads along the tensile direction. The wrinkling deformation in the membranes is modeled as an in-plane contraction as described by the tension field theory method. To enable a stable energy minimization for a membrane the relaxed strain energy is introduced by Pipkin [9,10], which leads to the required quasi-convexity of wðFÞ. Wrinkling of the membrane is taken into account using an in-plane contraction with zero strain energy. The minimum of the relaxed energy function is equal to the minimum of wðFÞ given that /ðEÞ is convex as assumed. 4. Variational principle The tension field theory method can be formulated with a relaxed strain energy function, which leads to a stable energy minimization as shown in Pipkin [10]. Based on Pipkin’s approach, we derive a closed form expression for the relaxation of a function in a slightly modified form in Appendix B. The relaxed strain energy function is reformulated in this section using the interior-point method. The proposed variational principle is then applied to minimize the total potential energy in the membrane. In addition to the Green–Lagrange strain tensor E, the wrinkling ^ is introduced to formulate the relaxed strain energy strain tensor E ^ represents the negative of in-plane contraction caused function. E by wrinkles and carries zero strain energy as proposed by the ten^ > 0, resulting in an elastic sion field theory. A wrinkle causes E

^ A third strain tensor is defined as the slack strain equal to E þ E. ~ ^ strain tensor E, this behaves like the wrinkling strain tensor E. The second Piola–Kirchhoff stress tensor is computed as:

S¼

 @/ @E EþE^

ð15Þ

The actual stress tensor in the membrane is denoted by the var^ iable S. 4.1. Relaxed strain energy function The quasi-convexification of wðFÞ is the largest quasi-convex function that does not exceed wðFÞ ¼ /ðEÞ. In this paper, the strain energy function /ðEÞ is assumed to be convex. It is shown in Proposition 3 of Appendix A that the quasi-convexification is then obtained by the relaxed strain energy function /] ðEÞ as given in (4):

h i ^  S^ : E ^ /] ðEÞ ¼ sup inf /ðE þ EÞ ^ E

^ SP0

h  i ^  S^ : E ~ ^E ¼ sup inf /ðE þ EÞ

ð16Þ

^ EP0 ~ E;

S^

~ is introduced in (16) to avoid the explicit conThe slack strain E ^ straint S P 0. This expression for the relaxed strain energy function will lead to a stable energy minimization due to its quasiconvexity in the deformation gradient F. The only explicit constraint in the formulation is the semi-def~ P 0. We propose to use an interior point method inite constraint E [23,25,26] to eliminate this constraint. A barrier function which equals 0 when the constraint is satisfied and 1 elsewhere, is added to the objective function to represent the effect of the constraint. A convex function approaching these properties is given by j:

8   <  lim l ln E ~  if E ~>0 ~ ¼ l!0þ jðEÞ : 1 elsewhere

ð17Þ

~ > 0, the in which l is called the penalty parameter. Assuming E relaxed strain energy function then becomes as in (18).

^ lÞ ^ E; ~ S; /] ðEÞ ¼ lim sup inf vðE; E; l!0þ

S^

ð18Þ

^E ~ E;

with:



 



^ lÞ ¼ /ðE þ EÞ ^ E; ~ S; ^  S^ : E ~  l ln E ~  ^E vðE; E;

ð19Þ

4.2. Potential energy minimization As shown in Section B the relaxed strain energy function /] ðEÞ is quasi-convex in the deformation gradient F and thus in the displacement u. The potential energy can be written with this relaxed strain energy to become:

PðuÞ ¼

Z

/] ðEðuÞÞdX 

X

Z

b  udX 

Z

X

t  ud@ X

ð20Þ

@X

With the definition of the relaxation in (18) and (19), the potential energy functional becomes:

^ lÞ ¼ ^ E; ~ S;u; ^ ðE; P

Z

^ lÞdX  ^ E; ~ S; vðE; E;

X

Z X

b  udX 

Z

t  ud@ X

ð21Þ

@X

and the minimum potential energy is then defined as:

^ u; lÞ ^ E; ~ S; ^ ðE; lim sup inf P

^ E;u ~ l!0þ S^ E;

ð22Þ

The equilibrium configuration of the membrane is represented by the values of the variables for which the sup and inf operators

33

R. de Rooij, M.M. Abdalla / Computers and Structures 151 (2015) 30–41

are obtained. These operators are satisfied if the first order variation of the potential energy equals zeros, which gives:

Z  Z     ^ X ^ Xþ ~ 1 : dEd ~ Xþ ^ E ~ : dSd E S  S^ : dEd S^  lE X X Z Z X Z @/ þ  dudX  b  dudX  t  dud@ X ¼ 0 ð23Þ X @u X @X

^¼ dP

Z 

5.2. Linearization The linearization of the governing Eqs. (24) and (25) requires the introduction of the material, wrinkling and effective stiffness tensors as shown in Table 1. The wrinkling ‘stiffness’ tensor H relates an increment in compressive stress to an increment in wrinkling strain such that ^

The underlined part of (23) is the variational form of standard potential energy function in (12) and it represents the weak form of equilibrium. By applying the localization theorem to the first ~ the governing equations three integrals in (23) and eliminating E, ^ are then obtained: for S^ and E

S  S^ ¼ 0 ^ S^ ¼ lI E

ð24Þ ð25Þ

Eq. (24) states that the actual stress S^ should be as given by the constitutive law of the material whereas the three states of the membrane are captured by (25). In the limit of l ! 0þ , this equation is ^ ¼ 0 in taut regions, S^ ¼ 0 in slack regions or when satisfied by E ^ and S^ are rank-one with perpendicular eigenvectors in wrinboth E kled regions. After applying consistent boundary conditions, Eqs. (24) and (25), and the underlined part of (23) uniquely describe the displacements, wrinkling strains and stresses in the membrane. The solution to this set of non-linear equations is obtained by discretization and linearization. 5. Discretization and linearization In this section we present the discretization of the membrane structure and linearization of the governing equations required to solve the governing equations for the unknown displacements, wrinkling strains and stresses. First the membrane structure is discretized into finite elements, which are widely used in mechanics as they can provide good approximations for the behavior of complex structures. 5.1. Discretization Discretizing the membrane structure is required to solve the governing equations over the complete membrane. The discretization is done using finite elements after which the initial configuration of the membrane r 0 and displacement u of the elements are represented by the shape functions N:

r0 ðnÞ ¼ N ðnÞq0

ð26Þ

uðnÞ ¼ N ðnÞv

ð27Þ

in which q0 and v are the nodal position and displacement vector of the finite element respectively. Substituting these relations into (7) yields the engineering strains in the membrane:

0

1

0

1

0

1

e1 qT0 P 1 v v T P1 v C B T C 1B T C e e¼B ¼ þ q P v v @ 2 A @ 0 2 A @ P2 v A 2 c12 qT0 P 12 v v T P 12 v

ð28Þ

in which:

H ¼  @@ES^ , as becomes apparent in (33). Table 1 shows that H approaches infinity in taut regions of the membrane where ^ ¼ 0. Physically, this means that incremental compressive stresE ses do not lead to wrinkles in taut regions of the membrane, as long as the regions remain taut. The effective membrane stiffness tensor C ¼ C in the taut regions of the membrane. In slack regions of the membrane, where S^ ¼ 0, the wrinkling ‘stiffness’ tensor H ¼ 0. Physically, it means that any increment in compressive stress would lead to an infinite increase in wrinkling strain. This shows that the membrane becomes unstable if a compressive force is applied. Conversely, it also shows that in case a compressive displacement is applied, the membrane stress does not become negative, i.e. the constraint S^ P 0 is satisfied. Linearization of (24) and (25) around any current state yields:

  ^  DS^ ¼ Rr C DE þ DE

ð30Þ

^ S^ þ E ^ DS^ ¼ Rl DE

ð31Þ

with Rr and Rl being the residuals of the governing equations at the ^ which ^ DS, current state; Rl also contains the non-linear term DE ^ operation. Solving the linearized equations ^ SÞ arises from the DðE for the stress and wrinkling strain updates gives:

DS^ ¼ C ðRw þ DEÞ   ^ ¼ H1 E ^ 1 Rl  DS^ DE

ð32Þ ð33Þ

in which:

^ 1 Rl  C1 Rr Rw ¼ H1 E Eq. (32) represents the effective constitutive law for a membrane. It relates the total membrane strain to the membrane stress through the effective stiffness tensor C , as defined in Table 1. Eq. (33) relates the wrinkling strain to the membrane stress through the wrinkling ‘stiffness’ tensor H, which was introduced in Table 1. Note that the residuals in (32) and (33) approach zero in the converged ^ 1 in (33) is well-defined due to the intesolution. Also note that E rior-point method implementation, see Section 6. This implementation ensures that the wrinkling tensor always contains two strictly positive eigenvalues. The final set of governing equations is obtained by linearizing the weak form of equilibrium, as represented by the underlined part of (23). Before linearizing, first the structure is discretized into finite elements and the integrals are evaluated using Gauss integration over integration points i. Let B ¼ @@Ev and f be the force vector containing the integrations over b and t, then the discretized weak form of equilibrium becomes:

X T Bi S i wi  f ¼ 0

ð34Þ

i

in which wi represents the Gauss weight of the integration point. Linearization of this equation about a current state yields: Table 1 Different fourth-order stiffness tensors.

P1 ¼

N T;1 N ;1

ð29aÞ

P2 ¼

N T;2 N ;2

ð29bÞ

C

Description  @S  Material stiffness tensor = @E ^ EþE

ð29cÞ

H

^¼E ^ 1 DE ^ S^ 8DE ^ Wrinkling ‘stiffness’ tensor s.t. HDE

C

Effective membrane stiffness tensor = C  CðC þ HÞ1 C

P 12 ¼ N T;1 N ;2 þ N T;2 N ;1 T The matrices B ¼ @@ve and @B S^ are derived from (28). @v

Parameter

34

R. de Rooij, M.M. Abdalla / Computers and Structures 151 (2015) 30–41

 X T DBi S i þ BTi DS i wi  Df ¼ Rv

ð35Þ

i

in which Rv is the residual of the governing equation often denoted as f ext  f int and Df is the incremental applied force which usually is non-zero in the first iteration only. Applying (32) and assuming DBTi S i  DBTi S^i , (35) can be rewritten to be:

K e Dv ¼ Df þ Re

ð36Þ

with:

! X @BT i ^ Ke ¼ S i þ BTi Ci Bi wi @v i X T  Re ¼ Rv  Bi Ci Rwi wi

ð37Þ ð38Þ

i

The subscript e is added to specify that this stiffness matrix K e and residual vector Re are defined at element level. Applying (32), (33) and (36) iteratively results in the wrinkling strains and stresses in the integration points and the nodal displacement vector.

These predictor wrinkling updates are used to update the Rl in (31) for the corrector step. In this step the penalty parameter is also reduced as:

lcorr ¼ bc l With bc ¼ max fbaux ; b2 g, in which b2 is a user-defined parameter yielding best convergence for b2 ¼ 0:1 and bc is computed as:

0

baux

  12 ^ ^ þ a e DE S^ þ ar DS^ : E A ¼@ ^ S^ : E

The updated residual vector Rl and lcorr are used to compute the corrected stress and wrinkling strains updates from (32) and (33) with their step sizes. The numerical procedure discussed so far is performed at integration point level within a finite element model. In a full model this procedure is applied for each integration point, except for computing the penalty parameter because l is defined at global level rather than at integration point level. Suppose ne is the number of finite elements and ng the number of integration points per element, then: 1

Pne Png

^ i : S^i Aj wi E

j¼1 i¼1 l¼2 P n e Pn g j¼1

6. Interior-point method ^ þ DE ^ and The updates in (32) and (33) do not guarantee E S^ þ DS^ are positive definite as is required, therefore the semidefinite programming algorithm (SDPA) [26] is applied in which step sizes are applied for both updates. The step size the for wrinkling strain update ae should satisfy:

^ þ ae DE ^ P 0; E

with 0 6 ae 6 1

ð39Þ

Note that the step size cannot exceed 1 since this would violate the linearization of the governing equations. The unconstrained step size ae is obtained by solving (40):

jE^ þ ae DE^ ¼ 0;

with 0 < j < 1

ð40Þ

^ þ DE ^ is strictly in which j is a safety parameter to guarantee E positive. Heuristically it has been found that j  0:9 yields best convergence of the solution. From (40), ae becomes:

1

ae

¼

1

j

  ^ 1 DE ^ eigmax E

ð41Þ

in which eigmax computes the maximum eigenvalue, which is most critical because it determines the smallest step size satisfying (40). Considering the constraint on ae in (39), the step size for wrinkling strains is given as:

ae ¼



ae if 0 < ae < 1

ð43Þ

ð44Þ

i¼1 Aj wi

Algorithm 1. Pseudo-code of numerical implementation. Step j, iteration k 1. Compute K e and Re from (37) and (38) for each finite element using results from iteration k  1. 2. Assemble global stiffness matrix K k and residual force vector Rk . 3. Obtain applied force Df which is non-zero iff k ¼ 1. 4. Predictor step: 4.1. Compute predictor displacement increment Dv pred ¼ k pred ¼ v k1 þ Dv pred . K 1 k ðDf þ Rk Þ and displacement vector v k k 4.2. Compute in each integration point i the total strain Epred i;k pred pred from v pred , and DEpred k i;k ¼ Ei;k  E i;k1 . 4.3. Only if (k ¼ 1 & j ¼ 1), see Algorithm 2. ¼ b1 lk1 . 4.4. Reduce penalty parameter lpred k pred 4.5. Compute residuals Rpred l;i and Rw;i at each integration

, and predict wrinkling updates point i based on lpred k ^ pred and DS^pred from (32) and (33). DE i;k i;k ^pred ar DS^pred . ^ pred ae DE^ pred 4.6. Update step sizes: DE i;k i;k and DS i;k i;k pred ^ ^ i;k1 þ DE ^ pred 4.7. Update predictor wrinkling strains:Ei;k ¼ E i;k and stresses:S^pred ¼ S^i;k1 þ DS^pred . i;k

i;k

5. Corrector step: ^ pred and DS^pred . 5.1. Update Rk based on DE i;k i;k

ð42Þ

5.2. Compute corrector displacement increment Dv k and consequently v k .

In the same way the step size for the stress update ar is computed. From Section 4.1 it can be seen that the penalty parameter l should approach zero in the limit. In the numerical implementation this is realized by decreasing its value in each iteration using the predictor–corrector approach presented in Mehrotra [27]. Sup^ and pose the current wrinkling strains and stresses are given by E ^ and the current penalty parameter is l. In computing DE ^ and DS^ S, one aims to reduce this penalty parameter, which is achieved by defining lpred to be:

5.3. Compute total strain and strain increment in each integration point i : Ei;k and DEi;k .

1

otherwise

lpred ¼ b1 l in which b1 is a ‘centering parameter’ which results in best convergence of the solution for b1 ¼ 0 [26]. Substitution of lpred in (32) and ^ pred and (33) gives the predictor wrinkling updates DS^pred and DE their corresponding step sizes can be computed.

5.4. Reduce penalty parameter 5.5. Compute, based on l ^ i;k and DS^i;k . DE

corr k ,

lcorr ¼ bc lk1 . k

the corrector wrinkling updates

ae DE^ i;k and DS^i;k ar DS^i;k . ^ ^ i;k1 þ DE ^ i;k and stresses: 5.7. Update wrinkling strains:Ei;k ¼ E S^i;k ¼ S^i;k1 þ DS^i;k . 6. Compute internal force vector f int;k . 7. Compute penalty parameter lk as in (44). 8. Check convergence criteria on l and strain energy / as: lk < l ¼ 106 and kD/k =D/1 k < / ¼ 104 . 9. Proceed to next iteration if not converged, proceed to next step/finish if converged. end ^ i;k 5.6. Update step sizes:DE

35

R. de Rooij, M.M. Abdalla / Computers and Structures 151 (2015) 30–41

In the finite element model the load is applied in several steps which all consist of iterations to reach convergence of the step results. Algorithm 1 shows the pseudo-code for step j and iteration k. The procedure in this algorithm is complete, except for the first iteration in each step. As is the case in any Newton-iteration method an initialization must be specified to start the optimization. ^ and For the membrane this means that the wrinkling strains E stresses S^ have to be initialized in each integration point. Theoretically the choice of initialization is arbitrary given that S^ > 0. However, to obtain an efficient and numerically robust method the initialization is based on the total strain E. Algorithm 2 shows the modifications in the pseudo-code for the first iteration of a step. Algorithm 2. Modifications in numerical implementation for first iteration. Fig. 1. Sketch of rectangular membrane strip under displacement loading.

if k ¼ 1 and j ¼ 1 then . First iteration in first step 1. Compute K e and Re from (37) and (38) assuming for each ^ i;k1 ¼ 0 and Rw ¼ 0. integration point S^i;k1 ¼ E i 4. Predictor step: ^ i;k1 in each integration point i from Ei 4.3.1. Initialize S^i;k1 ; E computed in item 4.2 in Algorithm 1. 4.3.2. Initialize lk1 from (44). else if k ¼ 1 then . First iteration in any other step 0. Initialization: ^ i;k1 in each integration point i from con0.1. Initialize S^i;k1 ; E verged Ei at the previous step j  1. 0.2. Initialize lk1 from (44). end if

The initialization of wrinkling strains and stresses is performed on integration point level and is based on the total strain E. Let S^0 be the stress tensor computed from the strain energy function of the material without any wrinkling effects:

 @/ S^0 ¼  @E E

which can be negative depending on the total strain. Let the initialization of the stress S^ be computed as:

S^ ¼ S^0 þ pI

ð45Þ

such that S^ > 0. Any value for p satisfying this constraint is allowed, however to increase efficiency and numerical robustness its value is obtained from:

p ¼ maxf0; k1  kp ðk1  k2 Þg

ð46Þ

in which k1 and k2 are the principal values of S^0 with k1 6 k2 . The value of kp in (46) indicates the margin of satisfying S^ > 0 and yields ^ can be comthe best convergence for kp  3. The initialization of E puted from (25) to be:

^ ¼ l S^1 E i

ð47Þ

in which li is the penalty parameter at integration point level, which is assumed to be a function of the strain energy:

li ¼ kl /ðEÞ

ð48Þ

The optimal value of kl was found heuristically to be kl  1 for the choice of b1 and b2 of the Mehrotra update presented in this section. 7. Numerical results 7.1. Sheared membrane strip The first example involved a rectangular membrane strip under shear and tension loading, which is one of the few relatively

complex cases for which an analytical solution has been developed by Mansfield [7]. This problem was therefore used to verify the proposed method in this paper. A convergence analysis for the different mesh types and loading conditions was performed. The membrane strip, as shown in Fig. 1, was clamped at the bottom and loaded with applied displacements u and v in the x and y directions respectively. The straight lines in the strip represent the straight tension rays as obtained from the Mansfield’s solution. Three of the different mesh types used for convergence analysis are shown in Fig. 2. The mesh in Fig. 2a was designed to be wellconditioned using equilateral triangles. A mesh with a second order polynomial bias towards the corners is displayed in Fig. 2b. The bias was applied to improve the accuracy near the corners which are singular for pure shear loading as predicted by the Mansfield solution. The irregular mesh in Fig. 2c was used to test the performance of the proposed method in an ill-conditioned mesh. In this section several figures are presented in which the stress distribution in a membrane after deformation is shown. These fig^, which are related to S^ ures all contain true stresses, denoted by s as:

1 J

^ T s^ ¼ F SF

ð49Þ

in which J represents the volume ratio between the current and reference configurations. The contour plot of the first principal true stress s1 in the integration points for two different load cases is shown in Fig. 3. These are pure shear in Fig. 3a and combined shear–tension loading in Fig. 3b. The stress directions are depicted by lines drawn in the integration points, and they form straight lines as predicted by the analytical solution. The stresses in the membrane were positive everywhere, which means that the main constraint was satisfied. In the regions where s1 approaches zero the membrane is close to slacking. As shown in Fig. 3 the stress directions in these regions were not computed as accurately since they should form straight tension rays but appeared to be more randomly oriented. One explanation is that these regions are close to the boundary of the feasible domain of the optimization problem, S^ > 0, where the barrier function, see (17), has a significant effect on the solution. The second contribution to this randomness in stress directions is the stress–strain relation which is close to singular in the slacked regions due to the low effective membrane stiffness. The finite numerical precision caused round-off errors in this case. For the second case in Fig. 3 including tensile displacement, the stress directions were modeled more accurately near the slacking corners. This could be explained by the fact that the tensile

36

R. de Rooij, M.M. Abdalla / Computers and Structures 151 (2015) 30–41

(a) Tria elements

(b) Quad elements, biased mesh

(c) Quad elements, irregular mesh

Fig. 2. Example meshes of rectangular membrane strip.

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Fig. 3. Contour plot of first principal true stresses

0

1

R X

ðDrÞT C 1 ðDrÞdX /MF

1

1.5

2

2.5

3

s1 for two different load cases with a ¼ 50 mm, L ¼ 100 mm.

deformation caused the stresses to increase and thus to move away from the boundary of the feasible region of the stresses, and the effective membrane stiffness increases. A convergence study was performed to assess the performance of the numerical model more thoroughly using the analytical solution presented in Mansfield [7]. In the Mansfield solution the linear elastic Saint–Venant material model is assumed and the Poisson’s ratio m does not appear because infinitesimal displacements are assumed. In the numerical model the same material law is implemented with a Poisson’s ratio of m ¼ 0:5. This choice was made to resemble the neo-Hookean material law, commonly used to model thin membranes, at small deformations. The convergence analysis was based on the error in total strain energy in the membrane computed as:

e¼2

0.5

ð50Þ

in which /MF is the total strain energy in the analytical solution and Dr ¼ r  rMF . Fig. 4 contains four mesh convergence plots for the rectangular membrane, each of which checks different element properties. It shows that convergence was obtained in all cases and thus the proposed element was verified. Each convergence plot contains a triangle indicating first order convergence with respect to the element size h. Fig. 4 shows OðhÞ convergence in all cases, which corresponds to the convergence rate of regular finite elements. The errors shown in Fig. 4 leveled off at a value of approximately 102 , whereas the tolerance on the penalty parameter l was 106 . This discrepancy was partly explained by the finite displacements in the numerical model and the assumption of infinitesimal displacements in the Mansfield solution. Other explanations were found in the limitations of integrating the stress singularities at the corners numerically and numerical difficulties in modeling the slacked regions accurately. Despite the relatively high level off error, the results in Fig. 4 show a positive verification for convergence.

From Fig. 4a it was observed that for pure shear deformation the proposed method contained a larger error than for the second load case including tension. This was expected due to the numerical difficulties in the regions close to slacking. The quadrilateral elements also yielded a smaller error than the triangular elements and a smoother convergence curve. Note that the convergence curves for triangular elements contained several peaks, which was explained by the symmetry properties of the problem. Both the structure and the loading were point-symmetric with respect to the center of the membrane. The triangular mesh, see Fig. 2a, contains the same symmetry properties when the number of rows of elements is even. The mesh does not have these same symmetry properties when the number of element rows is odd. In those cases the mesh was less accurate to model the symmetry properties, resulting in a higher error. Applying a mesh bias towards the corners of the membranes reduced the error as can be seen from Fig. 4b. The mesh bias was based on the analytical solution, in which stress singularities occur at two corners of the membrane. These stress singularities are of pffiffiffi order 0.5, meaning they are of the form 1= r , with r being the distance from the corner. The applied mesh bias, see Fig. 2b, was generated using the same order of 0.5, denoted by ‘bias = 0.5’ in Fig. 4b. This led to a better integration of the stress singularities and a lower error. The convergence behavior for square (AR = 1) and rectangular (AR = 2) quad elements is shown in Fig. 4c. The observed differences are small but as expected, since the better-conditioned elements performed better. Finally the convergence for different shapes of the quad elements is presented in Fig. 4d. In the skewed mesh all the elements were skewed uniformly and formed parallelograms, the irregular mesh is depicted in Fig. 2c. As expected the square mesh performed best whereas the irregular mesh ceased to converge at a higher error level. The latter was caused by highly distorted elements in the irregular mesh, which are elements with very sharp or very blunt corners. Numerical integration using four integration points is less accurate and contains more numerical errors in these elements.

37

R. de Rooij, M.M. Abdalla / Computers and Structures 151 (2015) 30–41 Different load cases

Different mesh biases

−0.4

−0.6 T, u=0.3 v=0 Q, u=0.3 v=0 T, u=0.3 v=0.1 Q, u=0.3 v=0.1

−0.6

−1

−0.8

−1.2

log10 (ε)

−1

log10 (ε)

T, bias=0 Q, bias=0 T, bias=0.5 Q, bias=0.5

−0.8

−1.2 −1.4 −1.6

−1.4 −1.6 −1.8 −2

−1.8

−2.2

−2

−2.4

−2.2

−2.6 0

0.5

1

1.5

2

2.5

3

3.5

0

0.5

1

1.5

2

log10 (ne)

log10 (ne)

(a) Different load cases.

(b) Different mesh biases.

Different mesh aspect ratios

3

Different element shapes

−0.6

−0.6 Q, square mesh Q, skewed mesh Q, irregular mesh

Q, AR=1 Q, AR=2

−0.8

−0.8

−1

−1

−1.2

−1.2

log10 (ε)

log10 (ε)

2.5

−1.4

−1.4

−1.6

−1.6

−1.8

−1.8

−2

−2

−2.2

−2.2 0

0.5

1

1.5

2

2.5

3

0

0.5

1

1.5

2

log10 (ne)

log10 (ne)

(c) Different quad mesh aspect ratios.

(d) Different element shapes.

2.5

3

Fig. 4. Mesh convergence for rectangular membrane strip with a ¼ 50 mm, L ¼ 100 mm.

7.2. Annular square The performance of the proposed element was further assessed using two examples considered by Haseganu et al. [12] to test their numerical method for membrane wrinkling based on a dynamic relaxation method. The first of these examples consists of an annular square membrane clamped at the outer edge and loaded at the inner edge as depicted in Fig. 5. The meshes in this and following examples were optimized to be well-conditioned.

The applied displacement to the inner square consisted of two parts. First an in-plane twist of 90° was applied after which the square was pulled out normal to the reference configuration by W=2 ¼ 25 mm. The neo-Hookean material law was used for this example, as well as for the examples in the next sections. Fig. 6 shows the deformed configuration of the membrane including the membrane states based on those given in Table 2. Note that the characteristics in Table 2 are based on the second Piola– ^ which is qualitatively equivalent to the true Kirchhoff stress S,

Fig. 5. Input mesh of annular square membrane with W out ¼ 50 mm and W in ¼ 17:5 mm.

Fig. 6. Membrane states in annular square membrane under axisymmetric loading.

38

R. de Rooij, M.M. Abdalla / Computers and Structures 151 (2015) 30–41

Table 2 Characteristics of membrane states based on principal second Piola–Kirchhoff stresses r1 and r2 . State

Exact

Approximated

Taut Wrinkled Slack

r1 > 0 & r2 > 0 r1 > 0 & r2 ¼ 0 r1 ¼ 0 & r2 ¼ 0

r1 > 0:01r1;max & r2 = maxf1; r1 g > 0:01 r1 > 0:01r1;max & r2 = maxf1; r1 g < 0:01 r1 < 0:01r1;max & r2 = maxf1; r1 g < 0:01

^. The exact characteristics of the wrinkled and slack states stress s never occurred in the numerical results because the barrier function in (17) led to strictly positive stresses. To account for this observation, the approximated characteristics were introduced as a post-processing tool to visualize the membrane states. Note that in the approximated characteristics, r1 was normalized with respect to maxf1; r1 g to prevent normalization with respect to r1  0 in the slack regions. From Fig. 6 it was observed that the complete annular square membrane under the given loading was either taut or wrinkled and the pattern of the different states was relatively complex. 7.3. Sheared membrane with hole The second example from Haseganu et al. [12] to be considered involved severe slacking of the membrane. This example consisted of a square membrane with a circular hole loaded in shear. The initial configuration with a triangular mesh is presented in Fig. 7. The bottom edge of the membrane was clamped and a shear deformation waspffiffiffi applied to the other outer edges with a shear angle of c ¼ 1= 3  65°. The states of the membrane in the deformed configuration are shown in Fig. 8.

Fig. 7. Input mesh of square membrane with a hole with W ¼ 100 mm and Rin ¼ 10 mm.

Fig. 9. First principal true stresses

Stress concentrations around the hole caused taut and slack regions at those locations expected from regular plates with holes. The detailed behavior of the membrane in the slack regions was obtained by refining the mesh in these regions, which yielded the first principal true stresses, shown in Fig. 9. A severe mesh distortion is displayed in Fig. 9 in the slack region due to the large shear angle that has been applied. This distortion was explained by the severity of slack and the low effective membrane stiffness in these slack regions. The former led to a large inplane contraction of the membrane due to wrinkles and the low effective membrane stiffness caused the stress–strain relation to become close to singular. Note that the results presented in this section were obtained using a two dimensional analysis, because the current example is numerically unstable in three dimensions. This instability is due to the very low out of plane stiffness of the membrane in slack regions, which is explained by the expression for the element stiffness matrix in (37). The first term within the brackets in (37) represents the geometric stiffness term, which is the sole contributor to the out of plane stiffness of the membrane. This contribution disappears in slack regions where S^  0. As a consequence, the slack regions of the membrane are numerically unstable in the out of plane direction. Note that recovery from the numerical instability is theoretically guaranteed by the convexity properties of the method, however, it is numerically expensive or unfeasible due to numerical errors. 7.4. Pressurized beam under bending In this example a cylindrical, pressurized beam with a lateral tip displacement is considered. Pressurization of the beam led to a fully taut membrane whereas the bending load compressed one side of the beam which caused wrinkling when the load exceeded a critical level. The wrinkling led to a softening of the structural bending stiffness of the beam, which was modeled and compared with experimental data. This example was included to serve as a validation of the element and it resembles practical applications such as a pressurized air balloon. To model the volume dependent pressure loading, several contributions needed to be added to the stiffness matrix and residual force vector. These are described in Rumpel [28] and are presented briefly here. The external virtual work done by pressure loading dPpress is given as:

dPpress ¼ Fig. 8. Membrane states in square membrane with a hole under shear loading.

s1 near the deformed hole.

Z

x1

Z

1

2

pn  dudx dx x2

ð51Þ

39

R. de Rooij, M.M. Abdalla / Computers and Structures 151 (2015) 30–41

in which p is the pressure, n ¼ a1  a2 according to the convention in Section 2, and x1 and x2 are the coordinates in the current configuration using this same convention. The governing equation including this pressure contribution then becomes:

dPstruct þ dPpress ¼ 0

ð52Þ

with dPstruct being the variation of potential energy of the membrane structure as in (20). The pressure p is determined by Poisson’s law as function of the enclosed volume V:

p0 V j0 ¼ pV j

ð53Þ

with p0 and V 0 being the pressure and enclosed volume in the reference configuration and j is the isentropy constant which takes the value j ¼ 1:4 in air. The volume V of any enclosed region is computed as:

V¼

1 3

Z

Z x1

1

2

r  n dx dx

ð54Þ

x2

The current volume of the pressurized cylinder is computed by evaluating (54) over all three bounding surfaces and adding their contributions. Linearization and discretization of the pressure contributions yields the set of governing equations:



 K struct  K press þ cd  d Dv ¼ Rstruct þ f press

ð55Þ

in which K press is the pressure load stiffness matrix, c ¼ j Vp ; d  d represents the volume change and f press is the pressure load vector. The cd  d-term in (55) ensures that the Poisson’s law in (53) is satisfied. For more detailed expressions, consult Rumpel [28].

The cylindrical beam considered in this example was clamped at the root and its circular end cap was modeled as a rigid part. The results were validated using experiments performed by Veldman [29], in which the applied material was PolyPhenyleneSulphide (PPS) with a stiffness modulus E ¼ 2850 MPa and a measured Poisson’s ratio between m ¼ 0:347 and m ¼ 0:351. The neo-Hookean material model was applied because it is numerically slightly more stable than the linear elastic model and yielded the same results for validation purposes. The cylindrical beam had a length of L ¼ 1000 mm, a radius of 50 mm and the membrane thickness is t ¼ 0:06 mm. The deformed configuration of a beam loaded in bending displacement, which was applied using a vertical deflection in the center of the rigid end cap, is shown in Fig. 10. The first principal stress was tensile everywhere due to the pressure loading, which means that slack did not occur. From Fig. 10b it can be seen that wrinkling occured in the compressed side of the beam where s2  0. The load–displacement curves of the pressurized beam for three different internal pressures are presented in Fig. 11. Initially the curves were linear, but after wrinkling initiated softening of the structural bending stiffness occurred. The load at which wrinkling occured was proportional to the internal pressure which caused the tensile stresses. The initial slope of the load–displacement curve was independent of the internal pressure in the beam. The numerical results in Fig. 11 were validated using experimental and analytical data obtained from Veldman [29]. In the linear regime the numerical results matched these data well, whereas the model underestimated the required force in the softening regime. This could be attributed to the assumption of zero bending

10

10

9

9

8

8

7

7

6

6

Force [N]

Force [N]

Fig. 10. First and second principal true stresses of cylindrical beam with p ¼ 0:02 MPa and a tip deflection of 100 mm.

5 4

5 4 3

3

2

2 Experiment Membrane Theory Proposed Method

1

Experiment Membrane Theory Proposed Method

1 0

0 0

20

40

60

80

Deflection [mm]

100

120

0

20

40

60

80

100

120

Deflection [mm]

Fig. 11. Load deflection curves for cylindrical pressurized membrane beam under bending used for validation of numerical results with experimental data. The beam has L ¼ 1000 mm, r ¼ 50 mm and t ¼ 0:06 mm.

40

R. de Rooij, M.M. Abdalla / Computers and Structures 151 (2015) 30–41

stiffness made in tension field theory, which is an underestimation of the actual, small bending stiffness. The slope of the numerical results in a softening regime did match the experimental data, which means that the ‘post-buckling’ stiffness was modeled accurately.

for 0 6 a 6 1. The last term in this argument is positive and, given that /ðEÞ is non-decreasing and convex, this expression can be written as:

wðð1  aÞF 1 þ aF 2 Þ 6 /ðð1  aÞE1 þ aE2 Þ 6 ð1  aÞ/ðE1 Þ þ a/ðE2 Þ ¼ ð1  aÞwðF 1 Þ þ awðF 2 Þ

8. Conclusion In this paper, an implementation method of the relaxed strain energy function for membranes in the presence of wrinkles was presented and applied based on the interior-point method. The relaxed strain energy function was reformulated as a constrained, convex optimization problem, suitable for the interior-point method implementation. This implementation was performed and the variational principle was applied to obtain the set of governing equations for membranes. By definition of the convexity properties, these equations uniquely determined the deformation of the membrane and the stresses and wrinkling strains in membranes, within a finite element context. The different membrane states, taut, wrinkled and slack, are captured implicitly in the governing equations and therefore no if-else statements are required to distinguish between them. An efficient predictor–corrector method was applied to obtain the solution to the governing equations. Verification of the solution was performed by modeling a rectangular membrane strip under shear loading and comparing the results with the analytical solution. An OðhÞ convergence rate of the energy norm was shown. Validation of the proposed method was achieved by modeling a pressurized cylindrical beam under bending displacement. The load–deflection curve was compared with experimental data and showed to agree in the linear regime of this curve. An underestimation of the structural bending stiffness in the presence of wrinkles was observed due to the zero-bending stiffness of the membrane assumed in the tension field theory method. It was shown that the proposed method efficiently, robustly and accurately describes the stress distribution in membranes in the presence of wrinkles. Appendix A. Convexity conditions The relaxed energy function, denoted /] ðEÞ ¼ w] ðFÞ, is the quasi-convexification of wðFÞ if it is the largest quasi-convex function that does not exceed /ðEÞ; here we present the conditions on /] ðEÞ to be satisfied for this quasi-convexification. The Pipkin’s [10] approach is followed. Shorter proofs are provided, including a more convenient expression of the relaxed energy functional. For the following, define a function f ðXÞ to be non-decreasing iff:

f ðX þ ZÞ P f ðXÞ 8Z P 0 Proposition 1. If /ðEÞ is convex and non-decreasing, then wðFÞ is convex.

which proves convexity of wðFÞ. In this expression the first inequality holds as /ðEÞ is non-decreasing and the second inequality exploits the convexity of /ðEÞ. h Proposition 2. If non-decreasing.

wðFÞ

is

rank-one

convex,

then

/ðEÞ

is

Proof. Rank-one convexity of wðFÞ means that for F 1 ¼ F and F 2 ¼ F þ ueT convexity conditions are satisfied, with u and e being any arbitrary vectors:

    w F þ aueT 6 ð1  aÞwðFÞ þ aw F þ ueT

ðA:4Þ

with 0 6 a 6 1. Now choose u such that F T u ¼ 0 which is possible since F forms a 3  2 matrix, as shown in Section 2. Choosing without loss of generality uT u ¼ 2 it can be shown similarly to (A.2) that:

    w F þ aueT ¼ / E þ a2 eeT Substitution into (A.4) then gives:

a

  / E þ a2 eeT  /ðEÞ

a2

  6 / E þ eeT  /ðEÞ

ðA:5Þ

Taking the limit a ! 0þ on both sides causes the left hand side of this expression to approach zero. This follows as:

lim

a!0þ

  / E þ a2 eeT  /ðEÞ

a2

represents the directional derivative of /ðEÞ in the direction e, which is bounded provided /ðEÞ is Lipschitz continuous. In the limit (A.5) then becomes:

  / E þ eeT P /ðEÞ Any semi-positive definite (spd) tensor DE P 0 can be written T as DE ¼ aaT þ bb , from which:

  T /ðE þ DEÞ ¼ / E þ aaT þ bb P /ðEÞ 

Proposition 3. If /ðEÞ is convex, then its lower, non-decreasing envelope /] ðEÞ is the quasi-convexification of wðFÞ. Proof. Propositions 1 and 2 together with convexity properties [24] lead to the following implications:

convex þ non-decreasing /ðEÞ

Proof. The following equality is applied:

) convex wðFÞ

ðð1  aÞF 1 þ aF 2 ÞT ðð1  aÞF 1 þ aF 2 Þ ¼ ð1  aÞF T1 F 1 þ aF T2 F 2  að1  aÞDF T DF

ðA:3Þ

ðA:1Þ

with DF ¼ F 1  F 2 . From the definition of the Green–Lagrange strain in (4), and (A.1) one gets:

h i

1 wðð1  aÞF 1 þ aF 2 Þ ¼ / ðð1  aÞF 1 þ aF 2 ÞT ðð1  aÞF 1 þ aF 2 Þ  I 2

1 ¼ / ð1  aÞE1 þ aE2  að1  aÞDF T DF 2 ðA:2Þ

) quasi-convex wðFÞ ) rank-one convex wðFÞ ) non-decreasing /ðEÞ The quasi-convexity of wðFÞ guarantees the existence and uniqueness of a displacement field carrying the minimum strain energy. The sufficient condition for this quasi-convexity is a convex and non-decreasing /ðEÞ and the necessary condition that /ðEÞ is non-decreasing. For convex /ðEÞ the lower non-decreasing envelope is both sufficient and necessary for quasi-convexity and thereby is exactly be the quasi-convexification of wðFÞ. h

R. de Rooij, M.M. Abdalla / Computers and Structures 151 (2015) 30–41

41

Appendix B. Relaxation formulation

References

In Proposition 3 it is shown that the quasi-convexification, or relaxation, of the strain energy is given by the lower, convex and non-decreasing envelope of /ðEÞ. In this section we will show how to obtain this relaxation for any arbitrary function f ðXÞ. The notation X : Y ¼ trX T Y is applied.

[1] Wagner H. Ebene blechwandträger mit dünnem stegblech. Z Flugtech Motorluftschif 1929;20:200–314. [2] Puntel E, Deseri L, Fried E. Wrinkling of a stretched thin sheet. J Elast 2011;105(1–2):137–70. [3] Vella D, Ajdari A, Vaziri A, Boudaoud A. Wrinkling of pressurized elastic shells. Phys Rev Lett 2011;107(17). [4] Flores FG, Oñate E. Wrinkling and folding analysis of elastic membranes using an enhanced rotation-free thin shell triangular element. Finite Elem Anal Des 2011;47(9). [5] Weinberg K, Neff P. A geometrically exact thin membrane model – investigation of large deformations and wrinkling. Int J Numer Methods Eng 2008;74(6):871–93. [6] Pagitz M, Abdalla M. Simulation of tension fields with in-plane rotational degrees of freedom. Comput Mech 2010;46(5):747–57. [7] Mansfield EH. Load transfer via a wrinkled membrane. Proc Roy Soc London. Ser A, Math Phys Sci 1970;316(1525):269–89. [8] Mansfield EH. Analysis of wrinkled membranes with anisotropic and nonlinear elastic properties. Proc Roy Soc London Ser A 1977;353(1975):475–98. [9] Pipkin AC. Relaxed energy densities for small deformations of membranes. IMA J Appl Mat (Inst Math Appl) 1993;50(3):225–37. [10] Pipkin AC. Relaxed energy densities for large deformations of membranes. IMA J Appl Math (Inst Math Appl) 1994;52(3):297–308. [11] Steigmann DJ. Tension-field theory. Proc Roy Soc London. Ser A, Math Phys Sci 1990;429(1876):141–73. [12] Haseganu E, Steigmann D. Analysis of partly wrinkled membranes by the method of dynamic relaxation. Comput Mech 1994;14(6):596–614. [13] Akita T, Nakashino K, Natori MC, Park KC. A simple computer implementation of membrane wrinkle behaviour via a projection technique. Int J Numer Methods Eng 2007;71(10):1231–59. [14] Jarasjarungkiat A, Wüchner R, Bletzinger K. Efficient sub-grid scale modeling of membrane wrinkling by a projection method. Comput Methods Appl Mech Eng 2009;198(9-12):1097–116. [15] Roddeman D, Drukker J, Oomens C, Janssen J. The wrinkling of thin membranes: Part i – theory. J Appl Mech 1987;54(4):884–7. [16] Schoop H, Taenzer L, Hornig J. Wrinkling of nonlinear membranes. Comput Mech 2002;29(1):68–74. [17] Raible T, Tegeler K, Löhnert S, Wriggers P. Development of a wrinkling algorithm for orthotropic membrane materials. Comput Methods Appl Mech Eng 2005;194(21-24 SPEC. ISS.):2550–68. [18] Epstein M. On the wrinkling of anisotropic elastic membranes. J Elast 1999;55(2):99–109. [19] Kang S, Im S. Finite element analysis of dynamic response of wrinkling membranes. Comput Methods Appl Mech Eng 1999;173(12):227–40. [20] Miyazaki Y. Wrinkle/slack model and finite element dynamics of membrane. Int J Numer Methods Eng 2006;66(7):1179–209. http://dx.doi.org/10.1002/ nme.1588. [21] Mosler J, Cirak F. A variational formulation for finite deformation wrinkling analysis of inelastic membranes. Comput Methods Appl Mech Eng 2009;198(27–29):2087–98. [22] Ball JM. Strict convexity, strong ellipticity, and regularity in the calculus of variations. Math Proc Camb Philos Soc 1980;87(03):501–13. http://dx.doi.org/ 10.1017/s0305004100056930. [23] Boyd S, Vandenberghe L. Convex optimization. New York, NY, USA: Cambridge University Press; 2004. [24] Schröder J, Neff P. Poly-, quasi- and rank-one convexity in applied mechanics. Springer; 2010. [25] Vandenberghe L, Boyd S. Semidefinite programming. SIAM Rev 1994;38:49–95. [26] Yamashita M, Fujisawa K, Kojima M. Implementation and evaluation of SDPA 6.0 (semidefinite programming algorithm 6.0). Optim Methods Softw 2003;18(4 II):491–505. [27] Mehrotra S. On the implementation of a primal–dual interior point method. SIAM J Optim 1992;2(4):575–601. [28] Rumpel T, Schweizerhof K. Volume-dependent pressure loading and its influence on the stability of structures. Int J Numer Methods Eng 2003;56(2):211–38. [29] Veldman S. Design and Analysis Methodologies for Inflated Beams. DUP Science; 2005. [30] Rockafellar R. Convex analysis. Princeton mathematical series. PRINCETON University Press; 1997.

Proposition 4. Let f ðXÞ be a scalar-valued function of a symmetric ] tensor X and let f ðXÞ be the largest, increasing and convex function ] with f ðXÞ 6 f ðXÞ 8X. Then: 

]

f ðXÞ ¼ sup½X : Y  f ðYÞ

ðB:1Þ

YP0 

with f ðYÞ being the complementary function of f ðXÞ defined by the Legendre–Fenchel transform: 

f ðYÞ ¼ sup½X : Y  f ðXÞ X ]

Proof. The function f ðXÞ is referred to as the relaxation of f ðXÞ. ] Convexity of f ðXÞ is obtained as a property of the Legendre–Fen] chel transformation [30]. In order to prove f ðXÞ is non-decreasing, let Z P 0 and apply the proposition: 

]



f ðXÞ ¼ sup½Y : X  f ðYÞ ¼ sup½Y : ðX þ ZÞ  f ðYÞ  Y : Z  YP0

YP0 

]

6 sup½Y : ðX þ ZÞ  f ðYÞ þ sup½Y : Z  ¼ f ðX þ ZÞ YP0

ðB:2Þ

YP0

From Proposition 4 it is evident that: ]





f ðXÞ 6 f ðXÞ ¼ sup½Y : X  f ðYÞ Y

Applying the Legendre–Fenchel transform and using:

sup inf hðxÞ 6 inf sup hðxÞ 

for any scalar valued function, it can be shown that f ðXÞ 6 f ðXÞ. ] Consequently, the third property follows: f ðXÞ 6 f ðXÞ. ] Finally, it is left to prove that f ðXÞ is the largest function with has these three proven properties. For this, let gðXÞ be any convex, nondecreasing function with gðXÞ 6 f ðXÞ 8X, meaning gðXÞ satisfies all proven properties and can be shown to give g ] ðXÞ ¼ gðXÞ. Using the definition of the Legendre–Fenchel transformation one gets: 

f ðYÞ ¼ sup½X : Y  f ðXÞ X

6 sup½X : Y  gðXÞ þ sup½gðXÞ  f ðXÞ 6 g  ðYÞ X

ðB:3Þ

X ]

Similarly it can be shown that g ] ðXÞ 6 f ðXÞ, which means that ] g ðXÞ ¼ gðXÞ 6 f ðXÞ 6 f ðXÞ. h ]

Proposition 5. Let f ðXÞ be a convex, scalar-valued function of a sym] metric tensor X and let f ðXÞ be the largest, increasing and convex ] function with f ðXÞ 6 f ðXÞ 8X. Then: ]

f ðXÞ ¼ sup inf ½f ðZÞ  Y : ðZ  XÞ

ð4Þ

YP0 Z

Proof. Applying the definition of the relaxation in Eq. (B.1) yields for a convex f ðXÞ:

 ] f ðXÞ ¼ sup Y : X  sup½Y : Z  f ðZÞ YP0



Z

¼ sup Y : X þ inf ½f ðZÞ  Y : Z  YP0



Z

¼ sup inf ½f ðZÞ  Y : ðZ  XÞ  YP0 Z

ð5Þ