Membrane wrinkling revisited from a multi-scale point of view

Membrane modeling in the presence of wrinkling is revisited from a multi-scale point of view. In the engineering literature, wrinkling is generally accounted at a macroscopic level by nonlinear constitutive laws without compressive stiffness, but these models ignore the properties of wrinkles, such as their wavelength, their size and spatial distribution. A new multi-scale approach is discussed that belongs to the family of Ginzburg- Landau bifurcation equations. By using the method of Fourier series with variable coefficients, several nonlinear macroscopic models are derived that couple the membrane response with equations governing the evolution of the wrinkles. Contrary to previous approaches, these macroscopic models are completely deduced from the “microscopic” shell model without any phenomenological assumptions. Some analytical and numerical solutions are discussed that prove the relevance of the presented modeling. A new class of models has been established. It permits to predict the characteristics of the wrinkles and their influence on membrane behavior.


Background
This paper deals with macroscopic modeling of very thin shells, also called membranes in the scientific literature and in everyday life. Membranes have increasing application fields such as spacecraft structures (antenna, telescope lenses, Gossamer structures…), civil structures, life jackets, electronics, biological tissues, textile composites [1][2][3][4][5][6]. The appearance of wrinkles is a major feature in the mechanical behavior of membrane, due to their vanishing bending stiffness.
Membrane wrinkling is a multi-scale phenomenon, but it is generally described by one-scale models that can be "microscopic" or "macroscopic". Microscopic models are simply based on elastic shell theory; see for instance [7][8][9][10][11]. Nowadays many commercial finite element codes permit to carry out such nonlinear shell computations. Shell analyses are able to describe the details of the membrane response: size, wavelength, orientation of the wrinkles, instability threshold, etc., but this leads to large scale numerical models that are moreover very difficult to be controlled in cases with a large number of wrinkles and therefore with a large number of equilibrium solutions. Indeed wrinkling can be seen as a small wavelength buckling phenomenon and, in such a case, there are several neighbor bifurcation modes and it is well known that this leads to compound response curves involving secondary bifurcations, see for instance [12,13] among an extensive literature. Even if the complexity of wrinkling patterns requires numerical treatments, one can mention several papers by Coman et al. [14][15][16] that were able to find analytically wrinkling modes in cases with inhomogeneous pre-bifurcation stresses.
The second group of numerical descriptions are pure membrane models that account for the effect of wrinkling on membrane behaviour. The bending stiffness is neglected and the compressive stresses are dropped: this yields nonlinear constitutive laws that are of unilateral type, the tensile behaviour being about linear and the compressive states being forbidden. This has led to many numerical studies; see for instance [17][18][19][20][21][22]. In most of these papers, the constitutive law relates simply membrane stress and strain and it distinguishes three states: slack, wrinkled and taut. Some authors prefer the method of Roddeman [23,24] that splits the deformation gradient into consistent membrane part and wrinkling part [25][26][27][28]. The partial differential equations deduced from these membrane models are not elliptic (or hyperbolic in the dynamical case) in the presence of a non-positive principal stress so that this problem is mathematically ill-posed. A well-posed problem can be obtained if the macroscopic model includes an internal length, for instance within Cosserat theory [29,30]. However this regularisation may be not necessary for an explicit dynamic computation. All these membrane models have a phenomenological character, the only mention to full shell models being the vanishing of the membrane compressive stress. A variant has been introduced by [31] and used for buckling problems of sheets under residual stresses generated by manufacturing processes [32].
The aim of this paper is to revisit membrane modelling from a multi-scale point of view. It will be clearly established that a consistent wrinkling model depends both on microscopic quantities (wavelength, orientation of the wrinkles…) and macroscopic data like boundary conditions, size and shape of the structure… In this respect, a new multi-scale membrane model was briefly presented [33] that couple a nonlinear 2D membrane model with an amplitude equation governing the evolution of wrinkles. Amplitude equations to describe spatial pattern formation follow generally from the asymptotic bifurcation analyses of Ginzburg-Landau type [34][35][36][37]. Here a slightly different method will be applied, where the nearly periodic fields are represented by Fourier series with slowly varying coefficients [38][39][40][41]. In other terms, we use a multi-scale modelling method whose result is a generalised continuum including an internal length and where the macroscopic stresses are Fourier coefficients of the microscopic stress.
In this paper, the new membrane model of [33] is discussed in detail. Some additional variants, several analytical and numerical solutions and the corresponding numerical schemes will be presented. The new membrane models including wrinkling will be deduced from Föppl-Von Karman plate theory. These analyses will be carried out by applying standard techniques of multi-scale bifurcation analyses: asymptotic multi-scale method for the linear bifurcation analysis and multi-scale Fourier approach leading to nonlinear models valid away from the wrinkling threshold. Note that, contrarily to most of macroscopic membrane models, the present models are deduced from the shell equations without any phenomenological assumptions. The bending stiffness effects are included not only to define the wrinkling wavelength, but also to predict the macroscopic evolution of the buckling pattern.
The paper is organized as follows. In section "About the origin of short wavelength instabilities", the origin of the wrinkling instability is briefly recalled, in order to underline its multi-scale origin and the generic character of this instability mechanism. section "Multi-scale nonlinear models for membrane wrinkling" is devoted to deducing new nonlinear wrinkling models from Föppl-Von Karman plate theory. Then two analytical solutions of these models will be presented (section "Two analytical solutions for clamped rectangular membranes"), the first one for defining the bifurcation features (bifurcation stress, wavelength and shape of the envelope) and the second one for a post-bifurcation analysis characterizing the evolution of the size of the wrinkles as a function of the applied load. The numerical analysis is presented in section "Numerical implementation" (techniques of implementation) and section Two numerical solutions where two examples of plate wrinkling will be solved. Finally, in section "Linear wrinkling analysis revisited by a double scale asymptotic analysis", a more classical asymptotic multi-scale analysis will be performed in a special case and the results will be compared with those arising from the multi-scale Fourier approach.

