# Thermal wrinkling of thin membranes using a Fourier-related double scale approach

Kodjo Attipou, Heng Hu, Foudil Mohri, Michel Potier-Ferry, Salim Belouettar

► To cite this version:

Kodjo Attipou, Heng Hu, Foudil Mohri, Michel Potier-Ferry, Salim Belouettar. Thermal wrinkling of thin membranes using a Fourier-related double scale approach. Thin-Walled Structures, 2015, 94, pp.532-544. 10.1016/j.tws.2015.04.034. hal-01514967

HAL Id: hal-01514967

https://hal.univ-lorraine.fr/hal-01514967v1

Submitted on 10 Mar 2025

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
Thermal wrinkling of thin membranes using
a Fourier-related double scale approach

Kodjo Attipou<sup>a,b,*</sup>, Heng Hu<sup>c</sup>, Foudil Mohri<sup>b</sup>, Michel Potier-Ferry<sup>b</sup>, Salim Belouettar<sup>a</sup>

a Centre de Recherche Public Henri Tudor, 29, avenue John F. Kennedy, L-1855 Luxembourg-Kirchberg, Luxembourg

b Université de Lorraine, Laboratoire d'Étude des Microstructures et de Mécanique des Matériaux (LEM3), UMR CNRS 7239, Ile du Saulcy, F-57045 Metz, France

c School of Civil Engineering, Wuhan University, 8 South Road of East Lake, 430072 Wuhan, PR China

The thermal wrinkling behavior of thin membranes is investigated in this paper. Wrinkles often occur at multiple length scales where induced compressive stresses are located during thermal loading. In the present study, the method of double scale Fourier series is used to deduce the macroscopic membrane wrinkling equations. The obtained equations account for the global and local wrinkling modes. Numerical examples are conducted to assess the validity of the approach developed. The present membrane's model needs only few degrees of freedom to recover more accurately the post-buckling equilibrium curves and the wrinkling shapes. Different parameters such as membrane's aspect ratio, wave number, and pre-stressed membranes are discussed from a numerical point of view and the properties of the wrinkles (critical load, wavelength, size, and location) are presented.

# 1. Introduction

Membrane structures are widely used in many engineering fields such as aeronautic, space, automotive and civil construction, as well as skydiving. For instance, gossamer structures are used in the space industry for satellites and spacecrafts. In civil construction, thin membranes such as PVC or aluminum are used in roof covering of large space area as in sport area and exhibition halls. Aeronautical industry uses thin polymer membranes to strengthen airplane wings while keeping them more lighter. Thin membranes are often subjected to stresses and strains due to their applications. For very thin membranes, it can lead to an out-of-plane deformation of the structures because of no flexural stiffness. In some applications, the membranes can be exposed to high temperature changes coupled with mechanical tensile and/or compressive stresses (satellites, PVC manufacturing, etc.). In relatively large thin membranes, thermal buckling occurs locally and repetitively and results into wrinkles. The buckling or wrinkling phenomena are instabilities that can damage the membrane structures when they are not controlled. Controlling the buckling phenomenon means assessing its occurrence conditions in order

to predict and simulate it. Buckling and wrinkling are caused by compressive stresses. They are all boundary conditions and imperfections dependent. Buckling modes involve length scale in the same order as the size of the structure. Wrinkling is observed especially in membrane structures with a modal length scale much smaller than the membrane's dimensions. This last instability is a complex non-linear phenomenon and its modeling is an interesting and up-to-date research topic that deserves attention. The scope of this study is the wrinkling of thin membranes under thermal loading.

The buckling and wrinkling phenomena in membrane have been studied in the literature. Among these papers, many are focused on the wrinkling of membrane caused by residual stresses. In [1], Cuong et al. have applied thermal stresses on a thin membrane to simulate the residual stresses that appear during the rolling process of a sheet in manufacturing. Residual stresses can also be induced by the growth of living tissues such as mucosa [2]. The most common case where residual stresses are induced is the case of thin membranes subjected to mechanical loads. For instance, Kim et al. [3] have conducted a modified Yoshida buckling test to investigate the wrinkling initiation caused by the mechanical properties as well as the loading case.

There are several techniques to study the wrinkling phenomena in thin membranes: experimental, analytical and numerical. Benchmarks tests have been developed by many authors in the literature. Blandino et al. [4] have proven through experimental

<footnote>* Corresponding author at: Université de Lorraine, Laboratoire d'Étude des Microstructures et de Mécanique des Matériaux (LEM3), UMR CNRS 7239, Ile du Saulcy, F-57045 Metz Cedex 01, France.</footnote>

<footnote>E-mail address: kodjo.attipou@yahoo.fr (K. Attipou).</footnote>
results that the number of wrinkles increases with the mechanical loads in isothermal conditions. Moreover, the wrinkles amplitude increases with the thermal loading. Besides benchmarks tests, analytical and/or numerical studies have also been conducted for comparison purpose [5,6]. The aim is to develop reliable methods to solve the wrinkling problems. These methods of resolution can be sort into two main groups: the full shell methods and the reduced methods. Full shell methods involve classical equations that are resolved via finite element analysis using shell elements [7]. Such microscopic models are simply based on elastic shell theory (see for instance [8–12]). The cons of the direct methods is that they required many degrees of freedom, then the computation cost becomes important. This is the reason why some alternative methods as the reduced methods are interesting where it applies. Some reduced methods are based on the modification of the constitutive equations of the wrinkling problems. It follows a simplification of the general equations and then a simplified model is built. In 1987, Roddeman et al. [13,14] have developed a wrinkling model for thin membranes under shear load by using a weighted residual method. Roddeman's model was just limited to the determination of wrinkled regions. No information on the wrinkles characteristics (size, wavelength) are available. Chiu et al. [15] have modified Roddeman's model to account for the thermal effect in thermo-mechanically loaded polymer membranes.

