# Numerical wrinkling prediction of thin hyperelastic structures by direct energy minimization

Yann Lecieux, Rabah Bouzidi

► To cite this version:

Yann Lecieux, Rabah Bouzidi. Numerical wrinkling prediction of thin hyperelastic structures by direct energy minimization. Advances in Engineering Software, 2012, 50, pp.57 - 68. 10.1016/j.advengsoft.2012.02.010. hal-01007136

HAL Id: hal-01007136

https://hal.science/hal-01007136

Submitted on 3 Dec 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical wrinkling prediction of thin hyperelastic structures by direct energy minimization

Yann Lecieux, Rabah Bouzidi

L'UNAM Université, Université de Nantes-Ecole Centrale Nantes, GeM, Institut de Recherche en Génie Civil et Mécanique, UMR CNRS 6183, 2 rue de la Houssinière, BP 92208 44322, Nantes cedex 3, France

This paper is concerned with an efficient algorithm for the wrinkling with finite strains of very thin structures made of hyperelastic material. In this work, the problem of wrinkling is solved by directly minimizing the total potential energy of the structure. The numerical solution is carried out by the means of an iterative method such as the conjugate gradient algorithm. Although the proposed approach is theoretically equivalent to the traditional finite element method, it proves to be an attractive alternative which is particularly efficient for thin wrinkled structures.

**Keywords:** Energy minimization, Wrinkling, Rotation-free triangular shell element Membranes experiments, Buckling, Biaxial load

## 1. Introduction

Over the past years, new structural concepts for large spacecraft applications involving thin film surface have been designed. Compared to traditional spacecraft structures, these "Gossamer structures" could provide many advantages such as reduced mass and package volume. However, the materials used in gossamer structures (as very thin Kapton® films) cannot support compressive stress because of their small bending rigidity. The result of compressive stress is that buckling occurs and wrinkles are formed. The existence of wrinkled regions may affect the performance and reliability of the flexible gossamer structures (as in an antenna or a reflector). Thus, the prediction of wrinkle patterns in a membrane surface is one of the many current technological interests in aerospace industry.

This paper deals with the use of the first order direct energy minimization approach to solve the wrinkling problem in thin structures. Although the proposed approach is theoretically equivalent to the traditional finite element method, it is an attractive alternative which is particularly efficient for thin wrinkled structures. In order to demonstrate it, a comparison is drawn with data resulting from two different studies [1,2], where the wrinkling problem has been studied experimentally as well as numerically using a bifurcation and a post-buckling analysis.

The main advantage of the first order minimization technique compared to those based on bifurcation analysis is that no buckling analysis is needed to achieve the wrinkling simulation. In order to place the structure out of the equilibrium path, a small random out

of plane displacement is imposed to the membrane at the beginning of the simulation. This is necessary to avoid staying on the fundamental path following a bifurcation point. The wrinkled shape is thus automatically obtained by minimizing the total potential energy of the structure.

The paper is organized as follows: Section 2 presents a review of previous numerical studies on the topic of thin elastic membrane wrinkling. In Section 3, the strain energy of thin shell is depicted. Then in Section 4, the prediction of wrinkling with a first order algorithm is explained. The finite elements used here are presented in Section 5. The performance of the proposed method to predict wrinkle are illustrated in Section 6. Then Section 7 concludes the paper.

## 2. Review of numerical procedures of wrinkling simulation

The analysis of the wrinkling behavior of membrane structures started with the works of Wagner [3], who initiated the method called Iterative Material Procedures (IMPs). This method only takes into account the in-plane stiffness. The bending stiffness is not evaluated. These wrinkling models are mainly based on the deformation gradient or the constitutive equation modification in order to avoid compressive stresses [4,5]. As a result, this kind of approach does not make it possible to represent wrinkle patterns and does not provide any information about the size of wrinkles. This aspect is a major disadvantage of the IMP procedure. However, the stress field and wrinkling zone are correctly determined. Mansfield in [6] and Pipkin in [7] reformulated the theory by using a suitable relaxed energy density as defined in the Tension Field Theory. The relaxed energy density represents the average energy per initial area unit over a region containing many wrinkles.
To enable the computational simulation to predict the wrinkle patterns, two methods are classically used. The first one involves the treatment of bifurcation points with asymptotic expansion or a similar process [8]. This method involves the detection of critical points based on the singularity of the tangent stiffness matrix, then the switching on a bifurcated branch. Its detailed developments can be found in [9–11].

The second one is the post buckling analysis [1]. To enable the computational method to predict the wrinkle patterns, the use of an extensively refined mesh with thin-shell elements possessing membrane and bending stiffness is necessary. This kind of elements enables to precisely analyze the buckling and post buckling response of the stretched membrane. This specific treatment of the buckling in the classical finite element method usually leads to heavily computational times due to the dense mesh.

A wrinkling analysis, using the buckling of shell element, is usually performed in three stages: the first consists in obtaining a stable initial state in the case of very thin shells. This is often achieved by applying slight initial pre-stress that increased the low bending stiffness of thin shells by the geometrical bending stiffness. The second step is an eigenvalue buckling analysis, which gives the membrane mode shapes, introduced as geometrical imperfections in the third step of the post buckling analysis.

Several recent computational studies have employed the post buckling analysis with geometrically nonlinear shell finite element models, see [1,12,2]. Wong and Pellegrino [1] performed a numerical analysis of membrane wrinkling in a rectangular membrane under shearing load and the wrinkling prediction in a square membrane under corner loads where the membrane was modeled as a thin shell with a near-zero stiffness.

The main difficulty in the post buckling analysis is the choice of a possible wrinkling mode. If the final result is approximately known, the convergence can be achieved by finding the buckling mode that matches most to the final pattern. This is the solution chosen by Wong and Pellegrino in [1]. Nevertheless, if the objective of the study is to find a possible wrinkle shape without knowing the final pattern, choosing to introduce a given buckling mode is difficult. Using the first wrinkling mode corresponding to the lowest level of energy may be considered as a solution to this issue. However the presence of slack zones can disturb the bifurcation modes and consequently the use of the post buckling method needs an excellent understanding of the physical problem to achieve some good simulation results.

In the present study, an alternative approach to those mentioned above is presented. It is based on the minimization of the energy by descent method. The main interest of this approach is that no bifurcation analysis is needed before performing the wrinkling computation. The method is sufficiently robust to find a possible wrinkle pattern after applying a small random displacement at each node of the membrane surface.

### 3. The approach of total potential energy minimization

This section is concerned with the presentation of the strain energy of a thin film submitted to external loads. Then, the fundamentals of the minimization of the total potential energy are explained.

#### 3.1. The internal potential energy

Engineering applications of thin film involves small strain and large deformations. In this study, only the elastic response of material is considered. Indeed, in wrinkling problems, the effects of large deformations are primarily due to large rotations. Thus, we model the elastic comportment of material by a Saint Venant–Kirchhoff material. Using this model and considering the kinematic of thin shell, the expression of the strain energy of a thin film is obtained.