About the origin of short wavelength instabilities
Spatially periodic instability patterns are encountered in a lot of physical problems. The instability often occurs in domains whose dimension is smaller than the others, typically in rectangles or parallelepipeds with a small aspect ratio. The instability wavelength can be of the same order as the smallest dimension of the domain as in Rayleigh-Bénard convection [42,43].
The instability wavelength is not always directly related to a structural dimension. A well known typical case of periodic patterns is the flexible beam resting on a Winkler foundation that is also referred as Swift-Hohenberg equation [44][45][46][47][48][49][50]. The governing equation has the following form in 1D linear case: The wrinkling modes v(X) and the wrinkling load λ are then given by: The instability wavenumber Q or the instability wavelength ℓ = 2π/Q corresponds to the minimum of the neutral stability curve Q → λ(Q): A physical example of such instability is the symmetric microbuckling mode of long fiber composite materials that are represented by an infinite stack of hard and soft layers (fiber, matrix), see Figure 1 and Rosen [51]. In this problem, the coefficients of (1) are related to Young moduli and thickness's of the layers (A≅E F h 3 F , B ≅ E M /h M ), the microbuckling wavelength is related to the width of the fiber and of the matrix by the following formula: ℓ≅h Quite similar wrinkling phenomena with a small wavelength can be found in studies about thin films coupled with compliant substrate, see for instance [52][53][54][55], or in the buckling of sheets due to residual stresses and forming processes [56,57]. In this second class of instability problems, the wavelength is larger than the microscopic lengths h F , h M because of a large stiffness ratio. Here is the origin of the multi-scale behavior: the characteristic distance in the width is smaller than the longitudinal wavelength and moreover these longitudinal waves are generally modulated on larger distances. The wrinkling of membrane belongs to a third class of instability problems that can be represented by a partial differential equation in a 2D domain [3]: The instability modes are in the form V(Y)cos(QX+φ) and their instability wavelength 2π/Q is generally small with respect to the characteristic length in the transverse direction. The Equations (1) and (4) can be related by assuming that the coefficient B is an eigenvalue of the differential operator − C∂ 2 /∂Y 2 , or equivalently B ¼ C=L 2 Y , where L Y is a macroscopic characteristic length in the transverse direction.
In the problem of a plane membrane of thickness h submitted to a tensile stress σ Y in the transverse direction, A≅ Eh 3 , C = σ Y h and σ X = − λ/h is a compressive stress in the X-direction so that the wrinkling wavelength is ℓ 3]. Hence this wavelength is much smaller than the transverse length L Y because a membrane is very thin (h < <L Y ) and it decreases when the tensile stress increases. From Equation (3), one can deduce the critical compressive stress σ X j j≅ ffiffiffiffiffiffiffiffi ffi Eσ Y p h=L Y that is much smaller than the prescribed tensile stress σ y if this stress σ y is not too small. It is worth to mention that the wrinkling stress depends both on the microscopic length h and on the macroscopic transverse one L Y .
Note that the Equation (4) seems to be more or less generic. It has also been established in the antisymmetric microbuckling of long fiber composite [51,58], with A≅E F h F 3 , C≅E M h F (this antisymmetrical mode occurs when the thickness of the fiber and of the matrix are of the same order h M ≅h F ). Hence the instability wavelength is ℓ where the transverse length L Y can be associated with the composite plate thickness or to the ply thickness [58]. Note that the nonlinear versions of these fiber microbuckling models permit to predict the compression failure of long fiber composite materials [59][60][61]. Hence a wrinkling mode follows from two stiffness's: first the bending stiffness A, that is very small for a thin membrane or a fiber in a composite material, second a transverse stiffness C that can be purely elastic in the antisymmetrical microbuckling or the so called geometric stiffness due to the tensile stress in the wrinkling of a membrane. Clearly, all these instability problems involve several length scales. That is why we propose several multi-scale approaches to analyze the behavior of membranes in the presence of wrinkling.

Methods
A family of nonlinear membrane models is presented that is based on a multi-scale approach. Unlike macroscopic models of the literature, there is no phenomenological assumption and the models are completely deduced from the full plate model. The deduction method relies on the concept of Fourier series with slowly varying coefficients [38,40], whose principle is to work with envelopes of the spatial oscillations. These envelopes are solutions of systems of nonlinear partial differential equations established in the next sections "Multi-scale nonlinear models for membrane wrinkling". These equations can be solved, sometimes analytically and in general numerically, as sketched in the section "numerical implementation".

