# Buckling and wrinkling of prestressed membranes

Adama Diaby, Anh Le Van, Christian Wielgosz

► To cite this version:

Adama Diaby, Anh Le Van, Christian Wielgosz. Buckling and wrinkling of prestressed membranes. Finite Elements in Analysis and Design, 2006, 42 (11), pp.992 - 1001. 10.1016/j.finel.2006.03.003. hal-01005295

HAL Id: hal-01005295

https://hal.science/hal-01005295

Submitted on 22 Oct 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
# Buckling and wrinkling of prestressed membranes

Adama Diaby, Anh Le van, Christian Wielgosz

GeM (Research Institute in Mechanics and Civil Engineering), Faculty of Sciences at Nantes, 2, rue de la Houssiniere, BP 92208, Nantes 44322 cedex 3, France

This paper deals with the numerical computation of buckles and wrinkles appearing in membrane structures by means of the total Lagrangian formulation, using genuine membrane finite elements (with zero bending stiffness) and a prestressed hyperelastic constitutive law. The bifurcation analysis is carried out without assuming any imperfections in the structure. The standard arclength method is modified by means of a specific solution procedure to cope with the occurrence of complex roots when solving the quadratic constraint equation. Applying the proposed formulation to a set of typical numerical examples shows its ability to correctly predict the wrinkling and buckling behaviour in membrane structures.

**Keywords:** Finite elements; Membranes; Hyperelasticity; Bifurcation; Buckling; Wrinkling

## 1. Introduction

The wrinkling of nonlinear elastic membranes has been investigated by many authors using a variety of approaches. Among those, two solution procedures are widely used: the bifurcation analysis and the tension field theory.

The tension field theory is generally applied on membrane elements, which have zero bending stiffness. In conjunction with wrinkling material models, the state (wrinkled, slack or taut) of each membrane element is assessed using different wrinkling criteria and the material properties are adjusted iteratively during the analysis to account for the behavior associated with the particular state of the element in hand. Wagner in 1929 [1] was the first to apply this theory in order to estimate the maximum shear load that can be carried by a thin web. In the early 1960's, Stein and Hedgepeth [2] predicted the stresses and strains in stretched membranes for both the wrinkled and unwrinkled regions, by employing the usual theory of elasticity. They dealt with the wrinkling of the membrane by modifying the constitutive equation of the material so that the membrane remains in a state of simple tension throughout the wrinkled region.

A wrinkling region is created when one of the two in-plane principal stresses becomes negative. Simple formulations and extensions of this theory were later proposed in Refs. [3–11]. The tension field theory can also be addressed from an energetic point of view, as presented in Refs. [12,13]. According to this, the relaxed energy density ensures that unstable compressive stresses should never appear in the solution. Barsotti et al. [14] considered the membrane as a von Karman plate with negligible bending stiffness. Their relaxed energy concept results in a consistent linear-wrinkle-elasticity theory, which allows identifying the boundary between taut, slack and wrinkled regions. The tension field theory provides correct stress distributions and predictions for wrinkled and slack regions, yet not wrinkle details such as amplitude, wavelength and number of wrinkles.

To take into account wrinkle details, another representation of membrane wrinkling can be obtained by using the bifurcation analysis coupled with shell elements with non-zero bending stiffness. The shell-based techniques make use of thin shell elements in conjunction with geometric imperfections embedded into the initial mesh to perform a postbuckling analysis. Wong and Pelligrino [15–17] used this approach to study the wrinkling formation and evolution in a sheared membrane. In the analytical approach, they assumed that a membrane carries a uniform compressive principal stress equal to the buckling
stress of an infinitely wide thin plate. On the basis of this assumption, they showed that the wavelength and the amplitude of the wrinkles are functions of the shear angle and of the geometric and material properties of the membrane. The results were then compared to both experimental measurements and finite element simulations. Other authors [18,19] investigated the wrinkle formation in the membrane using an updated Lagrangian schemes on similar shell elements.

Modelling the thin-film membranes can be made by using the cable method. Stanuszek [20,21] created a membrane/cable finite element based on the “natural approach” and subelement technique. He showed that the resulting element was able to predict the folds in several complex examples.

Few studies investigated wrinkles in membrane structures by coupling genuine membrane elements and the bifurcation analysis. Miyumura [22] studied the wrinkles on a stretched circular membrane under in-plane torsion experimentally and numerically. He formulated a four-node isoparametric membrane element without bending stiffness and used bifurcation analysis to study the wrinkles.

In this paper, we propose a finite element analysis of the buckling/wrinkling phenomenon in membranes in the spirit of the last above-mentioned approach. Use will be made of the total Lagrangian formulation and the resulting membrane element has zero bending stiffness with a prestressed hyperelastic constitutive law. The bifurcation analysis is carried out without introducing any imperfections in the structure. Simple efficient techniques will be shown to deal with the possible initial singularity of the stiffness and the occurrence of complex roots when solving the arclength method. A number of typical numerical examples will be presented in order to assess the validity of the proposed formulation. The outlined nonlinear procedure has been developed using a home made FORTRAN 90 program.

## 2. Finite element formulation

