Buckling and wrinkling of prestressed membranes

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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.


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][4][5][6][7][8][9][10][11]. The tension field theory can also be addressed from an energetical 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 postbucling analysis. Wong and Pelligrino [15][16][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.

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 . The reference and current positions of a material point on the membrane surface-X = X i e i and x = x i e i , respectively-are parametrized by the reference coordinates ( 1 , 2 ) as shown in the following chain: ( 1 , 2 ) ∈ e → X ∈ e ⊂ S 0 → x ∈ S, see Fig. 1. The mid-surface displacement field U is defined by U = x − X.
The covariant base vectors on the undeformed membrane are defined by The homologous vectors on the deformed membrane are defined by The deformation gradient tensor F relates the base vectors of the reference configuration to the base vectors of the deformed configuration as where 3 represents the through-thickness stretch. The strains are measured by the Green tensor E and its covariant components where Greek indices range from 1 to 2, and g and G are the metric tensor components in the deformed and undeformed configurations defined by Assuming a compressible hyperelastic behaviour of the membrane, the second Piola-Kirchhoff (symmetric) stress tensor is given by differentiating the strain energy function per unit volume w with respect to the Green strain E where Latin indices take the values 1, 2 and 3. The simplest constitutive model is described by the Saint-Venant Kirchhoff strain energy where and are the Lamé's constants and 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] where C = 2E + I is the right stretch tensor, J = √ det C, and and 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 and values.
Since the thickness of the membrane is small compared with other dimensions, one assumes the plane stress condition 13 = 23 = 33 = 0. With potentials (7) and (8), the conditions 13 = 23 = 0 are equivalent to E 13 = E 23 = 0, whereas the condition 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 as functions of in-plane covariant strains E . In the case of potential (8), the 2D stress-strain relationship is 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 U * , one has where = F is the first Piola-Kirchhoff (nonsymmetrical) stress tensor, f 0 the dead load per unit reference surface, p a possible pressure over the membrane and n the unit outward normal in the deformed configuration. The resulting discretized equations form a nonlinear matrix equation system with unknown {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 In the above, the integers a, b ∈ [1, NNE] and p, q ∈ [1,3] are determined from indices i, j ∈ [1, 3NNE] by i = 3(a − 1) + p, j = 3(b − 1) + q, and there is implicit summation over Greek indices , , , = 1, 2. The notation N a, means partial differentiation of the shape function N a with respect to . If the membrane is subjected to a pressure p which is a follower force, the following stiffness matrix is added to the previous one: ∀i, j ∈ [1, 3NNE], The integers a, b ∈ [1, NNE] and p, q ∈ [1, 3] are related to indices i, j ∈ [1, 3NNE] in the same way as in (11). The symbol x p∧q means x p∧q = pqm x m (the integers p and q being fixed, there is only one value for m so that pqm is not zero), is the unit outward normal to the reference element e , at a point on its boundary je . The surface integral over e is symmetric, whereas the line integral along je is skew-symmetric. Let the pressure be applied over a surface portion S p ⊂ 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 {U} is split into two parts: one denoted {Ũ} contains the unknown degrees of freedom, the other denoted {U} contains the prescribed degrees of freedom. The external force vector { } is split in the similar way: one part {˜ } corresponding to {Ũ} contains the prescribed force components, the other part { } corresponding to {U} contains the unknown reaction force components. Either the prescribed displacement or the external loading is assumed to be proportional In Relation (14b), the scalar C 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 [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.

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.

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)  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 × 3 Gaussian points in eightnode 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  of the tangent stiffness matrix K, which is factorized as where L is a lower triangular matrix with unit diagonal elements and D is a diagonal matrix. Since the number of negative eigenvalues of K is equal to the number of negative diagonal elements (pivots) of 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 is re-estimated several times using a dichotomylike 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 defined as [27] k whereŨ ref =K −1˜ ref ,K is the square matrix extracted from the tangent stiffness matrix K according to the partition U = (Ũ, 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 eigenvectorZ, solution ofK.Z = 0, is computed and the following predictions are used: 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).

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 = 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 a 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.  To make sure that the quadratic equation in x gives real roots, factor is chosen at 5% of | 2 − 1 , where 1 and 2 are the roots of another quadratic equation in . Fig. 7 shows the bifurcated solutions of the pinched hemisphere 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 R 2 0 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.

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 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 he 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 ration = 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.

Square airbag
Eventually, let us consider the problem of a pressurized airbag. Consider an initially flat square airbag inflated by an 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 = 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  shapes obtained with coarse meshes-4 × 4 and 5 × 5 eightnode 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  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.

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 imperfections 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 proposed 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   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.