Multi-scale nonlinear models for membrane wrinkling
The method of Fourier series with slowly variable coefficients A new multiple scale approach: the method of Fourier series with slowly variable coefficients We present here the methodology that will be used to deduce macroscopic nonlinear models of membrane wrinkling. For simplicity, this first discussion is limited to onedimensional case (space variable X ∈ IR) and to a one-dimensional beam model of Von Karman type that was studied in numerous papers [38][39][40][41][44][45][46][47][48][49]: Let us suppose that the instability wavenumber Q is known. Within this method, the unknown field U X Þ , whose components are axial displacement, transverse displacement, resultant stress, curvature and membrane strains, is written in the following form: The new macroscopic unknown fields U m (X) vary slowly on a single period X; X þ 2π Q h i of the pattern oscillation. As pictured in Figure 2, at least two functions U 0 (X) and U 1 (X) are necessary to describe nearly periodic patterns: U 0 (X) can be identified to the mean value and U 1 (X) represents the envelope or the amplitude of the spatial oscillations. The first envelope U 0 (X) is real-valued while the next ones are complex. In the whole paper, we limit ourselves to two envelopes U 0 (X) and U 1 (X). Most of the time, we shall assume that they are real.
The derivation operators are calculated exactly according to the following rule: It was established [38,40] that it is necessary to keep some spatial derivatives of the envelopes, despite of the assumption of slowly varying envelopes d/dX < <Q. This will be re-discussed in the membrane case. At first all the derivatives are kept as in (7) (8), but some ones can be dropped later.

Application to the membrane constitutive law
The application of the Fourier method to simple nonlinear equation is straightforward and it can follow from a simple identification of Fourier coefficients. For instance let us consider the 1-D membrane constitutive law (5-b) that relates the membrane stress n, the membrane strain γ and the displacement (u,v). From (5-b), one can deduce a macroscopic constitutive law for the m-th Fourier envelope: Especially the constitutive law for the mean stress n 0 (X) will be very useful. In the case of two real envelopes (u 0 , v 0 , n 0 , γ 0 , K 0 ) and (u 1 , v 1 , n 1 , γ 1 , K 1 ), it can be expressed as: Let us mention the last two terms are positive, what means that the membrane stretches out when it wrinkles.

Energetic approach
There is another manner to derive the full macroscopic model [38,40] by starting from the potential energy. In the elastic beam problem (5), the initial potential is given by: As a consequence of the assumptions of slowly varying envelopes, these envelopes are assumed to be constant on each period so that only terms corresponding to the harmonic zero of K 2 ,γ 2 ,v 2 and v 4 contribute to the approximated values of the potential energy. According to Parseval formula, the potential energy of the macroscopic model can be written as (the contribution of v 4 is omitted here for simplicity, full reduced models can be found in [40]): Further simplifications can be introduced in this model if one looks only at local bending instabilities from one-dimensional elastic states. In this respect, we can drop the mean deflection v 0 (X) and the envelope of the axial displacement u 1 (X). In this framework, the curvature and membrane strains of the wrinkled beam are approximated by (see (7) (9) and suppress the imaginary terms): A last approximation can be introduced by disregarding d 2 v 1 dX 2 in the potential energy. This latter approximation has been done in previous works [40] and it is justified by the slow spatial variations of the envelope. Finally, the potential energy can be approximated in the following form: Thus the extrema of the macroscopic energy (14) are solutions of the following system:

Comments
Hence the membrane law (15-b) accounts for the wrinkling oscillation v 1 (X) that is governed by a sort of bifurcation Equation (15-c) looking like a Ginzburg-Landau equation.
The full model is a nonlinear system coupling membrane behavior and evolution of wrinkles. It has been established [38,40] that the Fourier approach generalizes the asymptotic Ginzburg-Landau method, but it is not limited to the neighborhood of the bifurcation. This modeling by few Fourier envelopes can be applied with several levels of approximation. The system (15) is the simplest possible with an internal length. It can still be simplified by neglecting the derivatives of v 1 in (15-c): this leads to a nonlinear relation between membrane stress and strain, see [40] § 4.4. This pure membrane model will be extended in the 2D case in the following. It is also possible to include more harmonics, see [40] § 3.2, with obviously the cost of a greater complexity. The model with a complex envelope v 1 (X) has been evaluated in [62]: it permits for instance a better account of the phase in the bulk and improves the response near the boundary.
The previous approach can be considered as a multi-scale method or a computational homogenization technique. The account of the local behavior is very simplified by the assumption of a harmonic local variation and it is described only by the instability wavenumber Q. If the Equation (15-c) is reduced to its linear version without modulation of the envelope, ones recovers the classical approach (2) (3) that gives the wavenumber at the first bifurcation by minimizing the critical load. When solving a nonlinear macroscopic problem such as (15), the wavenumber Q has to be prescribed and this could be considered as a weak point of the present macroscopic approach. Nevertheless it is known that, in a cellular bifurcation problem, many solutions can exist [45,46,63], each one being characterized by its wavelength. Furthermore a model with a complex envelope permits to predict a wavelength slightly deviating from the one a priori prescribed [62]. Higher order harmonics could also be accounted: for instance, it was established in [40] that a rough account of the second harmonic is sometimes necessary to recover a consistent post-bifurcation behavior. But this will be not necessary in the case of membrane wrinkling.
In the next parts, the interaction membrane-wrinkling will be modelised in a bidimensional case within a framework similar to (14) (15), the starting model being the Föppl-Von Karman equations.