In this section, the finite element formulation is presented within the framework of the total Lagrangian formulation. In the 3D space with a fixed Cartesian coordinate system (O; $e_1$,$e_2$,$e_3$), consider a membrane with the reference configuration defined by the middle surface $S_0$ and the thickness $h_0$. Its current configuration is defined by the middle surface S and the thickness h. The reference surface $S_0$ is discretized into isoparametric elements, the current element will be denoted e and its reference (or parent) element is denoted $e_\xi$. The reference and current positions of a material point on the membrane surface—$\mathbf{X}=X^i e_i$ and $\mathbf{x}=x^i e_i$, respectively—are parametrized by the reference coordinates $(\xi_1, \xi_2)$ as shown in the following chain: $(\xi_1, \xi_2) \in e_\xi \mapsto \mathbf{X} \in e \subset S_0 \mapsto \mathbf{x} \in S$, see Fig. 1. The mid-surface displacement field U is defined by $\mathbf{U}=\mathbf{x}-\mathbf{X}$.

The covariant base vectors on the undeformed membrane are defined by

$$ \mathbf{G}_1 = \frac{\partial \mathbf{X}}{\partial \xi_1}, \quad \mathbf{G}_2 = \frac{\partial \mathbf{X}}{\partial \xi_2}, \quad \mathbf{G}_3 = \frac{\mathbf{G}_1 \wedge \mathbf{G}_2}{\| \mathbf{G}_1 \wedge \mathbf{G}_2 \|}. \quad (1) $$

The homologous vectors on the deformed membrane are defined by

$$ \mathbf{g}_1 = \frac{\partial \mathbf{x}}{\partial \xi_1}, \quad \mathbf{g}_2 = \frac{\partial \mathbf{x}}{\partial \xi_2}, \quad \mathbf{g}_3 = \frac{\mathbf{g}_1 \wedge \mathbf{g}_2}{\| \mathbf{g}_1 \wedge \mathbf{g}_2 \|}. \quad (2) $$

The deformation gradient tensor F relates the base vectors of the reference configuration to the base vectors of the deformed configuration as

$$ \mathbf{F} = \mathbf{g}_1 \otimes \mathbf{G}^1 + \mathbf{g}_2 \otimes \mathbf{G}^2 + \lambda_3 \mathbf{g}_3 \otimes \mathbf{G}^3, \quad (3) $$

where $\lambda_3$ represents the through-thickness stretch. The strains are measured by the Green tensor E and its covariant components

$$ \mathbf{E} = \frac{1}{2}(\mathbf{F}^T \cdot \mathbf{F} - \mathbf{G}), \quad E_{\alpha\beta} = \frac{1}{2}(g_{\alpha\beta} - G_{\alpha\beta}), \\ E_{\alpha3} = 0, \quad E_{33} = \frac{1}{2}(\lambda_3^2 - 1), \quad (4) $$

where Greek indices range from 1 to 2, and $g_{\alpha\beta}$ and $G_{\alpha\beta}$ are the metric tensor components in the deformed and undeformed configurations defined by

$$ g_{\alpha\beta} = \mathbf{g}_{\alpha} \cdot \mathbf{g}_{\beta}, \quad G_{\alpha\beta} = \mathbf{G}_{\alpha} \cdot \mathbf{G}_{\beta}. \quad (5) $$

Assuming a compressible hyperelastic behaviour of the membrane, the second Piola–Kirchhoff (symmetric) stress tensor $\Sigma$ is given by differentiating the strain energy function per unit volume w with respect to the Green strain E

$$ \Sigma = \Sigma^{ij} \mathbf{G}_i \otimes \mathbf{G}_j = \frac{\partial w}{\partial \mathbf{E}} \leftrightarrow \Sigma^{ij} = \frac{\partial w}{\partial E_{ji}}, \quad (6) $$

where Latin indices take the values 1, 2 and 3. The simplest constitutive model is described by the Saint-Venant Kirchhoff strain energy

$$ w(\mathbf{E}) = \mathbf{\Sigma}_0 : \mathbf{E} + \frac{1}{2}\lambda(\mathrm{tr}\mathbf{E})^2 + \mu \mathrm{tr}(\mathbf{E}^2), \quad (7) $$

where $\lambda$ and $\mu$ are the Lamé's constants and $\mathbf{\Sigma}_0$ the prestress. In the present work, all the numerical computations are performed with a compressible neo-Hookean model defined by the following strain energy [23]

$$ w = \mathbf{\Sigma}_0 : \mathbf{E} + \frac{\lambda}{2}(\ln J)^2 - \mu \ln J + \frac{\mu}{2}(\mathrm{tr} \mathbf{C} - 3), \quad (8) $$

where $\mathbf{C}=2\mathbf{E}+\mathbf{I}$ is the right stretch tensor, $J = \sqrt{\det \mathbf{C}}$, and $\lambda$ and $\mu$ are material constants. It can be checked that in small deformations, the strain energy (8) gives the same response as the linear isotropic model with the same $\lambda$ and $\mu$ values.