Recently, a new reduced method was introduced by Damil et Potier-Ferry [16,17] to analyze the wrinkling phenomenon in one-dimensional beam model. This method is based on a Fourier-related multiscale approach and it is able to couple the macroscopic response of the structure with the characteristics of the wrinkles. With this approach, it is possible to model the wrinkling behavior of a structure around the bifurcation point and away from where former techniques such as Landau-Ginzburg are no longer valid. More precisely the Landau-Ginzburg approach relies on two asymptotic assumptions: closeness to the bifurcation point and separation of scales. The method described in [16,17] is only based on the separation of scales and a truncation of higher Fourier modes and thus it can remain valid away from the bifurcation, what is confirmed by numerical tests in the literature. The other advantage of the Fourier method is that it requires only few degrees of freedom to recover the wrinkling behavior. Of course the advantages of such a reduced-order model are always counterbalanced by some more or less restrictive assumptions, the main one being here the need of choosing a priori the wavenumber of the wrinkling mode. Applications of Fourier method to sandwich structures were conducted by Liu et al. [18] and Yu et al. [19]. Our present study is based on the work of Damil et al. in [20,21] where the Fourier method has been successfully extended to two-dimensional membrane models. In our work, the thermal effect in thin membrane models is investigated. Thermal stresses are introduced in the classical equations within the framework of the elastic shell theory. The nonlinear problem is solved using an Asymptotic Numerical Method described in [22] for efficiency and reliability purpose.

This paper is divided into three main sections. In Section 2, the equilibrium equations are developed. Classical membrane model equations are recalled within the approximations of von Kàrmàn using the classical plate theory of Love-Kirchhoff. Then macroscopic membrane model's equations are deduced using the Fourier series. The generalized Duhamel-Neumann form of the Hooke's law is considered to account for the thermal stresses. In Section 3, an analytical study is performed to investigate the wrinkling critical load and wavenumber in some specific boundary and load conditions. Section 4 is dedicated to numerical examples where a validation test is performed and different parameters are studied to assess the reliability and efficiency of the Fourier series technique developed. The results of the examples come from the

implementation of the new membrane model equations in MATLAB code and the simulation of the classical membrane model in the finite element software ABAQUS for comparison.

## 2. Nonlinear equilibrium equations

### 2.1. Full model equations

The main objective of this paper is to develop a nonlinear macroscopic model of thin membranes subjected to thermomechanical loads. In this framework, the related equations of the macroscopic models are derived from the Love-Kirchhoff classical plate theory within von Kàrmàn approximations. An elastic isotropic membrane is considered as the reference model. The generalized Duhamel-Neumann form of the Hooke's law is applied to account for the thermal stresses. One recovers the well-known governing equations of the full membrane model as following:

$$ \mathrm{div} \, \mathbf{N} = 0 \qquad (1a) $$

$$ \mathbf{N} = \mathbf{L}^m \cdot (\boldsymbol{\gamma} - \boldsymbol{\alpha} \Delta T) \qquad (1b) $$

$$ 2\boldsymbol{\gamma} = \nabla \mathbf{u} + t \nabla \mathbf{u} + \nabla w \otimes \nabla w \qquad (1c) $$

$$ \mathbf{D} \nabla^4 w - \mathrm{div}(\mathbf{N} \cdot \nabla w) = 0 \qquad (1d) $$

where $\mathbf{u}=(u, v) \in \mathbb{R}^2$ is the in-plane displacement of the microscopic fields, $w$ the deflection, $\mathbf{N} \to t(N_x N_y N_{xy})$ and $\boldsymbol{\gamma} \to t(\gamma_x \gamma_y 2\gamma_{xy})$ are the membrane stresses and strains. The membrane and bending elasticity tensors are represented respectively by the matrices $\mathbf{L}^m$ and $\mathbf{D}$ in case of a single-layered membrane and its mid-plane set as reference of the displacement components:

$$ \mathbf{L}^m = \frac{Et}{1-\nu^2} \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1-\nu}{2} \end{bmatrix} \qquad (2) $$

$$ \mathbf{D} = \frac{Et^3}{12(1-\nu^2)} \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1-\nu}{2} \end{bmatrix} \qquad (3) $$

with $t$ being the membrane thickness, $E$ the Young's modulus and $\nu$ the Poisson ratio. The thermal properties of the membrane are represented by the expansion coefficient $\boldsymbol{\alpha} \to t(\alpha_x \alpha_y \alpha_{xy})$ and the temperature variation $\Delta T$. The material constants $(E, \nu, \boldsymbol{\alpha})$ are assumed to be independent towards the temperature change.

### 2.2. A simplified macroscopic model

The method of Fourier series consists in replacing the unknowns of the nonlinear full model (Eq. (1)) by the Fourier series which have slowly variable coefficients. Nonlinear macroscopic equations whose unknowns are the Fourier coefficients are recovered. The discussion in this paper is devoted to two-dimensional structures. The reader should refer to the literature [16, 17] for the one-dimensional development. First in our model, the wrinkles patterns are assumed to be in one single direction (the axial $x$ direction) and the wavenumber is assumed to be known. These assumptions are restrictive and limited to simplest case of membrane wrinkling and should be removed in future works when dealing with shear membrane wrinkling. Hence, to be more realistic, a discussion is first made around a multiscale approach with a given wavelength and direction of wrinkles. The local problem unknowns $U$ are sought in Fourier series terms in
the following form:

$$ \mathbf{U}(x, y) = \sum_{m=-\infty}^{+\infty} \mathbf{U}_m(x, y) \exp(imqx) \quad (4) $$

with $i^2 = -1$ being the imaginary number and $q$ the wavenumber. The number of harmonics is represented by the subscript $m$. The Fourier coefficients $\mathbf{U}_m(x, y)$ represent the unknowns of the macroscopic fields and vary slowly on a single period $[x, x+(2\pi/q)]$. The derivation operators are calculated exactly according to the following rules ($\mathbf{U}$ and $\mathbf{V}$ are of Fourier series type):