Föppl-Von Karman plate equations as a microscopic model
The main objective of the paper is to obtain nonlinear membrane models where the wrinkling is described by a bifurcation equation deduced from the initial plate model.
The method of Fourier series with slowly variable coefficients is applied to deduce the sought model from the well known Föppl-Von Karman equations for elastic isotropic plates that will be considered here as the reference model: where u = (u,v) ∈ IR 2 is the in-plane displacement, w is the deflection, N and γ are the membrane stress and strain. With the vectorial notations (N→ t (N X N Y N XY ), γ→ t (γ X γ Y 2γ XY )), the membrane elasticity tensor is represented by the matrix Eh . The corresponding energy E can be split into a membrane part E mem and a bending part E ben , as follows: Macroscopic modeling of membrane strain As explained previously, the unknown fields U(X,Y) are expressed in terms of two harmonics: the mean field U 0 (X,Y) and the first order harmonicsr, The second harmonic should be taken into account to recover the results of the asymptotic Ginzburg-Landau bifurcation approach, see [40]. Nevertheless the second harmonic does not contribute to the membrane energy in the present case, because the rapid onedimensional oscillations e iQX are inextensional so that N 2 = 0, w 2 = 0. Hence the second harmonic does not influence the simplest macroscopic models. For simplicity, the details of this calculation are omitted. A unique direction OX for the wrinkling oscillations is chosen in the whole domain. Of course this assumption is a bit restrictive and should be removed in the future. A true multi-scale approach should include two levels of modelisation as in the FE2 method [64], but with a basic cell that is not a priori known [65], what should require a rather intricate management. Thus a more realistic goal is to first discuss the multiscale approach with a given wrinkling wavelength and a given direction of the wrinkles.
The derivation rules (7) (8) can be extended easily in a bi-dimensional framework. For instance, the first Fourier coefficient of the gradient and the 0 th order coefficient of the strain (i.e. its mean value on a period) are given by: As in formula (10), the strain is split into a classical part γ FK that has the same form as the original Föppl-Von Karman model (16) and a wrinkling part γ wr which depends only on the envelope of the deflection w 1 .
Let us now turn to a simplified version of the displacement-strain law (19) (20), in the same spirit as for the 1D model (15). First the displacement field is reduced to a membrane mean displacement and to a bending wrinkling, i.e. u 1 = 0, w 0 = 0. This means that we account for the influence of the local buckling on the membrane behavior, but not for the coupling between local and global buckling as in [66,67]. Second the deflection envelope w 1 (X, Y) is assumed to be real, which disregards the phase modulation of the wrinkling pattern. Thus the envelope of the displacement has only three components u 0 = (u 0 , v 0 ) and w 1 that will be rewritten for simplicity as u; Þ . The simplified version of the strain field becomes: The simplified membrane strain formula (21) (22) is quite similar to that of the initial Von Karman model. It is split, first in a linear part ε(u) that is the symmetric part of the mean displacement gradient corresponding to the pure membrane linear strain, second in a nonlinear part γ wr (w) more or less equivalent to wrinkling stain of [23]. The main difference with the initial Föppl-Von Karman strain (16) is the extension Q 2 w 2 in the direction of the wrinkles. This wrinkling strain is always positive and this corresponds to a stretching. In the case of a compressive membrane strain, this wrinkling term leads to a decrease of the true strain.
As for the 1D model (14), we limit ourselves to the 0 th order harmonic to compute the reduced membrane energy:

Macroscopic bending energy
We saw two possible ways to get a reduced-order model via the technique of Fourier series with slowly variable coefficients: either identify the Fourier coefficient as in (9), or simplify the energy by keeping only the 0th order term in the energy, as in (12). Here we shall prefer the second approach which permits to provide a formulation easier to be managed for the numerical discretisation. The computation of the energy is based on the fact that only the 0th order harmonic φ 0 of a function φ has a non zero mean value: The identity (24) is applied to the two terms of the bending energy in the same framework as in (15) (21) (22) Because of the assumption of a real envelope w = w 1 , the first term of the bending energy is obtained: The second term of the bending energy φ B 0 is computed in the same way: To be consistent with the approach in the previous 1D case, the derivatives of order three or four in the differential equations are neglected. This has been justified in [39]: indeed spurious oscillations can appear in the response of the macroscopic model if these high order derivatives are kept. Finally the bending energy in the simplified framework is reduced as: Three macroscopic membrane wrinkling models In this section, one establishes the Partial Differential Equations of the macroscopic models associated with the strain energies previously presented. In fact, we have not defined a closed model, but a family of models that can depend on the number of harmonics and various other assumptions, such as the assumption of a real wrinkling envelope or the hypothesis w 0 = 0, which means no bending before wrinkling. Three cases will be considered. First the model presented in [33] couples a linear membrane behavior with a real envelope governed by a sort of Ginzburg-Landau equation. It corresponds to (15) in the 1D case and can be considered as a reference model because it is simple and able to account for the influence of wrinkling on membrane behavior. Next a pure membrane model will be derived by neglecting the dependence with respect to the derivatives of the envelopes. Last we shortly describe a model with a complex envelope, as the one studied in [62] in the 1D case that is a relatively simple improvement of the reference model.
A reference macroscopic membrane model The first proposed macroscopic membrane model follows from the minimum of the potential energy. In the absence of body forces, the sum of the membrane energy (23) and of the bending energy (27) is stationary at equilibrium. Let us recall that these energies associate 0th order harmonic for membrane quantities and a real first order harmonics for the deflection, as in the 1D model (15). In other words, it permits to couple a spatially modulated wrinkling with a linear membrane model: δE ben þ δE mem ¼ 0; for any virtual displacement that is zero at the boundary. This gives: After straightforward calculations, one obtains the partial differential equations of the macroscopic problems in the following form, the wrinkling membrane strain γ wr (w) being given in (22): The nonlinear model (30) (31) (32) couples nonlinear membrane equations with a bifurcation Equation (32) satisfied by the envelope of wrinkling pattern. It extends the previous analysis of the 1D case that couples a beam membrane with a one-dimensional Ginzburg-Landau equation. Hence the bifurcation Equation (32) is a sort of bi-dimensional Ginzburg-Landau equation, but it differs from the amplitude equation of Segel-Newell-Whitehead [42,43] who consider cases where the pre-bifurcation state is invariant under rotation. We shall see that finite element discretisation of (30) (31) (31) is straightforward. Two analytical solutions will be also presented in the section "Two analytical solutions for clamped rectangular membranes". Within the nonlinear model (30) (31) (32), one recovers two ideas of classical macroscopic membrane models. First the splitting between membrane and wrinkling strain of Roddeman theory [23] has been deduced from Föppl-Von Karman Equation, see (31), without any phenomenological assumption. Then the final bifurcation Equation (32) includes an internal length, which permits to retrieve the multi-scale instability analysis of Part 2. Finally this model can be qualitatively compared with the models of Banerjee et al. [29,30], where an internal length is introduced via Cosserat theory. Body forces and boundary forces can be introduced easily by the same procedure. This requires that these forces can also be put in the form of Fourier series with slowly variable coefficients. For instance if these forces vary slowly at the scale of the wrinkles, this leads to a classical body force in the membrane Equation (30), as well as in the corresponding boundary conditions. In the same way, transverse forces can also be accounted for within the more complete model entitled "A more sophisticated macroscopic membrane model with a complex envelope".