Since the thickness of the membrane is small compared with other dimensions, one assumes the plane stress condition $\Sigma^{13} = \Sigma^{23} = \Sigma^{33} = 0$. With potentials (7) and (8), the conditions $\Sigma^{13} = \Sigma^{23} = 0$ are equivalent to $E_{13} = E_{23} = 0$, whereas the condition $\Sigma_{33} = 0$ entails an expression for $E_{33}$ in terms of the in-plane strains. Accordingly, the 3D constitutive law should be modified in order to obtain in-plane contravariant stresses $\Sigma^{\alpha\beta}$ as functions of in-plane covariant strains $E_{\gamma\delta}$. In the case
<figure><img src="image_2.png" /><figcaption>Fig. 1. Parametrization of initial and current configurations of the membrane.</figcaption></figure>

of potential (8), the 2D stress–strain relationship is

$$ 
\begin{aligned}
\Sigma^{\alpha\beta} &= \mu(G^{\alpha\beta} - (C^{-1})^{\alpha\beta}) \\
&\quad + \frac{\lambda}{2} \ln \left[ \frac{\lambda q}{2\mu} W \left( \frac{2\mu}{\lambda q} \exp \left( \frac{2\mu}{\lambda} \right) \right) \right] (C^{-1})^{\alpha\beta}, \\
q &\equiv \frac{\det(G_{\alpha\beta})}{\det(G_{\beta\alpha})},
\end{aligned} 
\quad (9) $$

where W is the so-called Lambert W function giving the solution of equation $u \exp(u) = z$ as $u = W(z)$.

In order to derive the discretized equilibrium equations for the membrane, use is made of the principle of virtual work: for all virtual displacement field $\mathbf{U}^*$, one has

$$ - \int_{S_0} \mathbf{\Pi} : \mathbf{grad} \, \mathbf{U}^* h_0 \, dS_0 + \int_{S_0} \mathbf{f}_0 \cdot \mathbf{U}^* \, dS_0 \\ + \int_S \mathbf{U}^* \cdot \mathbf{p} \, \mathbf{n} \, dS = 0, \quad (10) $$

where $\mathbf{\Pi} = \mathbf{F}\mathbf{\Sigma}$ is the first Piola–Kirchhoff (nonsymmetrical) stress tensor, $\mathbf{f}_0$ the dead load per unit reference surface, $p$ a possible pressure over the membrane and $\mathbf{n}$ the unit outward normal in the deformed configuration. The resulting discretized equations form a nonlinear matrix equation system with unknown $\{\mathbf{U}\}$—the nodal displacement vector of the membrane, which is solved by means of the Newton iterative scheme. By denoting NNE the node numbers in element $e$, the element stiffness matrix due to the internal force is computed by

$$ \forall i, j \in [1, 3\text{NNE}], \\ K_{ij}^{e} = \int_{e} N_{a,\alpha} N_{b,\beta} \left[ \delta_{pq} \mathbf{\Sigma}^{\alpha\beta} + x_{p,\gamma} x_{q,\delta} \frac{\partial \mathbf{\Sigma}^{\gamma\alpha}}{\partial \mathbf{E}^{\delta\beta}} \right] h_0 \, dS_0. \quad (11) $$

In the above, the integers $a, b \in [1, \text{NNE}]$ and $p, q \in [1, 3]$ are determined from indices $i, j \in [1, 3\text{NNE}]$ by $i = 3(a - 1) + p, j = 3(b - 1) + q$, and there is implicit summation over Greek indices $\alpha, \beta, \gamma, \delta = 1, 2$. The notation $N_{a,\alpha}$ means partial differentiation of the shape function $N_a$ with respect to $\xi_\alpha$.

If the membrane is subjected to a pressure $p$ which is a follower force, the following stiffness matrix is added to the

previous one: $\forall i, j \in [1, 3\text{NNE}]$,

$$ (K_{\text{pressure}}^{e})_{ij} = \frac{1}{2} \int_{e_{\xi}} p [(N_a N_{b,1} - N_b N_{a,1}) x_{p \wedge q, 2} \\ -(N_a N_{b,2} - N_b N_{a,2}) x_{p \wedge q, 1}] d\xi_1 d\xi_2 \\ + \frac{1}{2} \oint_{\partial e_{\xi}} p N_a N_b (x_{p \wedge q, 2} v_1 - x_{p \wedge q, 1} v_2) \, ds. \quad (12) $$

The integers $a, b \in [1, \text{NNE}]$ and $p, q \in [1, 3]$ are related to indices $i, j \in [1, 3\text{NNE}]$ in the same way as in (11). The symbol $x_{p \wedge q}$ means $x_{p \wedge q} = \varepsilon_{pqm} x_m$ (the integers $p$ and $q$ being fixed, there is only one value for $m$ so that $\varepsilon_{pqm}$ is not zero), $\mathbf{v}$ is the unit outward normal to the reference element $e_\xi$, at a point on its boundary $\partial e_\xi$. The surface integral over $e_\xi$ is symmetric, whereas the line integral along $\partial e_\xi$ is skew-symmetric. Let the pressure be applied over a surface portion $S_p \subset S$, the line integral contribution at a point located on the boundary of an element inside $S_p$ cancels out with that of an adjacent element, so that at last there only remains the contribution at those points located on the boundary of $S_p$.

In the following numerical examples, whenever there is a pressure over the membrane, this pressure is uniform on $S_p$. Moreover, either the boundary of surface $S_p$ is fixed or $S_p$ is split by a symmetry plane. Consequently, the contribution of all the line integrals in (12) to the stiffness matrix is always zero.