$$ \left(\frac{d^n \mathbf{U}}{dx^n}\right)_m = \left(\frac{d}{dx} + imq\right)^n \mathbf{U}_m \quad (5) $$

$$ \left(\frac{d^n \mathbf{U}}{dy^n}\right)_m = \left(\frac{d}{dy}\right)^n \mathbf{U}_m \quad (6) $$

$$ (\mathbf{U} \cdot \mathbf{V})_m = \sum_{k=-\infty}^{+\infty} \mathbf{U}_k \cdot \mathbf{V}_{m-k} \quad (7) $$

where the superscript $n$ stands for the $n$th derivative of the associated function.

To recover the macroscopic version of the full membrane equations, the method of Fourier series is applied to the equations of the system (Eq. (1)) and a simple identification of the Fourier coefficients follows. The unknowns fields of the previous system that are sought in terms of Fourier series are $\mathbf{U}(x,y) = (\mathbf{u}(x,y), w(x,y), \gamma(x,y), \mathbf{N}(x,y))$. Only the components of $\mathbf{U}$ are discretized.

As depicted in Fig. 1, at least two macroscopic fields are needed to describe nearly periodic oscillations: a mean value and an amplitude of the oscillations. The first envelope is a real-valued (corresponding to harmonic 0) while the next ones are complex. For simplicity, the study in this paper is limited to two harmonics $\mathbf{U}_0(x,y)$, $\mathbf{U}_1(x,y)\exp(iqx)$ and $\mathbf{U}_{-1}(x,y)\exp(-iqx)$. The mean field $\mathbf{U}_0$ is recovered by the harmonic 0 and the first order harmonics $(\mathbf{U}_1, \mathbf{U}_{-1})$ represent the amplitude of the oscillation.

Some approximations are made to stay within the same framework as in Landau-Ginzburg asymptotic approach. For instance, the amplitudes are taken slowly variable and their variations are neglected: only the deflection $w_1$ is kept and all its derivatives of order four are dropped as well as the complex values. The membrane displacements and stresses are restricted to their non-oscillating part. It follows the definition of the new macroscopic unknowns fields:

$$ (u, v, w, \gamma, \mathbf{N}) \stackrel{\text{def}}{=} (u_0, v_0, w_1, \gamma_0, \mathbf{N}_0) \quad (8) $$

According to these principles, the macroscopic strain equation is deduced from Eq. (1c) and is written as follows:

$$ \{\gamma\} = \left\{ \begin{array}{l} \frac{\partial u}{\partial x} + \left(\frac{\partial w}{\partial x}\right)^2 + q^2 w^2 \\ \frac{\partial v}{\partial y} + \left(\frac{\partial w}{\partial y}\right)^2 \\ \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} + 2 \frac{\partial w}{\partial x} \frac{\partial w}{\partial y} \end{array} \right\} \quad (9) $$

The latter equation is similar to the classical well-known strain-displacement relation for the full model equations. The only relevant

term is the additional term $q^2 w^2$ that depends on the wavenumber $q$ in a quadratic form. This term appears in the $x$-direction as expected because of the assumption taken about the direction of the wrinkles patterns formation. It is also included in the nonlinear term of the strain-displacement relation and hence contributes for the out-of-plane displacement of the membrane. The linear and nonlinear terms of Eq. (9), defined respectively $\boldsymbol{\varepsilon}(\mathbf{u})$ and $\boldsymbol{\gamma}^{\text{wr}}(w)$ are written as follows:

$$ \{\boldsymbol{\varepsilon}(\mathbf{u})\} = \left\{ \begin{array}{l} \frac{\partial u}{\partial x} \\ \frac{\partial v}{\partial y} \\ \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \end{array} \right\} \quad (10) $$

$$ \{\boldsymbol{\gamma}^{\text{wr}}(w)\} = \left\{ \begin{array}{l} \left(\frac{\partial w}{\partial x}\right)^2 + q^2 w^2 \\ \left(\frac{\partial w}{\partial y}\right)^2 \\ 2 \frac{\partial w}{\partial x} \frac{\partial w}{\partial y} \end{array} \right\} \quad (11) $$

From the full membrane bending Eq. (1d) and by introducing the Fourier series definition, the macroscopic bending equation is recovered as follows:

$$ -6\mathbf{D}q^2 \frac{\partial^2 w}{\partial x^2} - 2\mathbf{D}q^2 \frac{\partial^2 w}{\partial y^2} + (\mathbf{D}q^4 + N_x q^2)w - \operatorname{div}(\mathbf{N} \cdot \nabla w) = 0 \quad (12) $$

The macroscopic version of the mechanical stresses includes the thermal strain components which is not sought in terms of Fourier series:

$$ \mathbf{N} = \mathbf{L}^m \cdot [\boldsymbol{\varepsilon}(\mathbf{u}) - \alpha \Delta T + \boldsymbol{\gamma}^{\text{wr}}(w)] \quad (13) $$

We admit that the temperature in the membrane thickness is constant. For this, the thermal strain contributes only to the in-plane displacements. The thermal stress applied to the membrane is written as

$$ \mathbf{N}^{\text{th}} = \mathbf{L}^m \cdot \alpha \Delta T \quad (14) $$

Finally the macroscopic system equations are