#### 3.1.1. Kinematics of deformation

The geometry of the structure is depicted using the natural coordinate system {$\xi_1, \xi_2, \xi_3$}. We limit our attention to the theory of thin shells which is based on the approximation that the displacement is assumed to vary linearly through the thickness. Then, we specialize this classical shell theory to Kirchhoff–Love theory by explicitly containing the shell director to remain normal to the deformed middle surface of the shell. Here, only the main results are presented. Furthers details and discussions may be found in [13,14].

We begin by considering a shell whose initial geometry (at time $t_0$) is characterized by a middle surface of domain ($\Omega(t_0)$). The deformed or actual configuration of the shell (at time $t$) is characterized by a middle surface of domain $\Omega(t)$ with ($\Omega(t_0)$, $\Omega(t) \subseteq \mathbb{R}^3$). $\vec{X}$ and $\vec{x}$ are the position vectors of the same material particle in the reference and deformed configuration. Curvilinear coordinates $\xi^1, \xi^2, \xi^3$ are used to describe the position of a thin shell characterized by its middle surface $\xi^1, \xi^2, 0$ and its thickness $\xi^3$, $(-\frac{h}{2} \le \xi^3 \le \frac{h}{2})$ assumed uniform as: (see Fig. 1)

$$ \vec{X}(\xi^1, \xi^2, \xi^3) = \vec{X}(\xi^1, \xi^2, 0) + \xi^3 \vec{A}_3(\xi^1, \xi^2) \quad (1) $$

$$ \vec{x}(\xi^1, \xi^2, \xi^3) = \vec{x}(\xi^1, \xi^2, 0) + \xi^3 \vec{a}_3(\xi^1, \xi^2) \quad (2) $$

$\vec{A}_3$ is the unit normal of the shell reference surface while the corresponding vector on the actual surface $\vec{a}_3$ is not a unit vector nor normal. This description undergoes the thin shell hypothesis for which the kinematics is linear over the thickness.

We note the curvilinear base $\vec{G}_i = \vec{X}_i$ in the initial configuration, and $\vec{g}_i = \vec{x}_i$ in the actual configuration. They can be written as:

$$ \begin{aligned} \vec{G}_{\alpha} &= \vec{A}_{\alpha} + \xi^3 \frac{\partial \vec{A}_3}{\partial \xi^{\alpha}} = \vec{X}(\xi_1, \xi_2, 0)_{,\alpha} + \xi^3 \vec{A}_{3,\alpha} \quad \text{and} \quad \vec{G}_3 = \vec{A}_3 \\ \vec{g}_{\alpha} &= \vec{a}_{\alpha} + \xi^3 \frac{\partial \vec{a}_3}{\partial \xi^{\alpha}} = \vec{x}(\xi_1, \xi_2, 0)_{,\alpha} + \xi^3 \vec{a}_{3,\alpha} \quad \text{and} \quad \vec{g}_3 = \vec{a}_3 \end{aligned} \quad (3) $$

where the Greek indices have the range {1, 2} while the Latin indices have the range {1, 2, 3}. $\vec{A}_i$ and $\vec{a}_i$ are the material basis vectors on the middle surface i.e.

$$ \begin{aligned} \vec{A}_i(\xi^1, \xi^2) &= \vec{G}_i(\xi^1, \xi^2, 0) \quad \text{and} \\ \vec{a}_i(\xi^1, \xi^2) &= \vec{g}_i(\xi^1, \xi^2, 0) \end{aligned} \quad (4) $$

The contravariant basis in initial configuration $\vec{G}^1, \vec{G}^2, \vec{G}^3$ and actual configuration $\vec{g}^1, \vec{g}^2, \vec{g}^3$ are defined as:

$$ \begin{aligned} \vec{G}^i \cdot \vec{G}_j &= \delta_j^i \quad \text{and} \quad \vec{g}^i \cdot \vec{g}_j = \delta_j^i \\ \vec{A}^i \cdot \vec{A}_j &= \delta_j^i \quad \text{and} \quad \vec{a}^i \cdot \vec{a}_j = \delta_j^i \end{aligned} \quad (5) $$

The metric tensors components are defined in the initial and actual configurations respectively by:

$$ \begin{aligned} G_{ij} &= \vec{G}_i \cdot \vec{G}_j \quad \text{and} \quad g_{ij} = \vec{g}_i \cdot \vec{g}_j \\ G^{ij} &= \vec{G}_i \cdot \vec{G}^j \quad \text{and} \quad g^{ij} = \vec{g}^i \cdot \vec{g}^j \end{aligned} \quad (6) $$

The components of the Green–Lagrange strain tensor are:

$$ E_{ij} = \frac{1}{2} (g_{ij} - G_{ij}) \quad (7) $$

The mixed components of this tensor are then:

$$ \begin{aligned} E_i^j &= E_{ik} G^{kj} = \frac{1}{2} (g_{ik} G^{kj} - G_{ik} G^{kj}) = \frac{1}{2} (g_i^j - \delta_i^j) \\ E_i^j &= \frac{1}{2} (C_i^j - \delta_i^j) \end{aligned} \quad (8) $$

where $C_i^j$ are the mixed components of the right stretch tensor.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Base vectors of the deformed system.</figcaption></figure>

### 3.1.2. Thermodynamic potential

The material is assumed elastic and governed by a quadratic potential of strain tensor $\bar{E}$ defined per unit volume as:

$$ \psi = \frac{1}{2} H^{ijkl} E_{ij} E_{kl} \ge 0 \quad (9) $$

where $\overline{\overline{H}}$ is the fourth-order tensor of elastic moduli which are constants.

The second Piola-Kirchhoff stress tensor $\overline{\overline{\Sigma}}$ is given by:

$$ \overline{\overline{\Sigma}} = \frac{\partial \psi(\overline{\overline{E}})}{\partial \overline{\overline{E}}} = \overline{\overline{H}} : \overline{\overline{E}} \quad (10) $$

In this study, materials are modeled as isotropic. Consequently, the fourth-order tensor of elastic moduli $\overline{\overline{H}}$ reduces to:

$$ H^{ijkl} = \lambda G^{ij} G^{kl} + \mu (G^{ik} G^{jl} + G^{il} G^{jk}) \quad (11) $$

where $\lambda$ and $\mu$ are the Lamé constants.

It follows from this relation, that the expression of the strain energy may be written as:

$$ \psi = \frac{\lambda}{2} (\operatorname{tr} \overline{\overline{E}})^2 + \mu \overline{\overline{E}} : \overline{\overline{E}} \quad (12) $$

which is the classical expression of the Saint-Venant Kirchhoff model.

### 3.1.3. Strain energy of a thin shell

In the case of the Kirchhoff Love theory, the deformed vector $\overrightarrow{a_3}$ in the actual configuration coincides with the unit normal to the deformed surface of the shell. This hypothesis and the kinematic linearization over the thickness lead to the following expression of the Green Lagrange strain tensor:

$$ E_{\alpha\beta} = \chi_{\alpha\beta} + \xi^3 \rho_{\alpha\beta} \quad (13) $$

where

$$ \chi_{\alpha\beta} = \frac{1}{2} (\overrightarrow{a}_{\alpha} \cdot \overrightarrow{a}_{\beta} - \overrightarrow{A}_{\alpha} \cdot \overrightarrow{A}_{\beta}) \quad (14) $$

$$ \rho_{\alpha\beta} = \overrightarrow{A}_3 \cdot \overrightarrow{A}_{\alpha,\beta} - \overrightarrow{a}_3 \cdot \overrightarrow{a}_{\alpha,\beta} \quad (15) $$

The components $\chi_{\alpha\beta}$ measure the plane strain of the middle surface and $\rho_{\alpha\beta}$ the curvature change. Here, the transversal shear and the normal stretching have been neglected, consequently, $G^{\alpha3} = G^{3\alpha} = g^{\alpha3} = g^{3\alpha} = 0$, $G^{33} = g^{33} = 1$.

The shell strain energy density per unit area is obtained by summing the volumetric density over the thickness:

$$ \Psi = \frac{1}{2} \int_{-\frac{h}{2}}^{\frac{h}{2}} H^{ijkl} E_{kl} E_{ij} d\xi^3 \quad (16) $$

which leads to the following expression:

$$ \Psi = H^{\alpha\beta\lambda\mu} \left[ \frac{h}{2} \alpha_{\alpha\beta} \alpha_{\lambda\mu} + \frac{h^3}{24} \rho_{\alpha\beta} \rho_{\lambda\mu} \right] \quad (17) $$

where the material constants at the mid-surface are defined by:

$$ \overline{\overline{H}}^{\alpha\beta\lambda\mu} = \lambda A^{\alpha\beta} A^{\lambda\mu} + \mu (A^{\alpha\lambda} A^{\beta\mu} + A^{\alpha\mu} A^{\beta\lambda}) \quad (18) $$

The kinematic hypothesis of thin shell enables to write the strain energy as the sum of membrane strain energy and a bending one. These two energies are thus decoupled even though a large deformation formulation is considered.

A finite strain formulation is considered here, although the deformations stay infinitesimal for the wrinkling problems studied in Section 6. This framework formulation allows to consider rubber-like materials which can be submitted to finite strains.

## 3.2. Total potential energy theorem

The equilibrium of the structure is stated by the principle of minimum total potential energy. The structure shall deform to a position that minimizes the total potential energy. We consider conservative structures for which no dissipative loads are present. All the loads applied to the structure derive from a potential $\Pi^{\text{ext}}$ and the constitutive law is hyperelastic, so that the stresses derive from the internal potential energy $\Pi^{\text{int}}$.

For a general membranous structure with a volume V and mid-surface S, the total potential energy, $\Pi$, is the sum of the internal energy of the body $\Pi^{\text{int}}$:

$$ \Pi^{\text{int}} = \int_V \psi dV = \int_S \Psi dS \quad (19) $$

and the external potential energy $\Pi^{\text{ext}}$:

$$ \Pi^{\text{ext}} = - \int_S \overline{T} \cdot \overline{u} dS - \int_V \overline{f} \cdot \overline{u} dV \quad (20) $$

Here $\overline{u}$ is the displacement field of the structure, $\overline{f}$ represents prescribed body forces and $\overline{T}$ pressures acting at S. So that:
$$ \Pi = \Pi^{\text{int}} + \Pi^{\text{ext}} \qquad (21) $$

For equilibrium configuration $\bar{x}$, the energy change $\delta\Pi$ should be stationary irrespective of a small displacement perturbation $\delta\bar{u}$.

The discrete formulation of the internal energy $\Pi^{\text{int}}$ is given in Eq. (47) and discussed in Section 5 while the discrete form of $\Pi^{\text{ext}}$ is obtained by:

$$ \Pi^{\text{ext}}(\{u\}) = -\{u\}^T \{R\} \qquad (22) $$

where $\{u\}$ and $\{R\}$ are respectively the vector of nodal displacements and the vector of reaction forces applied at each node of the whole structure. Supposing that we have a one-parameter load conservative system, then $\{R\} = \lambda\{q\}$ where $\{q\}$ is the reference load vector and $\lambda$ a scalar load multiplier.

Our attention is limited to a single control parameter. If multiple control parameters $\{\lambda\}$ are needed, the load process is staged and then linearized at each stage. Between two states $\{\lambda\}_i$ and $\{\lambda\}_{i+1}$, the variations of all the components of $\{\lambda\}$ are assumed to stay proportional to the control parameter $\lambda$ such as:

$$ \{\lambda\} = (1 - \lambda)\{\lambda\}_i + \lambda\{\lambda\}_{i+1} \qquad (23) $$

where $\lambda$ is a one stage control parameter and varies from 0 to 1.

### 4. Use of a first order method of energy minimization to predict the phenomenon of wrinkling

The energy minimization can be achieved either by first order methods, like descent methods, or second order methods, like quasi-newtonian ones. To explain the fundamentals of these two kinds of methods, we will focus on the series expansion of the change of discretized energy.

$$ \Delta\Pi = \left\{ \frac{\partial\Pi}{\partial\{u\}} \right\}^T \Delta\{u\} + \frac{1}{2} \Delta\{u\}^T \frac{\partial^2\Pi}{\partial\{u\}^2} \Delta\{u\} + O(\Delta\{u\})^3 \qquad (24) $$

$$ = \{g\}^T \Delta\{u\} + \frac{1}{2} \Delta\{u\}^T [K] \Delta\{u\} + O(\Delta\{u\})^3 \qquad (25) $$

$\{g\} = \frac{\partial\Pi}{\partial\{u\}}$ is the gradient of the total potential energy $\Pi$ while $[K] = \frac{\partial^2\Pi}{\partial\{u\}^2}$ is the Hessian matrix.

The first order minimization algorithm used here is a classical conjugate gradient algorithm proposed by Fletcher and Reeves in [15]. It is based on the descent method. Nevertheless, we insist on the fact that the gradient of the total potential, which gives the opposite direction of descent, must be exactly computed. Thus, it is necessary to provide the analytical expression of the gradient of the potential in order to obtain the sufficient accuracy and to correctly handle some phenomena like bifurcation.

The first order methods were used to seek for a minimum by successive line search. The descent direction is computed thanks to the gradient $\{g\}$ which gives the direction of greatest increase of $\Pi$. The global convergence of the conjugate gradient algorithm is demonstrated in the case of convex potentials. In wrinkling problems, the potential is nonconvex and the algorithm converges toward an existing minimum which may be either a local or global minimum. Further details about this subject will be found in [16].

The main asset of first order methods is that they allow to converge towards a local or global minimum even in the neighborhood of a bifurcation point. It is not the same for second order methods which require to compute the stiffness matrix $[K]$; it may be ill-conditioned when a significant loss of stiffness occurs. At a bifurcation point, the stiffness matrix is singular and ill-conditioned in their neighborhood. This is why corrective solution methods, like the Newton-Raphson method, may run into difficulties at or near critical points. This ill-conditioning may introduce noise which makes the solution procedure unstable. Sometimes traversing critical points may become computationally overwhelming and may