The usual arclength method [24] is used in order to proceed on the response curves. The nodal displacement vector $\{\mathbf{U}\}$ is split into two parts: one denoted $\{\bar{\mathbf{U}}\}$ contains the unknown degrees of freedom, the other denoted $\{\bar{\mathbf{U}}\}$ contains the prescribed degrees of freedom. The external force vector $\{\mathbf{\Phi}\}$ is split in the similar way: one part $\{\tilde{\mathbf{\Phi}}\}$ corresponding to $\{\bar{\mathbf{U}}\}$ contains the prescribed force components, the other part $\{\tilde{\mathbf{\Phi}}\}$ corresponding to $\{\bar{\mathbf{U}}\}$ contains the unknown reaction force components. Either the prescribed displacement or the external loading is assumed to be proportional

$$ \begin{aligned} \{\bar{\mathbf{U}}\} &= \{\bar{\mathbf{U}}\}^0 + \lambda\{\bar{\mathbf{U}}\}^{\text{ref}}, \\ \{\tilde{\mathbf{\Phi}}\} &= \lambda\{\tilde{\mathbf{\Phi}}\}^{\text{ref}}, \end{aligned} \quad (13) $$

where $\lambda$ is the control parameter, $\{\bar{\mathbf{U}}\}^{\text{ref}}$ and $\{\tilde{\mathbf{\Phi}}\}^{\text{ref}}$ denote reference prescribed quantities. Vector $\{\bar{\mathbf{U}}\}^0$ related to zero prescribed displacements does not change the value of $\{\bar{\mathbf{U}}\}$. The arclength method consists of moving forward on the response
curve by a given arclength $\Delta\ell$, see Wempner [WEM71] and Riks [RIK72], [RIK79]. The constraint equation is either of the following relations, depending on one has proportional prescribed displacement or loading:

$$ \begin{aligned} \|\Delta\tilde{\mathbf{U}}\|^2 + \Delta\lambda^2 \|\bar{\mathbf{U}}^{\text{ref}}\|^2 &= \Delta\ell^2, \\ \|\Delta\tilde{\mathbf{U}}\|^2 + \Delta\lambda^2 C_{\text{ref}}^2 &= \Delta\ell^2. \end{aligned} \quad (14) $$

In Relation (14b), the scalar $C_{\text{ref}}$ is a scale factor which makes the relation consistent dimensionally, it will be taken to be zero in the numerical examples. Combining relation (14) with the equilibrium equation leads to a quadratic equation in $\Delta\lambda$ [CRI91]. Eventually, the branch switching techniques have to be included in the numerical procedure in order to deal with limit and bifurcation points. The details of these techniques will be given for each examples treated in the next section.

### 3. Numerical examples

We consider some numerical examples to assess the capacity of the proposed approach to catch bubbles and wrinkles in membrane structures. In all examples the Newton-Raphson method is used for solving the matrix nonlinear equations of the problem. The path following is carried out either by displacement (14a) or force control (14b), use is also made of an extended version of the arclength method as proposed by Lam and Morley [25] to deal with the complex roots in the solution scheme.

#### 3.1. Inflated torus

In the first example, we consider the buckling of a toric membrane subjected to an internal pressure. This example is chosen in order to assess the ability of the numerical procedure to correctly compute both limit and bifurcation points.

In the reference configuration, which is assumed to be stress free, the middle surface of the toroidal membrane is generated by rotating the circle of radius $r_0 = 0.1$ m with its center at the distance $R_0 = 0.4$ m about the axis of symmetry, the z-axis (see Fig. 2). The membrane is of constant thickness $h_0=10^{-4}$ m. The numerical computation is carried out assuming a neo-Hookean material (8) with Young modulus $E = 4.0$ MPa and Poisson ration $\nu = 0.49$, where $E$ and $\nu$ are related to Lamé's constants $\lambda$ and $\mu$ by $\lambda = E\nu/(1 + \nu)/(1 - 2\nu)$, $\mu = E/2/(1 + \nu)$.

Taking into account the symmetry allows us to consider the fourth of the torus only. The mesh is made of 600 eight-node quadrilateral elements, generated by 50 elements along the large circumference (the circle of radius $R_0$) and 12 elements along the small one (the circle of radius $r_0$). Full integration is used throughout the paper: $3 \times 3$ Gaussian points in eight-node quadrilateral elements and 7 points in six-node triangular elements. Fig. 3 shows the displacement at point I defined in Fig. 2 versus the internal pressure.

The numerical computation shows that critical points appear on the fundamental branch, beyond a certain deformation level. The detection of such critical points is based on the singularity

<figure><img src="image_3.png" /><figcaption>Fig. 2. Reference geometry of the torus. Definition of point I.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Pressure versus displacement curves for the torus.</figcaption></figure>

of the tangent stiffness matrix $\mathbf{K}$, which is factorized as

$$ \mathbf{K} = \mathbf{L}\mathbf{D}\mathbf{L}^{\mathrm{T}}, \quad (15) $$