$$ \left\{ \begin{array}{l} \operatorname{div} \mathbf{N} = 0 \\ \mathbf{N} = \mathbf{L}^m \cdot \boldsymbol{\gamma} - \mathbf{N}^{\text{th}} \\ \boldsymbol{\gamma} = \boldsymbol{\varepsilon}(\mathbf{u}) + \boldsymbol{\gamma}^{\text{wr}}(w) \\ -6\mathbf{D}q^2 \frac{\partial^2 w}{\partial x^2} - 2\mathbf{D}q^2 \frac{\partial^2 w}{\partial y^2} + (\mathbf{D}q^4 + N_x q^2)w - \operatorname{div}(\mathbf{N} \cdot \nabla w) = 0 \end{array} \right. \quad (15) $$

When all the membrane's edges are restricted, only thermal stresses remain. A temperature raise would lead to compressive stresses that may cause global buckling of the membrane. In this paper, one wishes to investigate the wrinkling phenomenon in membrane under thermal conditions. Some special boundary conditions have been considered for such a study. For instance, boundary conditions consisting in two restrained opposite edges and a temperature drop induces localized compressive stresses (because to Poisson's ratio) in the direction perpendicular to the restrained edges. Detailed examples are presented in Section 4 for illustration. The next sections are dedicated to the resolution of the macroscopic system. An analytical solution is sought for the linearized version of the system in order to establish the multiple-scale character of the wrinkling phenomenon. Then a numerical analysis is conducted to solve the nonlinear problem (Eq. (15)) to show that this macroscopic model is able to describe the evolution of the wrinkles even with few degrees of freedom than the full model.

<figure><img src="image_2.png" /><figcaption>Fig. 1. At least two macroscopic fields needed to describe the oscillations.</figcaption></figure>
### 3. Analytical analysis

To solve analytically the system represented by Eq. (15), one must a priori determine the wavenumber $q$. Let us consider a rectangular membrane $L \times B$, depicted in Fig. 2, that is subjected to an uniform temperature drop $\Delta T < 0$ and a small uniform compressive stress loading $N_x = -|N_x| < 0$. As specific boundary conditions, the two edges at $x=0$ and $x=L$ are restrained in the

<figure><img src="image_3.png" /><figcaption>Fig. 2. Membrane dimensions and boundary conditions type for biaxial loading.</figcaption></figure>

<table><caption>**Table 1**
Membrane dimensions and material properties.</caption><thead><tr><th colspan="2">Dimensions and properties</th></tr></thead><tbody><tr><td>Length</td><td>L = 400 mm</td></tr><tr><td>Width</td><td>B = 200 mm</td></tr><tr><td>Thickness</td><td>t = 0.05 mm</td></tr><tr><td>Young's modulus</td><td>E = 70,000 N/mm<sup>2</sup></td></tr><tr><td>Poisson ratio</td><td>ν = 0.3</td></tr><tr><td>Thermal expansion</td><td>α = 23.10<sup>-6</sup>/K at 20 °C</td></tr></tbody></table>

<table><caption>**Table 2**
Biaxial load results table.</caption><thead><tr><th rowspan="2">Temperature ΔT (°C)</th><th rowspan="2">Wavelength l<sup>ωr</sup> (mm)</th><th colspan="2">Critical force σ<sub>x</sub><sup>ωr</sup> (N/mm<sup>2</sup>)</th><th rowspan="2">Error(%)</th></tr><tr><th>Analytical</th><th>Numerical</th></tr></thead><tbody><tr><td>-1</td><td>89.0466</td><td>0.1596</td><td>0.1640</td><td>2.75</td></tr><tr><td>-5</td><td>59.5490</td><td>0.3568</td><td>0.3600</td><td>0.89</td></tr><tr><td>-10</td><td>50.0746</td><td>0.5046</td><td>0.5075</td><td>0.57</td></tr><tr><td>-15</td><td>45.2475</td><td>0.6180</td><td>0.6201</td><td>0.34</td></tr><tr><td>-20</td><td>42.1075</td><td>0.7137</td><td>0.7150</td><td>0.18</td></tr><tr><td>-30</td><td>38.0485</td><td>0.8740</td><td>0.8742</td><td>0.02</td></tr></tbody></table>

y-direction. The thermal stress induced by the temperature drop is in the y-direction and is written as

$$N_y = -Et\alpha_y\Delta T = Et\alpha_y|\Delta T| \quad (16)$$

This stress is uniform throughout the membrane because of the free lateral displacement of the restrained edges. The instability wavelength will be deduced from the envelope equation (12) that is rewritten in a linearized version as

$$-6\mathbf{D}q^2\frac{\partial^2 w}{\partial x^2} - (2\mathbf{D}q^2 + N_y)\frac{\partial^2 w}{\partial y^2} + Dq^4w = |N_x| \left( q^2w - \frac{\partial^2 w}{\partial x^2} \right) \quad (17)$$

Since the boundary conditions type of all four edges is clamped-clamped, the deflection $w$ and its first spatial derivatives are assumed to be zero on the boundary, what leads to a linear mode in the form

$$w(x,y) = \sin\left(\frac{\pi x}{L}\right) \sin\left(\frac{\pi y}{B}\right) \quad (18)$$

This leads to a classical relation between compressive stress and wavenumber:

$$|N_x| = \frac{6\mathbf{D}q^2\frac{\pi^2}{L^2} + (2\mathbf{D}q^2 + N_y)\frac{\pi^2}{B^2} + \mathbf{D}q^4}{q^2 + \frac{\pi^2}{L^2}} \quad (19)$$

The wrinkling wavenumber is then determined by minimizing the stress as a function of $q$. For the sake of simplicity, orders of magnitude are taken into account to simplify Eq. (19). One assumes that the membrane dimensions are sufficiently large regarding the thickness to have the following statements: $1/L \ll q$, $2\mathbf{D}q^2 \approx |N_x| \ll N_y$. By applying those simplifications to Eq. (19), one obtains

$$|N_x(q)| = \frac{N_y\pi^2}{q^2B^2} + \mathbf{D}q^2 \quad (20)$$

<figure><img src="image_4.png" /><figcaption>Fig. 4. Membrane's boundary conditions and loading case.</figcaption></figure>

$$\sigma_x = -0.98806 \text{ N/mm}^2$$

<figure><img src="image_5.png" /><figcaption>Fig. 3. Membrane wrinkling patterns for $\Delta T = -30$ °C.</figcaption></figure>
By replacing Eq. (16) in Eq. (20), the wavenumber and the critical compressive stress, respectively $q^{wr}$ and $\sigma_x^{wr}$, are deduced from the minimum of the latter (Eq. (20)):