require either specialized techniques or intensive human interventions, such as the introduction of geometrical imperfections or the dissipation of a small part of the energy (see [1]).

In the present paper, we used the conjugate gradient algorithm to obtain the equilibrium path. This method uses only the gradient vector of the total potential energy. The numerical developments are implemented in the Surface Evolver code developed by K. Brakke and presented in [17].

#### 4.1. Energy variation along equilibrium path

Suppose that we have a one-parameter load conservative system. The total potential energy of the discrete model can be written as:

$$ \Pi(\{u\}, \lambda) = \Pi^{\text{int}} - \lambda\{q\}\{u\} \qquad (26) $$

$\{q\}$ is the reference load vector and $\lambda$ a scalar load multiplier. Then consider the power series expansion of the energy variation at an equilibrium point $(\{u\}, \lambda)$:

$$ \begin{aligned} \Pi(\{u\} + \Delta\{u\}, \lambda + \Delta\lambda) &= \Pi(\{u\}, \lambda) + \left\{ \frac{\partial\Pi}{\partial\{u\}} \right\}^T \Delta\{u\} + \frac{\partial\Pi}{\partial\lambda} \Delta\lambda \\ & \quad + \frac{1}{2} \Delta\{u\}^T \frac{\partial^2\Pi}{\partial\{u\}^2} \Delta\{u\} + \frac{1}{2} \frac{\partial^2\Pi}{\partial\lambda^2} (\Delta\lambda)^2 \\ & \quad + \frac{\partial^2\Pi}{\partial\{u\}\partial\lambda} \Delta\{u\}\Delta\lambda + O(\Delta\lambda^3, \Delta\{u\}^3) \\ &= \Pi(\{u\}, \lambda) + \{g\}^T \Delta\{u\} - \{q\}^T \{u\} \Delta\lambda \\ & \quad + \frac{1}{2} \Delta\{u\}^T [K] \Delta\{u\} \\ & \quad - \{q\}^T \Delta\{u\} \Delta\lambda + O(\Delta\lambda^3, \Delta\{u\}^3) \end{aligned} \qquad (27) $$

If $(\{u\}, \lambda)$ is a state of equilibrium, then $\{g\} = 0$. When the change $(\Delta\{u\}, \Delta\lambda)$ occurs along the equilibrium path, the displacement $\Delta\{u\}$ and the load factor $\Delta\lambda$ are bound by the following relation:

$$ \Delta\{u\} = \Delta\lambda[K]^{-1}\{q\} \qquad (28) $$

We note $\theta_i$ the eigenvalues of $[K]$ in ascending order and $\{z_i\}$ the corresponding eigenvectors. The components of displacement vector in the eigenvectors base are noted $A_i = \Delta\{u\}^T\{z_i\}$. The Eq. (28) is rewritten in the eigenvectors base as follows:

$$ \Delta\{u\}^T \{z_i\} = \Delta\lambda \frac{1}{\theta_i} \{q\}^T \{z_i\} \quad \Rightarrow \quad A_i = \frac{\Delta\lambda}{\theta_i} q_i \qquad (29) $$

From the general expression (27), the energy variation due to $\Delta\{u\}$, at a fixed load level $\lambda + \Delta\lambda$ can be simply obtained:

$$ \begin{aligned} \Pi(\{u\} + \Delta\{u\}, \lambda + \Delta\lambda) &= \Pi(\{u\}, \lambda + \Delta\lambda) \\ &= \frac{1}{2} \Delta\{u\}^T [K] \Delta\{u\} - \{q\}^T \Delta\{u\} \Delta\lambda + O(\Delta\lambda^3, \Delta\{u\}^3) \\ &= \frac{1}{2} [\theta_1 A_1^2 + \theta_2 A_2^2 + \dots + \theta_n A_n^2] - \Delta\lambda \{q\}^T \Delta\{u\} \end{aligned} \qquad (30) $$

The component $q_i$ of vector $\{q\}$ in the base of eigenvectors of $[K]$, is given by $\Delta\lambda\{q\}\{z_i\} = A_i\theta_i$. The product $\Delta\lambda\{q\}^T\Delta\{u\}$ is then equal to $\theta_1 A_1^2 + \theta_2 A_2^2 + \dots + \theta_n A_n^2$.

Then the energy variation can be expressed either by the eigenvalues of the stiffness matrix:

$$ \begin{aligned} \Pi(\{u\} + \Delta\{u\}, \lambda + \Delta\lambda) &= \Pi(\{u\}, \lambda + \Delta\lambda) \\ &= -\frac{1}{2} [\theta_1 A_1^2 + \theta_2 A_2^2 + \dots + \theta_n A_n^2] \end{aligned} \qquad (31) $$

Or as the dot product of load vector $\{q\}$ and displacement vector $\Delta\{u\}$

$$ \Pi(\{u\} + \Delta\{u\}, \lambda + \Delta\lambda) - \Pi(\{u\}, \lambda + \Delta\lambda) = -\frac{1}{2} \Delta\lambda \{q\}^T \Delta\{u\} \qquad (32) $$
As we can see from (31), the energy variation is negative for all displacement vectors $\Delta\{u\}$ obeying to the equilibrium Eq. (28).

## 4.2. Use of downhill methods

The expression (31) shows clearly that, when the loading parameter is incremented from $\lambda$ up to level $\lambda + \Delta\lambda$, the displacement vector $\Delta\{u\}$ accompanying the load increment will decrease the total potential energy for all possible stable equilibrium paths ($\lambda_i > 0, \forall i$). In these situation, $[K]$ is positive definite and a minimum energy exists. Algorithms based on descent method converge toward the minimum.

At critical point, the minimal eigenvalue of $[K]$ is equal to zero. The energy shows a neutral state for displacements proportional to the associated eigenvector: $\Delta\{u\} = A_1\{z_1\}$ i.e. $\{z_1\}^T[K]\{z_1\} = 0$. In these circumstances, a full investigation of the stability will require higher-order terms of the energy variation along the branch of bifurcation. The uses of downhill methods fail to switch on bifurcated branch since the gradient of the energy is orthogonal to the first eigenvector.

However this type of algorithms, based on the gradient of energy, has a linear convergence. The iterative process tends towards the equilibrium path without ever reaching it exactly, because the gradient tends to zero close to the equilibrium path. This lack of convergence becomes a benefit for the treatment of bifurcation as it acts as a permanent geometry perturbation during the load incrementation process.

In some situations, we need to perturb randomly the geometry in such a way that the underlying regularity intrinsic to a bifurcation point is destroyed. This is the case in this study when we introduce a very slight random perturbation of the actual geometry at the beginning of the iterative process. However the reference configuration is not perturbed. We fixed the position perturbation to $\alpha \times 10^{-3}$ mm, where $\alpha$ is a positive random number $\alpha \in [0, 1]$. We insist on the fact that the introduced perturbation is not an imperfection on the initial geometry but just a slight deviation of the actual configuration from the equilibrium path before loading, i.e. at time $t = 0$.