where $\mathbf{L}$ is a lower triangular matrix with unit diagonal elements and $\mathbf{D}$ is a diagonal matrix. Since the number of negative eigenvalues of $\mathbf{K}$ is equal to the number of negative diagonal elements (pivots) of $\mathbf{D}$, the critical points are determined by counting the negative pivot number.

Each critical point has to be isolated in order to determine its nature: limit point or bifurcation point. To do this, the current arclength $\Delta\ell$ is re-estimated several times using a dichotomy-like method, such as those proposed in [26].

A suitable way to distinguish a limit point from a bifurcation point is to calculate the current stiffness parameter
<figure><img src="image_5.png" /><figcaption>Fig. 4. Fundamental solution and buckling modes of the torus.</figcaption></figure>

defined as [27]

$$k = \frac{\tilde{\Phi}^{\text{ref}} \cdot \tilde{U}_{\phi \text{ref}}}{\tilde{U}_{\phi \text{ref}} \cdot \tilde{U}_{\phi \text{ref}}}, \qquad (16)$$

where $\tilde{U}_{\phi \text{ref}} = \tilde{\tilde{K}}^{-1} \tilde{\Phi}^{\text{ref}}$, $\tilde{\tilde{K}}$ is the square matrix extracted from the tangent stiffness matrix $K$ according to the partition $U = (\tilde{U}, \tilde{U})$ described in Section 2. The sign of parameter $k$ changes when passing a limit point, whereas it remains unchanged when passing a bifurcation point. On the fundamental branch of the considered example, it is found that the pressure reaches one limit point at $p = 1030$ MPa and all the bifurcation points occur after this point.

The switching on a bifurcated branch is performed by using the mode injection [28]: at the first step of a bifurcating branch, the eigenvector $\tilde{Z}$, solution of $\tilde{\tilde{K}} \cdot \tilde{Z} = 0$, is computed and the following predictions are used:

$$\delta\lambda = 0, \quad \delta\tilde{U} = \pm \Delta\ell \tilde{Z} / \|\tilde{Z}\|. \qquad (17)$$

Fig. 4 displays the four bifurcation modes corresponding to the appearance of asymmetrical bubbles on the torus.

It should be noted that no bifurcation is obtained if one considers a Saint-Venant Kirchhoff material defined by the strain energy (7).

## 3.2. Pinched hemisphere

Now consider the problem of a hemispherical membrane fixed along its rim and subjected to a pressure over one face and a vertical concentrated force at the center. An analytical study of axisymmetrical solutions was performed by Szyszkowski and Glockner [29], who obtained the deformed configuration by relaxing the constraint of inextensibility of the latitude circles so as to allow them to grow shorter in regions where the hoop stress vanishes. By definition, the axisymmetrical approach does not

exhibit wrinkles. Experimental investigations of this problem were also made by Szyszkowski and Glockner in [30], who showed that wrinkles are displayed for various pinching forces. A recent numerical computation was carried out by Stanuszek [21] who developed a new method for catching wrinkles using the cable analogy.

In this paper, we aim to compute the wrinkles using a more traditional bifurcation analysis. For the purpose of comparing with the results in [21,30], the present study uses the same geometry and material data as therein. Thus, in the reference configuration, the thickness of the membrane is taken as $h_0 = 25.4$ µm and the radius as $R_0 = 0.155$ m. Note that the radius-thickness ratio is $R_0/h_0 = 6102$, which means that the membrane is very thin. The membrane is made of a neo-Hookean material defined by (8), with Young's modulus $E = 2.7$ MPa and Poisson's ratio $\nu = 0.4$.

Due to the symmetry, only one quarter of the hemisphere is modelled with a mesh of 25 × 25 elements (see Fig. 5). All of them are eight-node quadrilateral elements, except for those at the center, which are six-node triangular ones.

The loading is applied in two stages: first, the membrane is inflated to a given pressure $p$, and then a concentrated force $F$ is applied at its vertex. Fig. 6 shows the axisymmetrical deformed shape corresponding to the fundamental solution, where there is no bifurcation. It is seen that the sinking at the apex does not show any fold.

When computing the bifurcated branches, one encounters severe computational difficulties due to complex roots which occur repeatedly when the quadratic equation in the arclength method is solved. It is found that an efficient way to cope with these complex roots is to modify the solution scheme according to Lam and Morley [26]. The main idea of the procedure is to project the residual force onto the external load vector. At a current iteration where complex roots occur, the residual force is split into one component in the load direction and another component orthogonal to this load.
<figure><img src="image_6.png" /><figcaption>Fig. 5. Mesh of the pinched hemispherical membrane.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Fundamental solution for a hemispherical membrane loaded by vertical force.</figcaption></figure>

The last component is mainly responsible for the complex roots and which should be eliminated. The algorithm is modified as summarized below.

To make sure that the quadratic equation in $x$ gives real roots, factor $\eta$ is chosen at 5% of $|\eta_2 - \eta_1|$, where $\eta_1$ and $\eta_2$ are the roots of another quadratic equation in $\eta$.

Fig. 7 shows the bifurcated solutions of the pinched hemispherical problem for two inflation pressures, $p=120$ and $5000$ Pa (the same values as in Refs. [30,21], respectively).