$$q^{wr} = \left( \frac{12(1 - \nu^2)\pi^2 \alpha_y}{t^2 B^2} \middle| \Delta T \right)^{1/4} \quad (21)$$

$$|\sigma_x^{wr}| = \pi E \left( \frac{\alpha_y}{3(1-\nu^2)} \right)^{1/2} (|\Delta T|)^{1/2} \frac{t}{B} \quad (22)$$

The wavelength is deduced from Eq. (21):

$$l^{wr} = \frac{2\pi}{q^{wr}} = \left( \frac{4\pi^2}{3(1-\nu^2)\alpha_y} \right)^{1/4} (tB)^{1/2} \left( \frac{1}{|\Delta T|} \right)^{1/4} \quad (23)$$

The multiple-scale character of the wrinkling phenomenon can be deduced from Eq. (23). Indeed, the microscopic quantity is represented by the wavelength $l^{wr}$ and the macroscopic quantity is represented by the membrane dimension term $(tB)^{1/2}$. Some comments can be made regarding Eqs. (22) and (23):

(a) When the temperature variation increases, the wavelength decreases.

(b) The membrane dimension $L$ has no effect either on the wavelength or the critical compressive stress; only the number of oscillations changes with $L$.

## 4. Numerical examples

In this section, several examples are presented for the wrinkling analysis of a thin membrane under thermal and thermo-mechanical loads. These examples are characterizations of the wrinkling behavior of the membrane. The scope is to assess the efficiency and the reliability of the Fourier series method developed in this paper. The macroscopic equations for the reduced membrane model are implemented in MATLAB code and resolved numerically using the Asymptotic Numerical Method for efficiency purpose. For validation and reference purposes, the full membrane model is simulated using commercial finite element code ABAQUS. In ABAQUS/Standard, only the full model equations are considered and the nonlinear equilibrium equations are solved according to the arc-length method (modified Riks method). The full membrane model studied using ABAQUS/Standard has been discretized by eight-node quadratic shell elements with reduced integration (S8R5). The reduced membrane model has been discretized by quadratic quadrilateral elements implemented in MATLAB (Q8 element-type).

Two main examples are conducted in this paper. The first example is related to the analytical resolution discussed in earlier (Section 3). For this example, a biaxial loading case is investigated. The second example is devoted to the analysis of a membrane under uniaxial loading case. It is characterized by a pure thermal

loading with specific boundary conditions that results in an uniaxial tensile stress.

### 4.1. An example of validation: biaxial loading

In this example, the numerical analysis of a thin membrane under biaxial loading case is conducted. This is an example of validation towards the analytical resolution of the wrinkling problem developed in Section 3. A fixed uniform temperature drop is applied in a first step followed by an increasing compressive force in a second step (Fig. 2). The objective here is to assess the validity of the Fourier double scale series method towards the analytical solution. The membrane dimensions and material properties are displayed in Table 1.

Different values of the temperature variation are set in the first step and the critical compressive force is computed in the second step. The results are compared to analytical ones. There will be no discussion about the mesh since the mesh convergence is discussed deeply in the next example. Nevertheless, a mesh convergence study has been conducted prior to this example to determine the optimized mesh to obtain accurate results.

The results of the biaxial simulation loading case are displayed in Table 2. The wavenumber used in the numerical simulation (Fourier method) is calculated from the analytical point of view defined by Eq. (23). An error percentage of the analytical solution on the numerical one is also deduced. Fig. 3 shows the wrinkles shape for a fixed temperature at a given compressive stress.

It appears that for high temperature variations, both numerical and analytical solutions are very close. However, for low temperature variation, the error is sufficiently small to be negligible.

### 4.2. Uniaxial loading

In this example, a thin isotropic membrane is subjected to a globally uniform temperature decrease. Table 1 displays the dimensions and properties of the membrane studied. In a general case, with all four-edges simply supported boundary conditions, such a thermal loading would lead to global buckling of the membrane [23]. To stay within the framework of membrane wrinkling, specific boundary conditions should be considered. Two opposite edges of the membrane are fully constrained in all

<table><caption>**Table 3**

Different meshes for both membrane models.</caption><thead><tr><th>n<sub>y</sub> × n<sub>x</sub> (Full model)</th><th>n<sub>y</sub> × n<sub>x</sub> (Reduced model)</th></tr></thead><tbody><tr><td>20 × 50 (15,705 dof)</td><td>4 × 10 (447 dof)</td></tr><tr><td>40 × 100 (61,405 dof)</td><td>6 × 15 (939 dof)</td></tr><tr><td>50 × 125 (95,505 dof)</td><td>8 × 20 (1611 dof)</td></tr><tr><td>100 × 250 (378,505 dof)</td><td>12 × 30 (3495 dof)</td></tr><tr><td>-</td><td>20 × 50 (9423 dof)</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Main stresses distribution near bifurcation point (full model simulation). (a) Global tensile stresses, $S_{xx}$. (b) Tensile and compressive stresses, $S_{yy}$.</figcaption></figure>
displacements. A uniform temperature decrease leads the membrane in a tensile state along the axis perpendicular to the restrained edges (Fig. 4). Small compressive stresses are induced because of restrained edges and Poisson's ratio and then cause wrinkles in the membrane (Fig. 5). Several tests are conducted to study the wrinkling phenomenon. The first one is a mesh convergence study to see the effect of meshing on the critical wrinkling temperature and wrinkles patterns. Further examples

involve parametric study such as preloaded membrane in traction and geometric aspect ratio study.

### 4.2.1. Mesh convergence study

The objective of this test is to assess the mesh effect on the membrane wrinkling and therefore set an optimized mesh for the next examples. The mesh convergence study involves the full