A pure membrane model
The reference membrane model (30) (31) (32) is not a pure membrane model because the Equation (32) includes spatial derivatives of the envelope of the wrinkles. It is natural to try to recover a pure membrane model, where the kinematic unknown is only the in-plane displacement, as suggested in [40], § 4.4. By dropping all the spatial derivatives in (32), one gets a bifurcation equation (N X + DQ 2 )w = 0 that can be transformed in a perturbed bifurcation equation as: From (33), one can obtain the deflection as a function of one component of the membrane stress. In (33), δ is a small perturbation parameter that transforms the perfect bifurcation equation into a perturbed bifurcation, what is more convenient for numerical path following calculations. If one simplifies the wrinkling strain (22) γ wr (w) = Q 2 w 2 e x ⊗ e x and if one combines (31) and (33), one can drop the deflection and deduce a nonlinear relation between membrane strain ε(u) and membrane stress N. With account of the balance of membrane forces, the obtained full model is restricted to: The model (34) is consistent with the pure membrane theories of the literature, see for instance [23], because N X cannot be lower than the wrinkling stress − DQ 2 .Generally, this wrinkling threshold is approximated by zero so that the membrane does not admit compressive stresses. In sections "Two analytical solutions for clamped rectangular membranes" and "Two numerical solutions", this model will be compared with our reference model (30) (31) (32).

A more sophisticated macroscopic membrane model with a complex envelope
More sophisticated models can be introduced in a rather easy way. For instance models with five envelopes (harmonics 0, ±1, ±2) have been presented in [39,40], starting models being respectively 2D hyperelasticity and the beam model of section The method of Fourier series with slowly variable coefficients. The reference model of section A reference macroscopic membrane model can be improved at least in two ways. First one can reintroduce the mean deflection w 0 to account for a coupling between local wrinkling and global buckling, as in [66]: in this manner, one will associate an envelope equation with the full Föppl-Von Karman plate equations. Moreover, this envelope can be a complex one for a better account of the phase field and of the boundary behavior [62]. To keep a relatively simple model, a more questionable hypothesis will be done: the membrane behavior will be described only by the 0 th order harmonic, i.e. u 1 = 0, N 1 = 0, γ 1 = 0. Hence the remaining unknown fields will be: i.e. the same variables as in the basic Föppl-Von Karman equations completed by the complex envelope w 1 of the wrinkles. This leads to the following system: This system is only presented as an example of the multi-scale procedure and it will not be discussed further in this paper.

Two analytical solutions for clamped rectangular membranes
In this Section, the reference macroscopic model (22) (30) (31) (32) is studied by seeking closed-form solutions in the case of clamped rectangular membranes: ω = [0, L X ] × [0, L Y ] as pictured in Figure 3. The membrane is submitted to a large uniform tensile stress N Y = hσ Y > 0 and to a small uniform compressive stress loading N X = hσ X = − λ < 0. If the plate is clamped, it is known [40] that the envelope w vanishes on the boundary. So the corresponding boundary conditions will be Classically, we start with a "linear stability analysis" that is quite simple due to the stress state depending linearly on the applied force. This will establish unambiguously the multi-scale character of membrane wrinkling. Next a classical nonlinear bifurcation analysis will be done, in the same manner as for a classical post-buckling computation [68,69]. This will provide a relation between the applied compressive stress and the size of the wrinkles.

An analytical solution for wrinkling initiation
The linearised version of the envelope Equation (32) is rewritten as Since the envelope vanishes at the boundary, the smallest eigenfunction of (41) is in the form: w(X, Y) = sin(πX/L x )sin(πY/L Y ), as shown in Figure 4. This leads to a classical relation between compressive stress and wavenumber.
Classically, the instability wavenumber Q is chosen by minimizing the stress as a function of the wavenumber. For simplicity, we take into account the orders of magnitude 1 < < QL X , 2DQ 2 /h ≈ |σ X | < < σ Y to simplify (42) in the following manner: The minimum of the latter yields values of the wavenumber and of the critical compressive stress that are consistent with the results of the literature [2,7]: This simple calculation brings out the multiple scale character of the wrinkling phenomenon: indeed, the wrinkling threshold depends on the wavelength that is a microscopic quantity, but this wavelength depends of the width of the plate that is a macroscopic length. Thus a full wrinkling analysis has to associate micro and macro scales.
Recent experimental results [9] have established that the wavelength increases when the thickness h increases and decreases when the tension σ Y increases. The scaling law Y seems qualitatively consistent with experimental results. The wrinkling stress can also be written in terms of the wavenumber by N X = − 2DQ 2 , to be compared with the bifurcation load N X = − DQ 2 of the pure membrane model (34), but of course the pure membrane model is not able to predict the wavenumber.