# 5. Finite element discretization

This section is concerned with the discretization of the shell potential given in (17) which is the expression used to find the possible wrinkle shape of thin film structures.

Because of the small thickness of the studied structures (between 25 µm and 125 µm), the strain energy of a thin shell finite element due to bending is negligible compared to the strain energy due to tension or compression. Consequently the shell element behavior can present a membrane locking mechanism that may affect the result of the wrinkling simulation. Even though the bending stiffness is very small in thin shells, it is fundamental in wrinkling deformation since, after buckling, it highly influences the shape of the wrinkles. Indeed, the bending stiffness, even if it is very low, is the parameter which, in conjunction with in-plane tension, determines the shape of the wrinkles in terms of amplitude and wavelength. The choice of a finite element is therefore an important aspect in wrinkling simulation and we will briefly discuss it.

The finite elements used to predict the shape and the size of wrinkles are classically thin shell finite elements formed by superimposition of a genuine membrane element and a plate element even if some authors computed it only with a genuine membrane element [18,8].

In [1] various finite elements have been tested. These studies emphasize that one should use a dense mesh of primitive elements

rather than a coarse mesh made up of higher order element to accurately predict the geometric nonlinearity of very thin structures.

In the study [1] the finite element S4R5 (code ABAQUS ® ) has been chosen for being the most adapted element in terms of computational efficiency and calculation time. This quadrilateral element is formed by the superimposition of a genuine membrane element and a plate bending element. The plate bending formulation DKQ used in the S4R derived from the DKT formulation proposed by Batoz [19]. The same element was also used in [12,2].

Numerical results obtained by using a membrane-DKT triangular element are presented in this paper. In this formulation, the plate bending behavior requires the introduction of rotational degrees of freedom, which may cause troubles with ill-conditioned stiffness matrices, such as numerical locking, when the thickness is low. In order to avoid this phenomenon, a rotation free shell element has been used. These elements include the bending behavior of thin shells without introducing any additional degree of freedom compared to a membrane element. The bending energy is discretized by using the out-of-plane displacements of a patch of four triangular elements. These kinds of elements have been successfully used in thin shells buckling problems like the simulation of an air-bag inflation [20] or in metal forming simulations [21].

In this study, we have used a triangular three nodes thin shell finite element formed by the superposition of a genuine membrane element and a plate element. For the plate element, two models have been tested: the DKT and the rotation-free element proposed by Sabourin and Brunet [22].

In Section 5.1 the expression of the membrane energy $\Pi^M$ is discretized and in Section 5.2 the expression of the bending energy $\Pi^B$ is detailed for Rotation Free element and briefly depicted for the DKT element. Then, the shell energy is obtained by summing the membrane energy and the bending one.

## 5.1. Discretization of the membrane energy

As shown in Fig. 2, the tangent basis vectors on the mid-surface are computed with the help of reference $X_i^j$ and actual $x_i^j$ vertices coordinates of triangle finite elements such as:

$$ \begin{aligned} \overrightarrow{A_1} &= \overrightarrow{X_2} - \overrightarrow{X_1}, \quad \overrightarrow{a_1} = \overrightarrow{x_2} - \overrightarrow{x_1} \\ \overrightarrow{A_2} &= \overrightarrow{X_3} - \overrightarrow{X_1}, \quad \overrightarrow{a_2} = \overrightarrow{x_3} - \overrightarrow{x_1} \end{aligned} \qquad (33) $$

The components of these vectors are then:

$$ \begin{aligned} a_1^j &= x_2^j - x_1^j \\ a_2^j &= x_3^j - x_1^j \end{aligned} \qquad (34) $$

<figure><img src="image_3.png" /><figcaption>Fig. 2. Finite element discretization of the triangular membrane element.</figcaption></figure>
where the subscripts indicate the node number (i = 1, ..., 3) and the superscript refers to the coordinate on the direction $\vec{e}_j$ with (j = 1, ..., 3).

The mixed components of the Green–Lagrange strain tensor are then:

$$E_{\alpha}^{\beta} = \chi_{\alpha}^{\beta} = \frac{1}{2} (\overrightarrow{a}_{\alpha} \cdot \overrightarrow{a}_{\beta} \overrightarrow{A}^{\alpha} \cdot \overrightarrow{A}^{\beta} - \delta_{\alpha}^{\beta}) = \frac{1}{2} (C_{\alpha}^{\beta} - \delta_{\alpha}^{\beta}) \quad (35)$$

The membrane energy $\psi^M$ density is evaluated, with plane stress assumption, by using the expression given in Eq. (12) which leads to:

$$\psi^M = \frac{\mu \lambda}{2(\lambda + 2\mu)} [3 - 3C_1^1 + (C_1^1)^2 + C_1^1 C_2^2 + C_1^2 C_2^1 - 3C_2^2 + (C_2^2)^2] \\ + \frac{\mu^2}{2(\lambda + 2\mu)} [2 - 2C_1^1 + (C_1^1)^2 + 2C_1^2 C_2^1 - 2C_2^2 + (C_2^2)^2] \quad (36)$$

Since the element is linear, the strain tensor is constant. The energy of a finite element (e) is simply computed by the product of the energy density by the volume of the element.

$$\Pi_e^M = Ah\psi_e^M \quad (37)$$

where A is the area of the finite element (e) and h its thickness.

The nine components of the finite element strain energy gradient are the partial derivatives of $\psi^m$ with respect to the coordinates of the vertices $x_i^j$ i.e.

$$(\nabla \Pi_e^M)_{ij} = Ah \frac{\partial \psi_e^M}{\partial x_i^j} = Ah \frac{\partial \psi_e^M}{\partial C_\alpha^\beta} \frac{\partial C_\alpha^\beta}{\partial x_i^j} \quad (38)$$

This is the exact formulation of the gradient of the membrane strain energy used for large scale problems.

## 5.2. Discretization of the bending energy

The rotation free triangular plate bending element used here was initially presented in [22]. Here we follow the formulation proposed in [23].

A flat patch of elements formed of a central triangular element and its three neighbors is considered (see Fig. 3). The Kirchhoff theory is used to approximate the curvatures tensor $\overline{\rho}$ from the out-of-plane displacements of the patch of four triangular elements. $\overline{n}_i$, ($i$ = 1, ..., 3) refers to the normal of the side i and $h_i$ are the orthogonal distance of vertex i, ($i$ = 1, ..., 6) to the opposite edge of the central triangle.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Discretization of a plate triangular patch element.</figcaption></figure>

The relative vertical displacements $w_i^*$ are also the out of plane displacements that cause the bending of the surface. The curvature tensor is assumed constant and consequently the surface is assumed parabolic irrespective of a relative vertical displacement $w^*$. The out of plane displacement along the normal $\overline{n}_1$, is then:

$$w^* = \frac{n(n - h_4)}{h_1(h_1 + h_4)} w_1^* + \frac{n(h_1 + n)}{h_4(h_1 + h_4)} w_4^* \quad (39)$$

where n is the coordinate along the normal direction. The curvature along the direction $\overline{n}_1$ is the second derivative of the out of plane displacement $w^*$ with respect to the coordinate n, i.e.

$$\rho \overline{n}_1 = \frac{d^2 w^*}{dn^2} = \frac{2}{(h_1 + h_4)} (\alpha_1 + \alpha_4) \quad (40)$$

where $\alpha_i = \frac{w_i^*}{h_i}$ are the rotation angles between two adjacent elements. The rotation angles are considered as rigid-body rotations. Then identical relations can be obtained for sides 2 and 3.

The curvatures along the directions $\overline{n}_1$ are used to obtain the components $\rho_{\alpha\beta}$ of the curvature tensor of the main element. Then it is possible to interpolate the rotation angle $\alpha_i$ from the patch element out of plane displacements vector {$w$}. It follows:

$$\{\rho\} = \begin{bmatrix} \rho_{11} \\ \rho_{22} \\ \rho_{12} \end{bmatrix} = [B^B] \cdot \{w\} \quad (41)$$

Further details about the formulation of the matrix $[B^B]$ may be found in [23].

As for the membrane strain tensor, the curvature tensor is constant over the central element of the patch. Thus the bending stiffness matrix of a particular finite element (e) is then written as:

$$[K_e^B] = A \cdot [B^B]^T \cdot \frac{h^2}{12} \overline{\overline{H}} \cdot [B^B] \quad (42)$$

where A is the area of the central element. The bending energy is a quadratic expression of the nodal displacement:

$$\Pi_e^B = \frac{1}{2} \cdot \{w\}^T \cdot [K_e^B] \cdot \{w\} \quad (43)$$

From which we deduce the gradient of $\Pi^B$:

$$\nabla \Pi_e^B = [K_e^B] \cdot \{w\} \quad (44)$$

The DKT formulation leads to similar expression for the energy and its gradient as in Eqs. (43) and (44). The vectors of degree of

<figure><img src="image_5.png" /></figure>
freedom and the interpolation matrix $[B^B]$ used in the patch element and DKT element are different. In the patch formulation, each node of the patch possesses only three translations as degree of freedom, while the DKT formulation uses three translations and three rotations for each node. The introduction of this rotational degree of freedom can possibly lead to ill-conditioned matrix.

The thin shell strain energy, for each element (e) of the wole structure, is the sum of the membrane strain energy and the bending one i.e.

$$ \Pi_e^{\text{int}}(\{u\}) = \Pi_e^{\text{M}}(\{u\}) + \Pi_e^{\text{B}}(\{u\}) \quad (45) $$

The gradient of the strain energy is then:

$$ \nabla \Pi_e^{\text{int}}(\{u\}) = \nabla \Pi_e^{\text{M}}(\{u\}) + \nabla \Pi_e^{\text{B}}(\{u\}) \quad (46) $$

The total internal energy of the structure is:

$$ \Pi^{\text{int}}(\{u\}) = \sum_{e} \Pi_{e}^{\text{int}}(\{u\}) \quad (47) $$

which is the discrete form of the internal energy given in Eq. (21) while the gradient is:

$$ \nabla \Pi^{\text{int}}(\{u\}) = \sum_{e} \nabla \Pi_{e}^{\text{int}}(\{u\}) \quad (48) $$

## 6. Numerical samples

This section shows the ability of a first order minimization algorithm to predict a possible wrinkle shape. The performance of the method proposed here is demonstrated by computing the wrinkles details of three thin film structures for which experimental data are available:

* The wrinkling prediction on a rectangular shell-membrane under transverse in-plane displacement.

* The wrinkling prediction on a square membrane under corner loads.

* The study of a cruciform membrane under biaxial load.

In these three problems, the wrinkling is the result of particular loads, membrane shapes or boundary conditions. Nevertheless, wrinkling may occurs even in the common case of a simple rectangular membrane longitudinally stretched. (see [24,25]). This is a consequence of the clamped boundary conditions preventing lateral displacements along the loaded edges as explained in [24].

The wrinkling problem of the cruciform structure has been investigated experimentally in [2]. In the same papers the experimental results have been compared with numerical wrinkling simulations performed using the post-buckling analysis presented in [1] (code ABAQUS ®). This simulation procedure starts with an eigenvalue buckling analysis which allows to obtain the possible wrinkling

modes of the membrane. The eigenvectors of the tangent stiffness matrix are possible wrinkling shapes of the membrane. Then, mode shapes are introduced as initial geometric imperfections in the structure before performing the post-buckling computation. It is carried out using the Newton Raphson method. Nevertheless when the formation of an additional wrinkle occurs, the equilibrium path may be unstable. To achieve convergence in spite of the presence of unstable paths, pseudo-viscous forces are added to the model. These forces are scaled by a stabilization factor which corresponds to a fraction of the dissipated strain energy. The stabilization parameter could affect the shape of the wrinkled membrane. Consequently, to achieve a good accuracy, it is generally preferable to set this parameter to the lowest possible value for which convergence can still be obtained (see [1]).

The main advantages of the post-buckling analysis are as follows:

It implies numerical methods usually available on the majority of finite elements codes. It allows to focus on the response of the structure on different stability paths.

The main drawback of the post-buckling analysis is it using complexity, because it requires to choose the wrinkling modes to introduce as geometrical imperfections and an appropriate value for the prescribed fraction of the dissipated strain energy. Moreover, crossing a critical point may still be difficult and some restarts of the analysis are often needed.

The first order minimization method, used to perform the three following wrinkling simulations instead of the post-buckling analysis, avoids these numerical difficulties. The conjugate gradient algorithm is able to converge towards a local minimum even in the neighborhood of a bifurcation point. Consequently, additional numerical techniques are not necessary to achieve the convergence of the computation. To obtain a bifurcated solution, a small random displacement has been applied at each node of the membrane surface. This allows to leave the fundamental equilibrium path. Then a bifurcated solution is automatically reached but a particular equilibrium path cannot be chosen to perform the analysis with this procedure.

### 6.1. The membrane in simple shear

Firstly, we present the study of the initially flat and stress free linear elastic membrane proposed in [26]. The dimensions of this structure are depicted in Fig. 4. The material that was used is the Kapton HN ®. Its material properties are listed in Table 1.

<table><caption>**Table 1**
Shear test of a membrane: material properties.</caption><tr><td>Thickness (μm)</td><td>25</td></tr><tr><td>Young's modulus E (MPa)</td><td>3500</td></tr><tr><td>Poisson's ratio ν</td><td>0.31</td></tr></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Shear test of a membrane: geometry.</figcaption></figure>
A shear displacement of 3 mm is imposed on the upper edge. The horizontal edges are held at fixed distance while the two side's edges are unconstrained.