<figure><img src="image_7.png" /><figcaption>Fig. 6. Full model's bifurcation curves for mesh convergence study.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Reduced model's bifurcation curves for mesh convergence study.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Wrinkling profiles of the full model at $\Delta T = -10$ °C (mesh convergence study).</figcaption></figure>
membrane and reduced models. Different meshes are considered in this test. The membrane aspect ratio is kept constant. Let $n_x$ be the number of elements in the x direction and $n_y$ the number of elements in the y direction for the two models studied. The following statement is assumed: $n_y/n_x = 0.4$ in all cases. Table 3 summarizes the different meshes for each analysis with the corresponding degrees of freedom (dof).

The results presented in Figs. 6 and 7 show the bifurcation curves of the temperature variation against non-dimensional term

of the maximum deflection (the amplitude over the thickness: w/t). The difference between the critical wrinkling temperature variation is very negligible (≈ 0.3%). It appears that only few degrees of freedom for the reduced model are needed to recover the bifurcation load path while a dense mesh is required for the full model to get accurate results. All the considered meshes in the reduced model are able to describe the bifurcation load path and the bifurcation point is determined accurately. The wrinkling profiles of the membrane at $x=L/2$ for a fixed temperature

<figure><img src="image_10.png" /><figcaption>Fig. 9. Wrinkling profiles of the reduced model at $\Delta T = -10$ °C (mesh convergence study).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Comparison of the bifurcation curves.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Comparison of the wrinkling profiles.</figcaption></figure>
variation (ΔT = −10 °C) near the bifurcation point are displayed in Figs. 8 and 9. One can see that the coarse mesh 4 × 10 in the reduced model is not suitable to recover the exact wrinkles patterns near the boundaries. However, the others meshes are very close. According to the results obtained, the meshes 40 × 100 (full model) and 8 × 20 (reduced model) give a good optimization time and accurate results. These meshes are then considered in the next calculations.

Figs. 10 and 11 are examples of validation of the Fourier series method towards the numerical resolution of the wrinkling problem using the full membrane model. The same critical temperature variation (≈ 9.3 °C) is obtained for both models. The wrinkles profiles are quite the same in the middle region of the membrane. The mismatch observed near the boundaries can be explained by the non-uniform wavelength in the full model, while in the reduced model, the wavelength is set prior to the analysis and remains constant during the computation. The wrinkles shapes of both models are quite similar (Fig. 12).

One advantage of the Fourier method is assessed in these results by considering the matter of the degrees of freedom. The reduced model requires less degrees of freedom than the full model. It follows a non-negligible gain of computation time and cost. Indeed, the reduced membrane model has required at least 30 times less degrees of freedom than the full model (Fig. 10). Then, the computation time for the reduced model in MATLAB is 364 s while the full model in ABAQUS takes 2460 s to complete.

### 4.2.2. Minimal wavelength

In this test, the wavenumber q is discussed from a numerical point of view. The Fourier series definition in Eq. (4) imposed that the determination of the wavenumber should be set prior to the

bifurcation analysis in the reduced model. Different values of the wavelength are set and a numerical simulation is run for each value of the wavelength. It appears a minimum wavelength through the multiple calculations conducted. The relation between the wavenumber defined in Eq. (4) and the wavelength lq is as follows:

$$l_q = \frac{2\pi}{q} \qquad (24)$$

An interval of the wavelength is set according to the width of the membrane and based on an assumption about the approximate wavenumber expected. In the present case, according to the membrane dimensions and properties (Table 1), the approximate interval in which the minimal wavelength is assumed to belong to is set from 50 mm to 60 mm. The interval step is chosen to be equal to 1 mm in order to get results rounded to unity. All the curves are not displayed in this paper for the sake of presentation and clarity. Yet, in Fig. 13, the bifurcation curves of 6 specific values of the wavelength varying from 55 mm to 59 mm are displayed. The bifurcation point decreases from the value of 51 to 58 and then starts increasing from 59.

The curves show that the minimum value of the bifurcation point (critical temperature variation) is obtained for lq = 58 mm. It appears that the minimum wavelength determined for the reduced model corresponds to the varying wavelength of the full model close to the bifurcation point. Indeed, in the full model, the wavelength varies along the load path continuously (Fig. 14(a) and (b)). Then, at the bifurcation point, the wavelength observed in the full model gives the minimum temperature variation in the reduced model calculation.

As said at the beginning of this example, a prediction of the wavelength is required in the reduced model calculations for

<figure><img src="image_13.png" /><figcaption>Fig. 12. Wrinkling shapes for both models at the same load point. (a) Full model membrane at ΔT = −10 °C. (b) Reduced model membrane at ΔT = −10 °C.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Bifurcation curves: wavelength variation (zoom around the bifurcation point).</figcaption></figure>
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Wrinkles evolution with respect to the loading temperature variation (full model). (a) Decay between two temperature variations. (b) Different temperature variations.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 15.** Membrane's boundary conditions in pretensioning step.</figcaption></figure>

uniaxial loading. Several methods can be used. Similar cases are studied by Friedl et al. [24] and Jacques and Potier-Ferry [25]. In those references, an approximated analytical solution is investigated for the wavenumber prior to the buckling analysis of a stretched strip mechanically loaded. Another method is to conduct standard nonlinear analysis (using ABAQUS for example) to find the critical wavelength. This latter method has been used in the present study. The use of this second method could be redundant since it already gives the results needed. However it allows to see the differences in use of the full model and the reduced model. There is a third method used by Wong and Pellegrino [7] which consists in conducting an eigenvalue analysis prior to the

<table><caption>**Table 4**
Pretensioning cases.</caption><thead><tr><td>Displacement <u>u</u> (mm)</td><td>Minimal wavelength (mm)</td></tr></thead><tbody><tr><td>0.01</td><td>54</td></tr><tr><td>0.05</td><td>47</td></tr><tr><td>0.10</td><td>43</td></tr></tbody></table>

bifurcation analysis. Then, the wavelength of the first mode corresponding to the lowest energy will be used in the bifurcation analysis. Since the objective here is to assess the Fourier method efficiency, the minimum wavelength study will not be fully developed. We limited ourselves to simple cases.