In both cases, the folds are regularly distributed in the vicinity of the concentrated load. In the low pressure case ($p = 120$ Pa), the membrane cannot bear a large pinching force, so that the folds are very fine and the depression is quite shallow. The computed deformed shape shown in Fig. 7a agrees very well with the experimental tests presented in [30]. In the high pressure case ($p = 5000$ Pa), the pinching force $F$ can be given larger values and the folds are more visible, as shown in Fig. 7b.

The dimensionless load $F/p\pi R_0^2$ is plotted versus the dimensionless deflection $w/R_0$ in Fig. 8. As expected, the bifurcated curve is very close to the fundamental curve since the deflection $w$ does not change very much.

### 3.3. Square thin film membrane subjected to in-plane shear loading

This section deals with a square membrane clamped along the bottom edge and subjected to a horizontal shear displacement on the top edge. The sheared membrane is known as one of the most difficult tests. Numbers of works have been devoted to this problem in the literature. Rossi et al. [31] used a membrane element with a ‘no-compression’ material based on a modification of the standard linear material and introduced a “dynamic relaxation” into the process. Löhnert et al. [32] also used membrane elements and split the strain tensor into a purely non-wrinkling part and a wrinkling strain to study shearing effect in orthotropic membranes. Wong and Pellegrino [15,16], studied a sheared membrane by using shell elements. Their numerical simulation is done in three main parts, namely (i) to apply a small uniform prestress to the membrane by moving the upper edge in the $y$-dimension; (ii) to predict buckling modes

Iteration loop

[...]

In the standard arclength method, solve the quadratic equation.

If the roots are complex, then

(i) Recast the displacement $\delta U$ as $\delta U = x\delta U_H + \eta\delta U_\Phi$, where $\delta U_H = -K^{-1}H$, $H$

is the component of the residual force orthogonal to reference load $\{\tilde{\Phi}\}^{\text{ref}}$

(see Eq. (13b)), $\delta U_\Phi = -K^{-1}\Phi$, $x$ and $\eta$ are scalar factors.

(ii) The new incremental displacement $\delta U$ is replaced into the constraint equation of the standard arclength method, which leads to a new quadratic equation to be solved for $x$. One chooses the root that enables the solution point to advance in the desired direction.

End if

[...]

End of iteration loop.
<figure><img src="image_8.png" /><figcaption>Fig. 7. Bifurcated solutions for a hemispherical membrane loaded by vertical force.</figcaption></figure>

by using a combination of eigenvectors to generate geometric imperfections; and (iii) to continue the post-buckling analysis. Tessler et al. [19] also investigated sheared membrane. They used an updated Lagrangian shell formulation and added to the structure an imperfection which varies as a function of the membrane thickness.

In our computations, one has to slightly stretch the membrane before shearing in order to give it a preliminary stiffness at the beginning of the process. Moreover, it should be noted that, contrary to the previous examples, here it is essential to apply the arc-length method with the displacement control (13a), not with the force control (13b); the arc-length parameter is then the shear displacement. As a matter of fact, it turns out that there are convergence difficulties when the arc-length parameter is related to the shearing force rather than the shear displacement. Use is also made of the treatment of complex roots described in Section 3.2.

In its reference configuration, the membrane side is 0.25 m, the thickness is $h_0=25$ µm. The material is assumed to be of the neo-Hookean type, with Young's modulus $E = 2500$ MPa and Poisson's ratio $v=0.34$. The computed mesh contains 40 × 40 eight-node quadrilateral elements. First, a stretching of 0.5% of the length is prescribed in the y-direction, and in the second stage the shear displacement is applied in the x-direction up to 5% of the side length. Whereas the fundamental solution which corresponds to an in-plane deformed shape is rather trivial, the bifurcated solution displays out-of-plane wrinkles which are much more difficult to be obtained. Fig. 9 gives the deformed shape of the bifurcated solution, where four large folds and two hardly noticeable smaller ones are observed in one diagonal direction.

Fig. 10 depicts the shear displacement versus the out of plane deflection at point A of coordinates (0, 12.5, 0). As the membrane was stretched before being sheared, the wrinkles appear only beyond a certain shear displacement level.

### 3.4. Square airbag

Eventually, let us consider the problem of a pressurized airbag. Consider an initially flat square airbag inflated by an

<figure><img src="image_9.png" /><figcaption>Fig. 8. Load-displacement curves for pinched hemispherical membrane.</figcaption></figure>

internal pressure. Although the airbag inflation is essentially a dynamic problem, here we confine ourselves in a static analysis framework.

The specific difficulty of this problem comes from the combination of two facts: (i) the whole boundary of the airbag is completely free to move in its plane, and (ii) the pressure is acting normally to the membrane surface. As a consequence, the stiffness matrix presents a high singularity at the very first computational step, when the pressure is still very low. To overcome this numerical difficulty, it is necessary to stretch the structure along x and y directions by means of dead forces applied on the outer boundary. Afterwards, those forces are gradually removed when keeping the internal pressure at a fixed value.