The finite element mesh used here consists of 26,660 triangular thin shell elements. Two different formulations have been tested: the membrane-DKT element and the membrane-patch element. Fig. 5 shows the wrinkle pattern of the membrane, computed by using the membrane-patch element, for a shear displacement $\delta_c = 3$ mm.

To achieve the simulation convergence, a small random displacement of magnitude $\alpha \times 10^{-3}$ mm has been applied at each node of the membrane surface before performing the shear displacement. This allows to leave the fundamental equilibrium path.

Fig. 5 shows that the wrinkle shape presents two slack regions near the unconstrained edges, while the wrinkles in the central regions are inclined at 45° to the upper and bottom edges of the structure. The wrinkle pattern is the same as that experimentally observed in [26].

The wrinkles shapes are further investigated by plotting the central cross section ($y = 64$ mm) as shown in Fig. 6. First, it appears that the simulations performed with the DKT element and the patch element lead to different results even if they are close. The wrinkle shape obtained using the DKT element is slightly noisy and one wrinkle starts to divide in two new wrinkles. It is probably a problem of membrane locking due to the very small thickness of the membrane.

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Wrinkle pattern of the membrane in simple shear computed using the membrane-patch element ($\delta_c = 3$ mm).</figcaption></figure>

The cross section shows 18 wrinkles for the DKT element against 19 with the patch element. That is the number of wrinkles observed by Wong and Pellegrino in [26]. Moreover, the wavelength $\lambda$ predicted by the numerical simulation is 16.6 mm and the average amplitude $A_w$ is 0.38 mm which is close to the experimental observations ($\lambda \approx 20$ mm and $A_w \approx 0.4$ mm).

## 6.2. The square membrane under corner load

The second validation test is the wrinkling prediction on a square membrane under corner load proposed in [26]. The membrane is loaded at the four corners by two diagonal pairs of equal and opposite forces. It was loaded up to 5 N at the two horizontal corners and up to 20 N at the two vertical corners. The geometry of this structure is depicted in Fig. 7. The material that was used is the same that for the previous sample, its material properties are listed in Table 1.

The corner of the structure were cut at 45° to a width of 25 mm. They are reinforced with Kapton® adhesive and steel pins over an area of 25 mm $\times$ 20 mm. Here, we have modeled these parts as rigid body. The finite element mesh used here consists of 62,300 triangular thin shell element. As for the previous sample, simulations have been performed using the membrane-DKT element and the membrane-patch element.

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** The square membrane under corner load: geometry and load case.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Central cross section of the membrane in simple shear ($\delta_c = 3$ mm).</figcaption></figure>
The numerical wrinkle shape (see Fig. 8) shows a large vertical wrinkle as observed in [26]. Here the use of the membrane-DKT or the membrane-patch element leads to very similar results. Wong and Pellegrino have observed a nonsymmetric shape for the large vertical wrinkle with a maximum of approximately 2.6 mm and a minimum of approximately -1.2 mm. The numerical simulation has predicted a maximum of 2.76 mm with the membrane-DKT element and 2.65 mm with the membrane-patch element. The minimum observed are respectively -0.96 mm with the membrane-DKT element and -1.066 mm with the membrane-patch element.

For the two previous samples, the numerical predictions show good agreement with the experimental data either using the membrane-DKT element or the membrane-patch element. Nevertheless, best results have been obtained with the membrane-patch element. Thus for the last sample, only this element has been used.

### 6.3. The cruciform membrane under biaxial load

The last numerical sample is the wrinkling prediction on a cruciform membrane under biaxial load. In [2] an experimental study of the formation and evolution of the wrinkle pattern that form in flat elastic and isotropic membranes under the action of in-plane tension is presented. The experiments were carried out on a cruciform specimen stretched along two uncoupled axes using various loading paths. The wrinkled shapes of the membrane were digitized by using a full-field measurement based on the fringe analysis method. This set of experiment allows us to discuss the quality of the numerical results for different thicknesses of the membrane, and the reproducibility of a kinematic configuration of wrinkles.

The biaxial tests were performed on cruciform specimens according to the geometry depicted in Fig. 9.

The membrane tested is the Kapton VN® polyimide film. The experimental Young's modulus and the Poisson's ratio corresponding to the experimental conditions (the average strain imposed upon the structure is less than 1%) have been identified in [2]. Table 2 summarizes the material characteristics of the specimens.

Various sets of displacements have been applied to the horizontal and vertical sides of the specimen. The displacement $δ_1$ is a negative displacement performed along the axis X while the displacement $δ_2$ is a tensile displacement imposed along the axis Y. Before performing the displacement ($δ_1$ and $δ_2$) the horizontal and vertical side edges are clamped while the circular edges are

<figure><img src="image_10.png" /><figcaption>Fig. 9. The cruciform membrane under biaxial load: geometry.</figcaption></figure>

<caption>**Table 2**</caption>

The cruciform membrane under biaxial load: material properties.

<table><tr><td>Thickness (μm)</td><td>25, 50 and 125</td></tr><tr><td>Young's modulus E (MPa)</td><td>3350</td></tr><tr><td>Poisson's ratio ν</td><td>0.3</td></tr></table>

unrestrained. We have computed here the wrinkle shape of 125 µm, 50 µm, and 25 µm, membranes for the load case $δ_1 = -3$ mm and $δ_2 = 3$ mm where the negative displacement $δ_1$ is performed before the tensile displacement $δ_2$.

Figs. 10–12 show the wrinkle shape obtained for respectively 125 µm, 50 µm and 25 µm membranes. The numerical results have been obtained using a mesh of 64,000 triangular membrane-patch elements.

The wrinkles shapes are further investigated by plotting the central cross section y = 0 mm and by superimposing the experimental and numerical sections on the same figure. For the speci-

<figure><img src="image_11.png" /><figcaption>Fig. 8. Wrinkle shape on the square membrane under corner load.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 10. Numerical out of plane displacements (Specimen thickness = 125 µm).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. Numerical out of plane displacements (Specimen thickness = 50 µm).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12. Numerical out of plane displacements (Specimen thickness = 25 µm).</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 13. Comparison of the numerical and experimental wrinkling results for $\delta_1 = -2$ mm and $\delta_2 = 0.5$ mm (Specimen thickness = 125 $\mu$m).</figcaption></figure>

men of 125 $\mu$m the shape of wrinkles (pattern, wavelength and amplitude) predicted by the numerical simulation appear to be close to those observed during the experiment. Nevertheless, for the 50 $\mu$m and 25 $\mu$m specimens, noticeable differences have been observed concerning the amplitude of the wrinkles in the central area. For instance, with the specimen of 25 $\mu$m the experimental value of the central wrinkle amplitude (at $x=0$) is about 0.27 mm against 0.18 mm for the numerical one.