### 4.2.3. Pre-stressed membrane

Another parametric example to study the thermal wrinkling phenomenon is the case of a mechanically pre-loaded membrane. In this test, the effect of the temperature variation on a pre-stressed membrane is investigated. The membrane is mechanically loaded in tension prior to the thermal loading. Fig. 15 shows the boundary conditions in the first step (step of pretensioning). The tension is applied by imposed displacement $\bar{u}$ ($\bar{u} > 0$).
Three cases in the pretensioning step are considered: $\bar{u} = \{0.01, 0.05, 0.1\}$. For each case, a prior computation is conducted in the same framework as previous test to determine the minimal wavelength. This is done to get a bifurcation curve corresponding to the minimal wavelength in each case. Table 4 displays the value of $\bar{u}$ and its corresponding minimal wavelength.

In Figs. 16–18, a comparison between the full model and the reduced model results are showed. The legend shows the wavelength in each model. $L_q^{min}$ corresponds to the wavelength required

to get the minimal critical temperature load in the reduced model while $l_q$ corresponds to the wavelength observed from the full shell model simulation. Figs. 19 and 20 display the evolutions of the bifurcation curves and the associated wrinkling profiles in the post-bifurcation state only for the reduced model (with $\Delta T = -20$ °C). The results obtained are in a good agreement of what expected. When the membrane is pre-stressed, its stiffness increases and the bifurcation point as well (Fig. 19). However, when the imposed displacement is too high, it can lead the

<figure><img src="image_18.png" /><figcaption>Fig. 16. Bifurcation curves for $\bar{u} = 0.01$ mm (full and reduced model).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 17. Bifurcation curves for $\bar{u} = 0.05$ mm (full and reduced model).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 18. Bifurcation curves for $\bar{u} = 0.1$ mm (full and reduced model).</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Fig. 19. Evolution of bifurcation curves for the reduced model.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 20. Evolution of wrinkles profiles at $\Delta T = -20$ °C for the reduced model.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 21. Bifurcation curves for $\bar{u} = 0.5$ mm (variation of the wavelength).</figcaption></figure>

membrane into a state close to the bifurcation point. Indeed, for the membrane studied, a pure mechanical analysis has been performed in [24], that consists in applying traction on the two short edges until the wrinkling occurs. The critical stress calculated is $\sigma_{crit} = 236$ N/m$^2$, corresponding to a displacement of $\bar{u} = 0.6743$ mm. Hence, for our pretensioning test, a displacement of $\bar{u} = 0.5$ mm is applied in the first step. The smallest critical temperature load is obtained for a value of $l_q = 30$ mm (Fig. 21). The bifurcation point decreases where it was expected to increase

according to Fig. 19. This can be explained by the fact that a displacement of 0.5 mm puts the membrane into a state close to the bifurcation point where the critical displacement is 0.6743 mm.

#### 4.2.4. Geometry aspect ratio

The effects of the geometry aspect ratio on the membrane wrinkling behavior are assessed in this test. The ratio between the
<figure><img src="image_24.png" /><figcaption>**Fig. 22.** Wrinkling mode for both models at the same load point (r=3). (a) Full model (ΔT = -20 °C). (b) Reduced model (ΔT = -20 °C).</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Fig. 23.** Wrinkling mode for both models at the same load point (r=5). (a) Full model (ΔT = -30 °C). (b) Reduced model (ΔT = -30 °C).</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Fig. 24.** Wrinkling mode for both models at the same load point (r=7). (a) Full model (ΔT = -30 °C). (b) Reduced model (ΔT = -30 °C).</figcaption></figure>

membrane's length and width is written as $r=L/B$. The width is fixed at $B=200$ mm and the length varies according to the ratio value. Three values of the ratio are studied: $r = \{3, 5, 7\}$ which implies $L = \{600, 1000, 1400\}$ (in mm).

According to the results obtained, the phenomenon observed is a variable location of the wrinkled patterns as the ratio increases. Indeed, referred to Figs. 22–24, it appears that the wrinkles move from the middle of the membrane to both ends symmetrically. This could be explained by means of the compressive stresses in the y-direction (described in Fig. 5b) that are located near the restrained edges and which cause the wrinkles apparition. As expected, the critical temperature variation increases with the ratio $r$. In long membranes, the wavelength is larger in the center region than near the boundaries where the maximum compressive stresses are located. Since in the reduced model the wavelength is set to a constant value throughout the membrane, the deformation shape is not well recovered like in the full model. Yet, the mode deformation at the maximum deflection region can still be well recovered (Figs. 23 and 24).

Fig. 25 shows the maximum deflection location according to each value of the wavelength at the bifurcation load. It appears that for an aspect ratio up to 2.75, the maximum deflection is located at the center of the membrane. This same result has been obtained by Zheng in his thesis in 2008 [26] where a dielectric elastomer was subjected to uniform tensile load. Indeed in this uniaxial tensile test, the longitudinal length of the wrinkles is known to be related to the size of the compressive zone that is more or less independent of the plate length [25].

## 5. Conclusion

In this paper, the thermal wrinkling of thin membrane has been analyzed from a multiscale point of view. The method of double scale Fourier series is used to recover a macroscopic membrane model that couples both global and local wrinkling. The membrane parameters studied show a good agreement of the results with the full model, which prove the reliability of the Fourier series approach on 2-D thermal wrinkling problems. The post-bifurcation investigations can also be conducted, however, the results obtained show boundaries effects past the bifurcation point. To deal with this issue, one can use full membrane model discretization at the boundary and reduced model discretization far from the boundary. Both membrane models are then related by coupling techniques such as Arlequin method [27]. This technique has been applied for beams lying on a non-linear elastic foundation [28,29] and sandwich structures [19] and should be extended to membrane structures.

<figure><img src="image_27.png" /><figcaption>**Fig. 25.** Membrane off-center distance of the maximum deflection.</figcaption></figure>
Acknowledgments