Initial post-bifurcation analysis
Now we try to connect the applied loads and the size of the wrinkles. As in any plate buckling, the deflection is proportional to the square root of the difference between the current load and its critical value. This post-bifurcation analysis will be done by starting from the new reference membrane model (22) (30) (31) (32) in order to illustrate its ability to characterize the evolution of the wrinkles beyond the instability threshold.
The membrane model (22) (30) (31) (32) is a macroscopic one already simplified by a multi-scale analysis. Therefore there is no need to come back to a multiple scale bifurcation analysis of Ginzburg-Landau type. A standard post-bifurcation analysis as in [68,69] will be sufficient to predict the amplitude of wrinkles.
According to [69], the solution of the symmetric bifurcation problem (22) (30) (31) (32) can be solved by seeking the unknowns and the compression load as Taylor series with respect to a scalar bifurcation parameter "a": To avoid any ambiguity with Fourier expansion (U i ), the terms of the Taylor expansions have been denoted with a superscript U (i) . This straightforward bifurcation analysis leads to Partial Differential Equations at each order.
At order 0, the solution is a membrane state: By introducing the linear operator the PDE's at orders 1, 2, 3 are: Of course the first order Equation (48) corresponds exactly to the bifurcation Equation (41) that has been solved previously. It is not necessary to re-discuss this point. The smallest buckling mode is Let us solve the second order problem (49). The equilibrium Equation (49-a) can be solved by introducing a stress function: The membrane constitutive law (49-b) can be detailed as: The in-plane displacement u (2) ,v (2) can be eliminated from the combination ∂ 2 The resulting equation is: In what follows, this first term will be neglected. To account for the second order term λ (2) in (52-c), one can introduce a new stress function ψ(X, Y) by: so that ψ(X, Y) is solution of the following boundary value problem: The Dirichlet boundary conditions in (55) are deduced from (49-c) (49-d) by integration along the boundary, by a quite well known calculation when dealing with a stress function.
Because of the boundary conditions, it seems difficult to get a closed form solution of (55). We propose to build an approximate equation by using a Galerkin approximation. Let us transform the domain ω = [0, Let us consider a shape function that satisfies the boundary conditions (55-b) (55-c) (55-d), for instance The PDE (55-a) can be rewritten as (with α ¼ L Y L X ): Its Galerkin approximated solution is: The corresponding axial stress field is: Finally we take into account the third order Equation (50-a), but we only need to write an orthogonality condition between the mode w (1) and the r.h.s. of (50-a): ZZ For simplicity we disregard the second term of (59) that is of order O 1 L 2 while the first one is of order O Q 2 À Á ¼ O 1 ℓ 2 (ℓ is the small instability wavelength), which leads to a simpler solvability condition: ZZ This leads to the value of the second order load λ (2) , that is written in terms of four adimensional integrals I 1 , I 2 , I 3 , I 4 This identity (60) permits to evaluate the increment of the load due to the wrinkles. It is convenient to compare this load to the approximated critical load λ (0) = 2DQ 2 : Curiously, this ratio depends rather little on the macroscopic lengths L X , L Y and the wrinkling wavelength ℓ, except via the aspect ratio α ¼ L Y L X . Obviously the critical load depends strongly on ℓ or on L Y by λ 0 ð Þ Eh ≈ h l À Á 2 ≈ h L Y : This means that for a given aspect ratio α ¼ L Y L X , there is an universal law connecting the ratio "amplitude of wrinkles over the thickness" to the ratio "loading over critical wrinkling load", but clearly the critical wrinkling load depends on the traction load. Finally by noting; and by accounting for (44), we obtain the law connecting the wrinkling amplitude as a function of the compressive stress: Note that this analytical formula comes from a structural analysis. Hence the results depend on the domain via the coefficient α ¼ L Y L X and on the boundary conditions. The formula (62) is illustrated by numerical values in the case of a square membrane in Table 1 and Figure 4, for four values of the tensile force. Let us underline that the previous bifurcation analysis depends strongly on the spatial derivatives in the envelope Equations (32) and therefore on the boundary conditions satisfied by the envelope w. For instance, if one considers a simply supported membrane, the macroscopic boundary conditions are of Newman type ( ∂w ∂n ¼ 0 ) and the mode is w (1) = 1 so that the bifurcation analysis from (32) should imply a bifurcation branch with a constant load, as in straight beam buckling. In other terms, the stable postbuckling response obtained in (62) is strongly related to the spatial evolution of the wrinkling mode. The non-linearity of this response is not negligible, as shown by the last column of Table 1 and Figure 5, which means that a significant wrinkling amplitude requires a significant increasing of the compression.

Numerical implementation
We describe now the discretisation techniques and the solving algorithm for the reference model (30) (31) (32). Because this system is a second order partial differential system, any classical C 0 finite element is acceptable. In the applications, quadratic quadrilateral finite elements (Q8) will be used. The resulting nonlinear discrete equations will be solved by the Asymptotic Numerical Method [70]. Some numerical tests will be performed with the pure membrane model (34); its discretisation being straightforward, the details will be omitted.

Matrix form of the bending virtual work
The bending energy (27) depends only on the envelope of the deflection that can be represented by the vector: After introducing the matrix The bending energy can be written as: Matrix form of the membrane virtual work The membrane energy (23) is a quadratic function of the membrane strain γ given in (21) (22) that is a quadratic function of the displacement. The displacement field and its useful derivatives are described by the vector: Note that the two displacement vectors (63) (66) are connected by a transform matrix [T β ]: The strain {γ} can be expressed by a constant matrix [H] and by a matrix [A(θ)] that depends linearly of the displacement vector: With these notations, the membrane constitutive laws and the membrane energy can be written as: Discretisation A classical 2D-Q8 finite element that is defined by eight nodes and three degrees of freedom per node is used. The displacement (u,w) and the full vector on each element are discretised in the form: We need also few matrices to connect the discrete displacement and the membrane strain: Hence we obtain the discrete form of the full internal virtual work: The external virtual work is quite classical and the details will be omitted.

Continuation procedure
The full discrete problem (54) (55) is solved by a continuation procedure called Asymptotic Numerical Method [70] that is a step by step continuation algorithm. In each step, the unknowns are expanded into series with respect to a real path parameter "a": To compute the term of order "n" in (56), we followed a procedure described in numerous papers of the literature [66,70], which is not developed in this paper.

Two numerical solutions
The procedure described in the previous section has been applied to two numerical tests: the wrinkling of a rectangular membrane submitted to a compressive-tensile loading studied previously in section "Clamped rectangular membrane submitted to tension and compression" and the wrinkling of long rectangular membrane submitted to a uniaxial traction that was studied by several authors in the literature. The main question is to evaluate the validity range of the reduced model (30)(31)(32) and (34).

Clamped rectangular membrane submitted to tension and compression
In this section, we study again the example of a clamped rectangular membrane under bi-axial tension-compression load (see Figure 3). The side lengths L X and L Y are respectively 400 mm and 200 mm, and the thickness h is 0.05 mm. The applied tension is N Y = 10 N/mm and the compression force N X increases in the path following procedure. In the macroscopic model, a wavenumber Q has to be chosen and we take the one predicted by the analytic formula (44). The Figure 6 presents the wrinkling pattern just after the bifurcation. Clearly the envelope is nearly sinusoidal in the two directions OX and OY, as predicted by the analytical solutions. It appears also that the envelope model is able to describe correctly the wrinkling shape predicted by the full shell model. A more quantitative picture is provided by Figure 7, where the results of the pure membrane model (34) are also reported. First the finite element study corroborates the analytical bifurcation study. Second the pure membrane model (34) underestimates very much the wrinkling load. Last the macroscopic model (30) (31) (32) is able to describe quite perfectly the initial post-bifurcation response. Let us underline that the macroscopic model requires less degrees of freedom than the full shell model.
As it is well known [43,46] with the Ginzburg-Landau equation, the post-bifurcation profile evolves rapidly. Just after the bifurcation (w/h = 0.16), the envelope is nearly sinusoidal as predicted analytically, cf Figure 8. Next it changes rapidly for a hyperbolic tangent shape (w/h = 0.2, Figure 9), what is also predicted by Ginzburg-Landau theory. The numerical results show that the profile changes gradually and becomes more and more localized, with two peaks at the end of the boundary layers (w/h = 2, Figure 10).
The Figure 11 evaluates the ability of the macroscopic model to predict large sizes of wrinkles. It appears that correct amplitudes can be obtained for rather large wrinkles (w/h = 1). Nevertheless one has to be careful in the analysis of the responses for large deflections because these problems can have many stable and unstable solutions and it is not obvious that a path following method with a full shell model provides always the most relevant solution. In the same Figure 11, we reported the result predicted by the approximate analytical solution (61). It gives correctly the beginning of the bifurcated branch up to about w/h ≈ 0.3.

A rectangular membrane submitted to uniaxial traction
A thin rectangular membrane whose dimensions are h = 0.05 mm, L X = 1400 mm, L Y = 200 mm, is submitted to a uniaxial load (see Figure 12). The long sides are stress-free. Along the short sides, an increasing displacement is applied in the X-direction, the OY-displacement being locked. Linear bifurcation was studied in [56,71] with the same boundary conditions and a nonlinear bifurcation as here in [33], but with a prescribed stress at the boundary. In this case, the wrinkling occurs for rather high axial stresses since it is caused by small transverse compressive stresses due to the boundary effects (for a finer analysis, see [71]). Comparing with the previous tests, it permits to evaluate the ability of the macroscopic model in a case of non uniform pre-buckling stresses.
In Figure 13, one sees that the post-bifurcation patterns obtained by the new reduced model (30) (31) (32) are quite similar as those provided by the full shell model. This  establishes the relevance of this new reduced model to represent the wrinkling modes in a case with a non uniform pre-buckling stress field.
In Figure 14, a response curve is plotted for these two models: the macroscopic model (30) (31) (32) and the full shell model. As in the previous case, the new reduced model gives about the same bifurcation point as the reference model as well as the post-bifurcation response. The number of degrees of freedom is much smaller with the envelope models because they do not need to describe explicitly the full details of the wrinkles. Then a wrinkling pattern along the width is plotted in Figure 15 for the full and the reduced model. Globally they are quite similar but with slight differences for the amplitude and for large sizes of wrinkles (up to w/h = 4). One can wonder if these results could still be improved by keeping a complex envelope as in (36)(37)(38)(39) or by keeping higher order harmonics.

Linear wrinkling analysis revisited by a double scale asymptotic analysis
Last we re-discuss the problem of membrane wrinkling from the point of view of the traditional asymptotic double scale method [72,73] that is slightly different from the technique of Fourier series with slowly variable coefficients. In both methods, one distinguishes two levels of spatial evolution. Within the asymptotic double scale method, the considered problem involves further small parameters and the solutions are also expanded with respect to these small parameters. It will be interesting to compare shortly the two approaches. This will done for the problem of wrinkling initiation in a rectangular membrane 0 < X < L X ,0 < Y < L Y studied previously. A uniform ompression-traction membrane field (N X = hσ X , N Y = hσ Y , σ X < 0 < σ Y , σ Y > > |σ X |) is applied. The wrinkling mode and the wrinkling compressive stress are assumed to solve the linearized Fôppl-Von Karman equation: A small parameter η is introduced that can be identified to the stress ratio: The unknown wrinkling mode is assumed to depend on three independent spatial variables: the starting space coordinates X, Y and another rapidly varying coordinate x ¼ X η in the compressive direction. The eigenpair w(x, X, Y), N X is sought in the form of an asymptotic expansion with respect to η. In agreement with the order of magnitude (81), we define an unknown control parameter λ and another fixed load parameter N −4 Y by: in such a way that λ and N −4 Y will be O(1) in the formal asymptotic expansion. The classical rules of the asymptotic expansion are as follows: If the deflection w and its Y-derivative are assumed to be zero on the sides Y = 0, Y = L Y , the amplitude A(X,Y) must satisfy the boundary conditions ([62] § 2.4): This leads to a first bifurcation pair: If one goes back to the corresponding initial quantities Y hη 4 , one recovers exactly the classical wrinkling stress and wavelength of (44). The wrinkling wavelength is defined at this level by minimizing the critical load λ (2) (q) with respect to the wavenumber. Now Equation (85) corresponding to the second level O 1 η 3 is discussed. Because of (87), (85) is rewritten as: The first term of the right hand side of (91) vanishes because λ (0) = 2Dq 2 . Since the operator L(.) is singular, the existence of a bounded solution x → w (1) (x, X, Y) requires the elimination of "secular terms". This solvability condition leads to λ (1) = 0. Then the solution of (91) has the same shape as the first order term: Last we look at the Equation at the third level (86). By taking into account the previous results, we get The elimination of the secular terms leads to a PDE satisfied by the slowly varying envelope A(X, Y) = B(X)sin(πY/L Y ): The resulting Equation (94) is a one-dimensional eigenvalue problem that permits to define a correction λ (2) η 2 to the critical stress and the variation of the amplitude A(X, Y) = B(X)sin(πY/L Y ) in the compressive direction. In this respect, boundary conditions along the sides X = 0,X = L X are needed. For instance with Dirichlet boundary conditions, the smallest eigenvalue and eigenmode of (94) are: Hence the analysis of the corrective Equation (86) permits to define the evolution of the mode shape in the membrane, especially in the X-direction. This equally permits to correct the buckling stress according to (95-b). With account to this value of λ (2) , the critical stress is corrected as σ wr2 Nonetheless the correction is generally small for very thin membranes, i.e. h L X Þ. This correction depends mainly on the bending effects, as this can be seen from (95) or (96). It depends also on the boundary condition on the four sides, while the first order term σ wr1 X depends only on boundary conditions on X = 0, L.

Last comments about the asymptotic modeling
The asymptotic modeling splits the modal analysis into several steps. First a harmonic fast evolution is obtained at the scale of the wrinkle wavelength, see (87). Then the differential Equation (88) yields a first approximation of the wrinkling load, the wavelength and the variation of the envelope in the transverse direction. Finally the mode is completely characterized if one takes into account its evolution along OX via (94), the latter equation being an eigenvalue problem that provides a small correction λ (2) to the wrinkling load. This splitting of the governing envelope equation corresponds to a consistent account of the different levels of applied stresses and of length scales. This cumbersome multiscale analysis does not exist within the Fourier series that only distinguishes the Fourier coefficients of the response. Unlike the purely asymptotic approach, Fourier method leads to consistent bi-dimensional models as (30) (31) (32) that are easier to handle numerically. Nevertheless the two approaches are not contradictory and an asymptotic procedure could be applied by starting from the reduced models (30-32) (34) or (36)(37)(38)(39).

Conclusion
Membrane wrinkling has been re-discussed from multi-scale approaches. Basically, a local wrinkling wavelength is chosen and macroscopic models are deduced, whose unknowns are envelopes of the fastly varying oscillations of the wrinkles. This amounts to coupling a membrane model or a shell model with a sort of bifurcation equation that governs the appearance of wrinkles. This new class of models leads to well-posed partial differential systems that are easy to be solved numerically. These models have been studied in details, by computing several numerical and analytical solutions. This new approach is rather different from those classically used in the literature, but it permits to recover their main features: a quite different behavior of the membrane with or without wrinkling, a splitting into macroscopic strain and wrinkling strain and the presence of an internal length that proves to be necessary to get a consistent distribution of the wrinkles. A main originality is the ability to predict the main characteristics of the wrinkles from macroscopic models and coarse meshes. Likely the weakest point is the necessity to choose a priori the wavelength in the case of a nonlinear analysis, but this could be improved by coming back to a true multi-scaled approach as in FE 2 method.