The reference thickness is $h_0 = 10^{-4}$ m and the reference diagonal length is 1.2 m (Fig. 11), the material is of neo-Hookean type (8) with Young's modulus $E = 588$ MPa and Poisson's ratio $v = 0.4$ (all the numerical values are those used in [5]). The airbag is subjected to an internal pressure $p = 5000$ Pa. By taking into account symmetry considerations, one can reduce to studying one upper quarter of the airbag only. The deformed
<figure><img src="image_10.png" /><figcaption>Fig. 9. Sheared membrane.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. The shear displacement versus the out of plane deflection.</figcaption></figure>

shapes obtained with coarse meshes—4 × 4 and 5 × 5 eight-node quadrilateral elements—as shown in Fig. 12 are in good agreement with Contri’s results [5].

Fig. 13 shows a more realistic solution of the airbag using 25 × 25 eight-node quadrilateral elements. The deformed shape is close to that computed by Stanuszek [21]. Furthermore, beside two deep central folds near point B, other small wrinkles are found on the sides, this result is in very good agreement with experimental results.

Table 1 compares the following displacements obtained with
several different meshes: (i) the deflection $w_M$ at the center

point M of the airbag, (ii) the displacement $u_A$ in the x-axis at
the corner A, and (iii) the displacement $u_B$ in the x-axis at the
mid-point B of the airbag (see Figs. 12 and 13). The number of
elements in Table 1 is related to one upper quarter of the airbag.

Table 1 shows that the results from the present work agree quite well with the others. It should be noted that in this example, the wrinkles appear naturally while the dead forces are gradually removed. No bifurcation is detected during the computation, so that here one need not have recourse to the bifurcation analysis.

**4. Conclusions**

The objective of this paper has been to propose a robust tool
to study the wrinkling in membrane structures. To this end, a
membrane element has been presented and the pure bifurcation
analysis has been carried out without introducing any imper-
fections in the structure. Also, a simple yet efficient technique
as proposed by Lam and Morley [24] has been incorporated in
the solution procedure to deal with the possible complex roots
when solving the arclength method, making thus possible the
treatment of the wrinkles by the bifurcation. The whole pro-
posed formulation has been programmed in the FORTRAN 90
language by the authors.

The selected numerical examples have shown the ability of
the proposed formulation to correctly predict the critical values
for the wrinkles to appear as well as the wrinkled regions. The
usual singularity of the stiffness matrix at the beginning of the

<figure><img src="image_12.png" /><figcaption>Fig. 11. Initial geometry of the airbag.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 12. Deformed airbag using: (a) 4 × 4 elements; (b) 5 × 5 elements in one fourth of the upper surface.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Deformed airbag using 25 × 25 elements.</figcaption></figure>

<table><caption>Table 1
Displacements at some particular points of the airbag and comparison with the results in the literature</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Bauer (quoted in [5])</th><th colspan="2">Contri and Schrefler [5]</th><th colspan="2">Kang and Im [33]</th><th colspan="3">Present work</th></tr><tr><th>16 Elements</th><th>25 Elements</th><th>16 Elements</th><th>25 Elements</th><th>16 Elements</th><th>25 Elements</th><th>625 Elements</th></tr></thead><tbody><tr><th>w<sub>M</sub> (m)</th><td>0.205</td><td>0.209</td><td>0.217</td><td>0.215</td><td>0.216</td><td>0.2145</td><td>0.2144</td><td>0.2245</td></tr><tr><th>u<sub>A</sub> (m)</th><td>0.033</td><td>0.04</td><td>0.045</td><td>0.043</td><td>0.042</td><td>0.0282</td><td>0.0265</td><td>0.0307</td></tr><tr><th>u<sub>B</sub> (m)</th><td>0.13</td><td>0.102</td><td>0.11</td><td>0.117</td><td>0.117</td><td>0.126</td><td>0.1207</td><td>0.1158</td></tr></tbody></table>

process has been avoided by preloading the structure, either with an artificial internal prestress or a real load or displacement prescribed on the boundary.

Although the wrinkle patterns have been correctly described, further analyses should be carried out in order to determine more quantitative results such as the amplitude, the wavelength and the number of wrinkles.

## References

[1] H. Wagner, Flat sheet metal girder with very thin metal web, Z. Flugtechnik Motorlurfetschifffahrt 20 (1929) 200–314.

[2] M. Stein, J.M. Hedgepeth, Analysis of partly wrinkled membranes, NASA TN D-813 (1961) 1–23.

[3] R.K. Miller, J.M. Hedgepeth, V.I. Weingarten, P. Das, S. Kahyai, Finite element analysis of partly wrinkled membranes, Comput. Struct. 20 (1985) 631–639.

[4] M. Hideki, O. Kiyoshi, A study on modelling and structural behaviour of membrane structures, in: Proceedings of the IASS Symposium on Shells, Membranes and Space Frames, Osaka 2 (1986) 161–167.

[5] P. Contri, B.A. Schrefler, A geometrically nonlinear finite element analysis of wrinkled membrane surface by no-compression material model, Commun. Appl. Numer. M. 4 (1988) 5–15.

[6] M. Fujikake, O. Kojima, S. Fukushima, Analysis of fabric tension structures, Comput. Struct. 32 (1989) 537–547.

[7] B. Tabarrok, Z. Qin, Nonlinear analysis of tension structures, Comput. Struct. 45 (1992) 973–984.