The first author is supported by the Fonds National de la
Recherche Luxembourg (FNR) via the AFR Grant PHD-08-069. This
project has been carried out in the framework of the FNR project
WRINKLE (C10/MS/784868). Pr. Heng Hu gratefully acknowledges
the National Natural Science Foundation of China (Grant no.
11372231) and the Science and Technology Agency of the Hubei
Province of China (Grant no. 2011CDA047). F. Mohri and M. Potier-
Ferry acknowledges the support of French National Research
Agency ANR (Labex DAMAS, Grant no. ANR-11-LABX-0008-01).

References

[1] Cuong TD, Tardif N, Limam A. Experimental modeling of flatness defects in strip cold rolling. In: Rolling 2013, 2013, p. 1-11.

[2] Li B, Cao Y-P, Feng X-Q, Gao H. Surface wrinkling of mucosa induced by volumetric growth: theory, simulation and experiment. J Mech Phys Solids 2011;59(April (4)):758-74.

[3] Kim JB, Yoon JW, Yang DY. Wrinkling initiation and growth in modified yoshida buckling test: finite element analysis and experimental comparison. Int J Mech Sci 2000;42:1683-714.

[4] Blandino JR, Jonhston JD, Miles JJ, Soplop JS. Thin film membrane wrinkling due to mechanical and thermal loads. Am Inst Aeronaut Astronaut 2001;1345 (April).

[5] Blandino JR, Johnston JD, Dharamsi UK. Corner wrinkling of a square membrane due to symmetric mechanical loads. J Spacecr Rockets 2002;39(5).

[6] Blandino JR, Jonhston JD, Miles JJ, Dharamsi UK. The effect of asymmetric mechanical and thermal loading on membrane wrinkling. Am Inst Aeronaut Astronaut 2002;1371.

[7] Wong YW, Pellegrino S. Wrinkled membranes part III: numerical simulations. J. Mech. Mater. Struct. 1 January (2006), 2006.

[8] Wang CG, Tan HF, Du XW, Wan ZM. Wrinkling prediction of rectangular shell-membrane under transverse in-plane displacement. Int J Solids Struct October 2007;44(20):6507-16.

[9] Wang CG, Du XW, Tan HF, He XD. A new computational method for wrinkling analysis of gossamer space structures. Int J Solids Struct March 2009; 46(6):1516-26.

[10] Lecieux Y, Bouzidi R. Experimental analysis on membrane wrinkling under biaxial load-comparison with bifurcation analysis. Int J Solids Struct September 2010;47(18-19):2459-75.

[11] Lecieux Y, Bouzidi R. Numerical wrinkling prediction of thin hyperelastic structures by direct energy minimization. Adv Eng Softw 2012;50 (August):57-68.

[12] Healey TJ, Li Q, Cheng R-B. Wrinkling behavior of highly stretched rectangular elastic films via parametric global bifurcation. J Nonlinear Sci April 2013; 23(Aril (5)):777-805.

[13] Roddeman DG, Drukker J, Oomens CWJ, Janssen JD. The wrinkling of thin membranes: part I—theory. J Appl Mech 1987;54(December)).

[14] Roddeman DG, Drukker J, Oomens CWJ, Janssen JD. The wrinkling of thin membranes: part II—numerical analysis. J Appl Mech 1987;54(4).

[15] Chiu HC, Benson RC, Fiscella MD, Burns SJ. Mechanical and thermal wrinkling of polymer membranes. J Appl Mech March 1994;61(1):67-70.

[16] Damil N, Potier-Ferry M. A generalized continuum approach to describe instability pattern formation by a multiple scale analysis. C. R. Méc. 2006;334:674-8.

[17] Damil N, Potier-Ferry M. Influence of local wrinkling on membrane behaviour: a new approach by the technique of slowly variable fourier coefficients. J Mech Phys Solids 2010;58:1139-53.

[18] Liu Y, Yu K, Hu H, Belouettar S, Potier-Ferry M, Damil N. A new fourier-related double scale analysis for instability phenomena in sandwich structures. Int J Solids Struct 2012;49(November (22)):3077-88.

[19] Yu K, Hu H, Chen S, Belouettar S, Potier-Ferry M. Multi-scale techniques to analyze instabilities in sandwich structures. Compos Struct 2013;96 (February):751-62.

[20] Damil N, Potier-Ferry M, Hu H. New nonlinear multi-scale models for wrinkled membranes. C. R. Méc. August 2013;341(8):614-24.

[21] Damil N, Potier-Ferry M, Hu H. Membrane wrinkling revisited from a multi-scale point of view. Adv Model Simul Eng Sci 2014;1(1):6.

[22] Cochelin B, Damil N, Potier-Ferry M. Méthode asymptotique numérique, Lavoisier edition. Paris, Londres: Hermes science; 2007.

[23] Jones RM. Buckling of bars, plates and shells. Blacksburg, VA: Bull Ridge Corporation; 2006.

[24] Friedl N, Rammerstorfer FG, Fischer FD. Buckling of stretched strips. Comput Struct 2000;78:185-90.

[25] Jacques N, Potier-Ferry M. On mode localisation in tensile plate buckling. Comptes Rendus Mécanique 2005;333(October):804-9.

[26] L. Zheng. Wrinkling of dielectric elastomer membranes [PhD thesis]. Passadena, CA: California Institute of Technology; September 2008.

[27] Dhaia HB, Rateau G. The Arlequin method as a flexible engineering design tool. Int J Numer Methods Eng 2005;62(11):1442-62.

[28] Hu H, Damil N, Potier-Ferry M. A bridging technique to analyze the influence of boundary conditions on instability patterns. J Comput Phys 2011; 230(10):3753-64.

[29] Xu F, Hu H, Potier-Ferry M, Belouettar S. Bridging techniques in a multi-scale modeling of pattern formation. Int J Solids Struct 2014;51(18):3119-34.