As shown in Figs. 10–12 the accuracy of the numerical predictions decreases when the thickness reduces. Two phenomena could explain the differences between the experimental and the numerical data. The first one is the finite element discretization. The shell elements used here are probably sensitive to numerical locking when the thickness is very low. The second one concerns experimental considerations. Before performing the experiment on the biaxial set up, the cruciform specimen depicted in Fig. 9 is slightly stretched and then assumed flat (see [2]). In the experiment, the initial shape of the membrane may show a small deflection which could affect the pattern of wrinkles when performing the displacements.

As explained in Section 4.2, the conjugate gradient algorithm can converge toward different local minima if they exist. Thanks to this phenomenon, different possible wrinkle shapes of the structures can be obtained for a given load parameter. These different wrinkle shapes are associated with different buckling modes of the structure. Here, for the load case $\delta_1 = -2$ mm and $\delta_2 = 0.5$ mm we have obtained different wrinkle shapes by applying different random out of plane displacements $\alpha$ at the start of the simulation (see Fig. 13). The different solutions are associated with possible kinematic configurations of the physical system. Nevertheless the algorithm doesn't allow choosing a particular bifurcated branch to achieve convergence.

## 7. Conclusions

The wrinkling simulations have been performed by directly minimizing the total potential energy using a first order descent method, instead of satisfying the equilibrium equations as done with the usual finite element method. It proves to be an attractive alternative which is particularly efficient for thin wrinkled structures. Indeed, the proposed solving procedure has been first validated through three problems: the wrinkling prediction of rectangular membrane under transverse in-plane displacements

and the wrinkling prediction in a square membrane under corner loads. Then comparison has been done with a biaxial experiment performing on cruciform membrane specimens. All the wrinkles patterns predicted during the numerical computations showed a good agreement with the experimental results.

Moreover the main interest of the minimization algorithm is that no specific buckling analysis is carried out during the calculation. The wrinkling procedure is simply started by applying a small random displacement at each node of the membrane surface.

## Acknowledgments

Financial support from CNES and EADS Astrium is gratefully acknowledged. The author thanks C. Talec for useful corrections and K. Brakke for his helpful code "Surface Evolver".

## References

[1] Wong YW, Pellegrino S. Wrinkled membranes Part III: Numerical simulations. J Mech Mater Struct 2006;1:61–93.

[2] Lecieux Y, Bouzidi R. Experimental analysis on membrane wrinkling under biaxial load – comparison with bifurcation analysis. Int J Solids Struct 2010;47(18–19):2459–75. doi:10.1016/j.ijsolstr.2010.05.005.

[3] Wagner H. Flat sheet metal girders with very thin metal web, Tech. Rep. 604, National Advisory Committee for Aeronautics; 1929.

[4] Kang S, Im S. Finite element analysis of wrinkling membranes. J Appl Mech 1997;64:263–9.

[5] Epstein M, Forcinito MA. Anisotropic membrane wrinkling: theory and analysis. Int J Solids Struct 2001;38(30–31):5253–72.

[6] Mansfield EH. Tension field theory: a new approach which shows its duality with inextensional theory, in: XII Int Cong Appl Mech; 1968, p. 305–20.

[7] Pipkin A. The relaxed energy density for isotropic elastic membranes. IMA J Appl Math 1986;36:8599.

[8] Diaby A, van AL, Wielgosz C. Buckling and wrinkling of prestressed membranes. Finite Elem Anal Des 2006;42(11):992–1001. doi:10.1016/j.finel.2006.03.003.

[9] Riks E. An incremental approach to the solution of snapping and buckling problems. Int J Solids Struct 1979;15(7):529–51.

[10] Crisfield M. Non-linear finite element analysis of solids and structures: advanced topics. New York, NY, USA: John Wiley & Sons, Inc; 1997.

[11] Magnusson A. Treatment of bifurcation points with asymptotic expansion. Comput Struct 2000;77(5):475–84. doi:10.1016/S0045-7949(00)00004-3.

[12] Wang C, Tan H, Du X, Wan Z. Wrinkling prediction of rectangular shell-membrane under transverse in-plane displacement. Int J Solids Struct 2007;44(20):6507–16. doi:10.1016/j.ijsolstr.2007.02.036.

[13] Simo JC, Fox DD. On a stress resultant geometrically exact shell model. Part I: Formulation and optimal parametrization. Comput Methods Appl Mech Eng 1989;72(3):267–304. doi:10.1016/0045-7825(89)90002-9.

[14] Simo JC, Fox DD, Rifai MS. On a stress resultant geometrically exact shell model. Part II: The linear theory; computational aspects. Comput Methods Appl Mech Eng 1989;73(1):53–92. doi:10.1016/0045-7825(89)90098-4.
[15] Fletcher R, Reeves CM. Function minimization by conjugate gradients. The Comput J 1964;7(2):149-54. doi:10.1093/comjnl/7.2.149. arXiv: http://comjnl.oxfordjournals.org/cgi/reprint/7/2/149.pdf.

[16] Grippo L, Lucidi S. A globally convergent version of the polak-ribière conjugate gradient method. Math Program 1997;78(3):375-91. doi: http://dx.doi.org/10.1007/BF02614362.

[17] Brakke KA. The surface evolver. Exp Math 1992;1:141-65.

[18] Tabarrok B, Qin Z. Nonlinear analysis of tension structures. Comput Struct 1992;45(5-6):973-84. doi:10.1016/0045-7949(92)90056-6.

[19] Batoz J, Dhatt G. Modélisation des structures par éléments finis, vol. 2 poutres et plaques, HermFs; 1990.

[20] Flores FG, Oate E. Improvements in the membrane behaviour of the three node rotation-free bst shell triangle using an assumed strain approach. Comput Methods Appl Mech Eng 2005;194(6-8):907-32. doi:10.1016/j.cma.2003.08.012.

[21] Brunet M, Sabourin F. A simplified triangular shell element with a necking criterion for 3-d sheet-forming analysis. J Mater Process Technol 1995;50(1-4):238-51. doi:10.1016/0924-0136(94)01384-D.

2nd Int Conf Numer Simulat
3-D Sheet Metal Form Proc.

[22] Sabourin F, Brunet M. Analysis of plates and shells with a simplified three node triangular element. Thin-Walled Struct 1995;21(3):209-23. doi:10.1016/0263-8231(94)00001-G. Finite Elem Thin-Walled Struct.

[23] Gärdsback M, Tibert G. A comparison of rotation-free triangular shell elements for unstructured meshes. Comput Methods Appl Mech Eng 2007;196(49-52):5001-15. doi:10.1016/j.cma.2007.06.017.

[24] Friedl N, Rammerstorfer FG, Fischer FD. Buckling of stretched strips. Comput Struct 2000;78(1):185-90. doi:10.1016/S0045-7949(00)00072-9. <http://www.ingentaconnect.com/content/els/00457949/2000/00000078/00000001/art00072>.

[25] Cerda E, Ravi-Chandar, Mahadevan LK. Thin films: wrinkling of an elastic sheet under tension. Nature 2002;419:579-80.

[26] Wong YW, Pellegrino S. Wrinkled membranes Part I: Experiments. J Mech Mater Struct 2006;1:1-23.