[8] W.W. Schur, Development of a practical tension field material model for thin films, 32nd Aerospace Sciences Meeting, Reno, NV, AIAA-94-0636, 1994.
[9] S. Kang, S. Im, Finite element analysis of wrinkling membranes, J. Appl. Mech. 64 (1997) 263-269.

[10] A.L. Adler, Finite element approaches for static and dynamic analysis of partially wrinkled membrane structures, Ph.D. Thesis, Departement of Aerospace Engineering, University of Colorado, Boulder, CO, 2000.

[11] A.L. Adler, M.M. Mikulas, J.M. Hedgepeth, Static and dynamic analysis of partially wrinkled membrane structures, 41st AIAA Structures, Structural Dynamics and Materials Conference, Atlanta, GA, AIAA-2000-1810, 2000.

[12] X. Li, D.J. Steigmann, Finite deformation of a pressurized toroidal membrane, Int. J. Non-linear Mech. 30 (1995) 583-585.

[13] X. Li, D.J. Steigmann, Point loads on a hemispherical elastic membrane, Int. J. Non-linear Mech. 30 (1995) 569-581.

[14] R. Barsotti, S.S. Ligaro, G.F. Royer-Carfagni, The web bridge, Int. J. Solids Struct. 38 (2001) 8831-8850.

[15] Y.W. Wong, S. Pellegrino, Amplitude of wrinkles in thin membranes, in: New Approaches to Structural Mechanics Shells and Biological Structures, Kluwer Academic Publishers, Dordrecht, The Netherlands, 2002, pp. 257-270.

[16] Y.W. Wong, S. Pellegrino, Computation of wrinkle amplitudes in thin membranes, 43rd AIAA Structures, Structural Dynamics and Materials Conference, Denver, Co, AIAA-2002-1369, 2002.

[17] Y.W. Wong, Analysis of wrinkling patterns in prestressed membrane structures, M.Phil. Thesis, Department of Engineering, University of Cambridge, Cambridge, UK, 2000.

[18] X. Su, F. Abdi, B. Taleghani, J.R. Blandino, Wrinkling analysis of a Kapton square membrane under tensile loading, 44th AIAA/ASME/ASCE/AHS Structures, Structural Dynamics and Materials Conference, Northfolk, Virginia, AIAA 2003-1985, 2003.

[19] A. Tessler, D.W. Sleight, J.T. Wang, Nonlinear shell modelling of thin membranes with emphasis on structural wrinkling, 44th AIAA/ASME/ASCE/AHS Structures, Structural Dynamics and Materials Conference, Northfolk, Virginia, AIAA 2003-1931, 2003.

[20] M. Stanuszek, J. Orkisz, Finite element analysis of membrane cable reinforced structures at large deformations, Adv. Non-linear Finite Element Methods 4 (1994) 67-76.

[21] M. Stanuszek, Finite element analysis of large deformations of membranes with wrinkling, Finite Elements Anal. Design 39 (2003) 599-618.

[22] T. Miyamura, Wrinkling on stretched circular membrane under in-plane torsion: birfurcation analyses and experiments, Eng. Struct. 23 (2000) 1407-1425.

[23] T. Belytschko, W.K. Liu, B. Moran, Nonlinear Finite Elements for Continua and Structures, Wiley, England, 2000.

[24] M.A. Crisfield, Non-linear Finite Element Analysis of Solids and Structures, vol. 1, Wiley, England, 1991.

[25] W.F. Lam, C.T. Morley, Arc-length method for passing limit points in structural calculation, J. Struct. Eng. 118 (1992) 19-185.

[26] J. Shi, G.F. Moita, The post-critical analysis of axisymmetric hyperelastic membranes by the finite element method, Comput. Method. Appl. Mech. Eng. 135 (1996) 265-281.

[27] P.G. Bergan, G. Horrigmoe, B. Krakeland, B. Soreide, Solution techniques for non-linear finite element problem, Int. J. Numer. Methods Eng. 12 (1978) 1677-1696.

[28] W. Wagner, P. Wriggers, A simple method for the calculation of post-critical branches, Eng. Comput. 5 (1988) 103-109.

[29] W. Szyszkowski, G. Glockner, Finite deformation and stability behaviour of spherical inflatables under axisymmetric concentrated loads, Int. J. Non-linear Mech. 19 (1984) 489-496.

[30] W. Szyszkowski, G. Glockner, Spherical membranes subjected to vertical concentrated loads: an experimental study, Int. J. Eng. Struct. 9 (1987) 183-192.

[31] R. Rossi, M. Lazzari, R. Vitaliani, E. Onate, Convergence of the modified material model for wrinkling simulation of light-weight membranes structures, in: E. Onate, B. Kröplin (Eds.), Textile Composites and Inflatable Structures, CIMNE, Barcelona, 2003, pp. 148-153.

[32] S. Löhnert, K. Tegeler, P. Wriggers, T. Raible, A simple wrinkling algorithm for orthotropic membranes at finite deformations, in: E. Onate, B. Kröplin (Eds.), Textile Composites and Inflatable Structures, CIMNE, Barcelona, 2003, pp. 119-122.

[33] S. Kang, S. Im, Finite element analysis of dynamic response of wrinkling membranes, Comput. Method. Appl. Mech. Eng. 173 (1999) 227-240.