Prof Hamid Zahrouni

Prof Salim Belouettar

Prof Ahmed Makradi

Dr Mohamed Ei Hachemi

Dr Qian Shao

Dr Yin Liu

Dr Kun Yu

Dr Jie Yang

Siyu Cheng

Yanchuan Hui

Wei Huang

Rui Xu

Jian Liu

Shenyang Liang

Zhengze Cui

Keywords: Membrane, Film/substrate, Sandwich, Instability, Buckling, Wrinkling, Fourier series, Asymptotic Numerical Method Membrane, Film/substrat, Sandwich, Instabilité, Flambage, Rides, Série de Fourier, Méthode numérique asymptotique, Méthode de pontage de domaine biaxial load . . . . . . . . . . . 20 uniaxial load . . . . . . . . . . 22

The main aim of this thesis is to develop advanced and efficient multi-scale modeling and simulation techniques to study instability phenomena in three common engineering structures, i.e., membrane, film/substrate and sandwich structures, by combining the Technique of Slowly Variable Fourier Coefficients (TSVFC) and the Asymptotic Numerical Method (ANM). Towards this end, based on the Von Karman plate equations, the TSVFC has been firstly used to develop a two-dimensional (2D) Fourier double-scale model for membrane, which has also been implemented into ABAQUS via its subroutine UEL. Then a 2D Fourier model is constructed for film/substrate. Further, making use of deformation features of the film/substrate, a 1D Fourier model is developed by using both the TSVFC and the Carreras Unified Formulation (CUF). Subsequently, based on high-order kinematics belonging to Zig-Zag theory, a 2D Fourier model is deduced for sandwich plates. The governing equations for the above models are discretized by the Finite Element Method, and the resulting nonlinear systems are solved by the efficient and robust nonlinear solver ANM. These models are then adopted to study instabilities in these structures. Results show that the established models could accurately and efficiently simulate various instability phenomena. Besides, its found that the membrane instability is very sensitive to boundary conditions, and there exists a dimensionless parameter that is almost constant near bifurcation point for various loading cases and geometric parameters, which may be helpful for fast predicting the occurrence of wrinkles.

Résumé

L'objectif de cette thèse est de développer des techniques de modélisation et de simulation multi-échelle avancées et efficaces pour étudier les phénomènes d'instabilité dans trois structures d'ingénierie courantes: membrane, film/substrat et structures sandwich, en combinant la technique des coefficients de Fourier lentement variables (TSVFC) et la méthode numérique asymptotique (ANM). À cette fin, basée sur les équations de la plaque de Von Karman, la TSVFC été utilisée pour développer un modèle de Fourier à bidimensionnel (2D) qui a également été implémenté dans ABAQUS via sa sous-routine UEL. Ensuite, un 2D modèle de Fourier est construit pour le film/substrat. En outre, en utilisant leurs caractéristiques de déformation, un 1D modèle de Fourier est développé en utilisant à la fois le TSVFC et le CUF.

Par la suite, sur la base d'une cinématique Zig-Zag d'ordre supérieur, un 2D modèle de Fourier est déduit pour une plaque sandwich. Les équations directrices pour les modèles ci-dessus sont discrétisées par la méthode des éléments finis, et les systèmes non linéaires résultants sont résolus par le solveur non linéaire efficace et robuste ANM. Ces modèles sont ensuite adoptés pour étudier les instabilités dans ces structures. Les résultats montrent que les modèles établis peuvent simuler avec précision et efficacité divers phénomènes d'instabilité. En outre, on constate que l'instabilité membranaire est sensible aux conditions aux limites et qu'il existe un paramètre sans dimension presque constant près du point de bifurcation pour différents cas de charge et paramètres géométriques, ce qui peut être utile pour prédire rapidement l'apparition des rides.

Chapter 1

Introduction

Instability phenomena are widely observed in engineering (such as aerospace, automotive, smart materials, biomedical, etc.). The occurrence of wrinkles may pose a limit on the performance of materials or structures and is often thought to be avoided. For example, in the rolling process of thin metal sheet [START_REF] Jacques | Buckling and wrinkling during strip conveying in processing lines[END_REF]; [START_REF] Damil | Influence of local wrinkling on membrane behaviour: a new approach by the technique of slowly variable Fourier coefficients[END_REF]) (as shown in Fig. (1.1)), the compressive stress induced by the strong tensile stress causes the central fold near the bite, while the residual stress produced by the non-uniform plastic strain causes wrinkles in the edge of the thin metal sheet. If the amplitude of the wrinkles is too large, tens of meters of the sheet should be cut and thrown out, which may lead to at least one day off for the production line. When stretching the metal film (see Fig.

(1.2)), compressive stresses perpendicular to the stretching direction, caused by the Poisson effect, may also lead to local wrinkles in the film. Similar instability phenomena can also be found in the following structures: the overall buckling and local wrinkling of the sandwich structures (Hu et al. (2009a); [START_REF] Liu | A new Fourierrelated double scale analysis for instability phenomena in sandwich structures[END_REF]), long-fibre microbuckling phenomena of composite structures [START_REF] Drapier | A structural approach of plastic microbuckling in long fibre composites: comparison with theoretical and experimental results[END_REF]), etc.. In recent years, the wrinkles with the periodic nature have found some applications such as functional morphology [START_REF] Mei | Buckling modes of elastic thin films on elastic substrates[END_REF]), the measurement of mechanical properties of materials [START_REF] Stafford | A buckling-based metrology for measuring the elastic moduli of polymeric thin films[END_REF]; [START_REF] Chung | Surface wrinkling: A versatile platform for measuring thin-film properties[END_REF]) and flexible electronic devices design [START_REF] Kim | Stretchable Electronics: Materials Strategies and Devices[END_REF]). For example, the film/substrate system can be obtained by heating the soft substrate, then depositing the nano-film, and cooling the system (Bowden et al. (1998), see Fig. (1.3)). Due to the mismatch in thermal expansion coefficient between the soft substrate and the nano-film, the residual stress produced during the cooling process will cause the micro-bending of the nano-film. By measuring the wavelength of the wrinkles, the elastic modulus of the nano-film can be predicted easily [START_REF] Stafford | A buckling-based metrology for measuring the elastic moduli of polymeric thin films[END_REF]).

No matter hundreds of nanometers for graphene sheet or tens of meters for metal sheet, the wrinkling phenomena have two common characters: 1) the instability pattern is spatially nearly periodic, and 2) often the modal wavelength is very small compared to the membrane size.

Compared to experiments (Wong and Pellegrino (2006a); [START_REF] Wang | A new computational method for wrinkling analysis of gossamer space structures[END_REF]; [START_REF] Lecieux | Experimental analysis on membrane wrinkling under biaxial load-Comparison with bifurcation analysis[END_REF]; [START_REF] Nayyar | Stretch-induced wrinkling of polyethylene thin sheets: Experiments and modeling[END_REF]), numerical simulation may be [START_REF] Bowden | Spontaneous formation of ordered structures in thin filmsofmetals supported on an elastomeric polymer[END_REF] the best option to predict these instabilities in the view of time and economic cost. For example, to find out and avoid the main causes that lead to wrinkles during the rolling process of thin metal, the nonlinear shell theory (see [START_REF] Wang | A new computational method for wrinkling analysis of gossamer space structures[END_REF]; Wong and Pellegrino (2006b); [START_REF] Lecieux | Numerical wrinkling prediction of thin hyperelastic structures by direct energy minimization[END_REF]; [START_REF] Healey | Wrinkling Behavior of Highly Stretched Rectangular Elastic Films via Parametric Global Bifurcation[END_REF]) could be used, which could accurately describe the details of the thin metal response: size, wavelength, orientation of the wrinkles, instability threshold, etc..

Then by parameter study and analyses that nowadays could be realized with the aid of powerful and stable commercial software, it could be more easier, faster and cheaper than experiments. However, the traditional numerical strategies have two difficulties: 1) it is difficult to take into account both the computational efficiency and accuracy, and 2) it is difficult to control the nonlinear calculation in cases with a large number of wrinkles and therefore with a large number of equilibrium solutions. Thus, there is a need to develop advanced multi-scale modelling and simulation techniques to accurately predict the wrinkling phenomena in a computationally efficient way. The developed methodology will be extended to the following structures that are widely used in the engineering: (i) wrinkling of thin metal membranes; (ii) buckling and postbuckling phenomena of nano-film/substrate systems; (iii) buckling and wrinkling of sandwiches.

The instability phenomena can be modelized by bifurcation analysis according to the Landau-Ginzburg theory [START_REF] Wesfreid | Cellular Structures in instabilities[END_REF]). This famous Landau-Ginzburg equation follows from an asymptotic double scale analysis. At the local level one accounts for the periodic nature of the wrinkles, while the slow variations of the envelope are described at the macroscopic scale. Few elements are required to simulate the slow variation of the envelope at the macroscopical scale. However, it has several drawbacks: 1) the solution is not valid far away from the bifurcation, 2) it could not capture the coupling between a global buckling and local wrinkling and 3) the boundary effects are neglected.

A slightly different approach has been proposed recently, where the nearly periodic fields are represented by Fourier series with slowly varying coefficients. This leads to macroscopic models, that are generalized continua, the macroscopic field being defined by Fourier coefficients of the microscopic field (Damil andPotier-Ferry (2006, 2008)). In this sense, the technique is similar to homogenization theory, where a double scale analysis permits to deduce macroscopic generalized continua from microscopic classical ones [START_REF] Forest | Cosserat overall modeling of heterogeneous materials[END_REF]; [START_REF] Kouznetsova | Multi-scale second-order computational homogenization of multi-phase materials: a nested finite element solution strategy[END_REF]).

For example, in order to describe the response of the system which is the sum of a slowly varying mean field and a fluctuation that is nearly periodic in one spatial direction, at least two slowly varying functions (mean and amplitude) are needed to model the phenomenon. These above double scale methods are generic and could be applied in all cases of instabilities with spatially periodic pattern. For instance, the asymptotic approach has been used for Rayleigh-Benard convection [START_REF] Newell | Finite band width, finite amplitude convection[END_REF]; [START_REF] Segel | Distant side walls cause slow amplitude modulation of cellular convection[END_REF]) and for many fluid instabilities [START_REF] Wesfreid | Cellular Structures in instabilities[END_REF]; [START_REF] Cross | Pattern formation out of equilibrium[END_REF]), for the buckling of a beam on foundation [START_REF] Amazigo | Asymptotic analyses of the buckling of imperfect columns on non-linear elastic foundations[END_REF]; [START_REF] Pomeau | Wavelength selection in one-dimensional cellular structures[END_REF]; [START_REF] Potier-Ferry | Amplitude modulation, phase modulation and localization of buckling patterns[END_REF]), for plate buckling [START_REF] Damil | Wavelength selection in the postbuckling of a long rectangular plate[END_REF]) and cylindrical shell buckling [START_REF] Amazigo | Buckling under external pressure of cylindrical shells with dimple shaped initial imperfection[END_REF]; [START_REF] Abdelmoula | Influence of distributed and localized imperfections on the buckling of cylindrical shells[END_REF]). In the same way, such cellular instabilities appear in the buckling of thin elastic film bound to compliant substrate [START_REF] Chen | Herringbone buckling patterns of compressed thin films on compliant substrates[END_REF]; [START_REF] Wang | Local versus global buckling of thin films on elastomeric substrates[END_REF]; [START_REF] Audoly | Buckling of a stiff film bound to a compliant substrate-Part II: A global scenario for the formation of herringbone pattern[END_REF] which will be studied in this thesis.

The second approach, based on the Fourier series, is also able to account for the coupling between local and global buckling, as in the buckling of stiffened plate [START_REF] Sridharan | Stiffened plates and cylindrical shells under interactive buckling[END_REF]) or of sandwich structures (Hu et al. (2009a); Léotoing et al. (2002a)) or to account for the influence of wrinkles on the behaviour of membranes [START_REF] Wang | A new computational method for wrinkling analysis of gossamer space structures[END_REF]; [START_REF] Nayyar | Stretch-induced wrinkling of polyethylene thin sheets: Experiments and modeling[END_REF]; Wong and Pellegrino (2006c)).

In this thesis, the technique of slowly varying Fourier coefficients will be applied to obtain the so called "macroscopic model" or "reduced model", with the aim to represent the influence of local buckling on the thin walled structures, e.g., thin metal sheet, film/substrate systems and sandwich structures behaviour in a computationally efficient way. Moreover, the macroscopic model allows one to select the wavelength of wrinkles that makes it easier to pilot the nonlinear simulation.

The established nonlinear systems in this thesis will be solved by the Asymptotic Numerical Method (ANM), see [START_REF] Damil | A new method to compute perturbed bifurcations: Application to the buckling of imperfect elastic structures[END_REF]. The ANM is an alternative method which falls into the category of numerical perturbation techniques. By introducing the power series expansions into the equilibrium equation, the nonlinear problem is transformed into a sequence of linear problems and solved by the standard finite element method. Many studies [START_REF] Zahrouni | Computing finite rotations of shells by an asymptotic-numerical method[END_REF]; [START_REF] Cadou | ANM for stationary Navier-Stokes equations ans with Petrov-Galerkin formulation[END_REF]; [START_REF] Cochelin | Méthode asymptotique numérique[END_REF]) presented in the literature show that the ANM is more reliable and less time consuming than the iterative classical methods, for instance the modified Newton method and the Newton-Raphson method. This thesis will be structured as follows, each chapter corresponding to a published research work with a self-consistent notation:

• in Chapter 2: [Int J Solids Struct, 64-65: 246-258 (2015)], the technique of slowly varying Fourier coefficient is introduced in detail to deduce a macroscopic model for membrane. The model is first implemented in a home-made code, the nonlinear system being solved by the Asymptotic Numerical Method (ANM), which has advantages of efficiency and reliability for stability analyses. It is also implemented as a user element in a commercial software to evaluate the effectiveness of the reduction technique. Various loading cases are considered and the numerical tests show that the macroscopic model can accurately and fast predict the wrinkling patterns. The numerical results highlight the strong influence of a dimensionless parameter for wrinkling • in Chapter 3: [Thin Wall Struct, 114: 144-153 (2017)], the technique of slowly varying Fourier coefficient is extended for the study of local wrinkling phenomena in thin film/substrate system. The established model, in which the film is modelled by Euler-Bernoulli beam and the substrate by continuum solid model, is demonstrated to be efficient and accurate in predicting the critical load and wrinkling pattern with several numerical examples. Besides, the importance of using the first harmonic of Fourier series in approximating the axial displacement in substrate is, in particular, discussed, and a spurious stiffening effect is pointed out. To overcome this phenomenon, modifications on either the Fourier series or the constitutive equation of the substrate are proposed.

• in Chapter 4: [Compos Struct, 160: 613-624 (2017)], a novel one-dimensional Fourier finite element is proposed for the thin film/substrate structure, in which the substrate is modelled by the Carrera's Unified Formulation (CUF) and the film still by the Euler-Bernoulli beam as in Chapter 3. The numerical results show that the new Fourier model yields accurate results with lower computational cost. Besides, several harmonics of Fourier series are taken into account in this new model, which proves to efficiently trace the wrinkling pattern corresponding to the lowest critical load. The established nonlinear system is solved by the Asymptotic Numerical Method (ANM).

• in Chapter 5: [Comput Method Appl M, 318: 270-295 (2017)], based on the layer-wise kinematic proposed by [START_REF] Yu | A novel two-dimensional finite element to study the instability phenomena of sandwich plates[END_REF], a Fourier-based double scale model is developed for the local wrinkling of sandwich plate. Both antisymmetrical and symmetrical wrinkling for sandwich plates under different kinds of loads are investigated and the obtained results demonstrate that the Fourier-based finite element model can accurately yet efficiently predict wrinkling patterns and critical loads, especially when dealing with wrinkling phenomena with extremely large wavenumbers.

Finally, some conclusions are made by summarizing main outcomes of the proposed theoretical and numerical approaches, and possible future works are also given.

the wrinkling patterns, even when there are few wrinkles. The numerical results highlight the strong influence of a dimensionless parameter for wrinkling initiation.

Present Chapter corresponds to the published research paper [START_REF] Huang | Macroscopic simulation of membrane wrinkling for various loading cases[END_REF] [International Journal of Solids and Structures, 64-65: 246-258, 2015.]. A selfconsistent notation is adopted.

Keywords: Fourier series, Membrane, Instability, Wrinkling, Asymptotic Numerical Method

Introduction

Membranes have recently seen a surge of interest in spacecraft structures (sunshields, solar sails, solar panels, reflector antennas and inflatable membranes), civil structures (pneumatic structures for roofs), biomedical materials. Due to the almost negligible bending stiffness and the associated inability to withstand compressive loading, wrinkles arise in very thin membranes, which has an adverse effect on the static and dynamic characteristics and longevity of membrane structures. Hence, to predict and avoid wrinkles, modeling and numerical analysis of membrane wrinkling has been a subject of interest.

The pioneer work on membrane wrinkling belongs to [START_REF] Wagner | Flat sheet metal girders with very thin metal web[END_REF] who developed a tension field theory, which assumes that out-of-plane flexural and in-plane compressive stress are negligible and bending stiffness is neglected. On the basis of tension field theory, many approaches have been put forward to address issues in membrane wrinkling. In most of them, see Stein and Hedgepeth (1961); [START_REF] Liu | Large deflection analysis of pneumatic envelopes using a penalty parameter modified material model[END_REF]; [START_REF] Rossi | Simulation of light-weight membrane structures by wrinkling model[END_REF]; [START_REF] Jarasjarungkiat | A wrinkling model based on material modification for isotropic and orthotropic membranes[END_REF], the stress-strain relationship is modified to eliminate compressive stresses and a wrinkling criterion is established by distinguishing three membrane states: taut, wrinkled and slack.

Another method by Roddeman et al. (1987a,b) splits the deformation tensor into two parts, the membrane part and wrinkling part. It has also attracted many researchers, see for instance [START_REF] Miyazaki | Wrinkle/slack model and finite element dynamics of membrane[END_REF]; [START_REF] Akita | A simple computer implementation of membrane wrinkle behaviour via a projection technique[END_REF]; [START_REF] Shaw | Analyses of wrinkled and slack membranes through an error reproducing mesh-free method[END_REF]; [START_REF] Pimprikar | New computational approaches for wrinkled and slack membranes[END_REF]. As the partial differential equations deduced from these membrane models are not elliptic (or hyperbolic in the dynamical case), the presence of slack regions may result in near singular stiffness matrices leading to difficulties in numerical solution. This problem can become well posed if an internal length is included, for instance within Cosserat theory [START_REF] Pimprikar | New computational approaches for wrinkled and slack membranes[END_REF][START_REF] Banerjee | The theory of Cosserat points applied to the analyses of wrinkled and slack membranes[END_REF], but perhaps this regularization is not necessary when us-ing an explicit dynamic computation. However, all these membrane models just characterize the stress distribution and the wrinkled region, but cannot identify the details of wrinkles such as the amplitude of the wrinkles, their wavelength, the sensitivity to boundary conditions and the instability critical load.

In general, shell elements accounting for very small bending stiffness are preferred to qualitatively characterize the details of membrane wrinkles. In recent papers, Wong and Pellegrino (2006b) presented a general procedure for simulating the onset and growth of wrinkles with shell element using the commercial finite element package ABAQUS and obtained consistent results with experiments in Wong and Pellegrino (2006a,c). [START_REF] Wang | A new computational method for wrinkling analysis of gossamer space structures[END_REF] proposed a new Modified Displacement Component (MDC) method to eliminate the singularity of stiffness matrix and apply it to a finite element code ANSYS. [START_REF] Flores | Wrinkling and folding analysis of elastic membranes using an enhanced rotation-free thin shell triangular element[END_REF] used a rotation-free linear strain shell triangle element within explicit time integration strategy that avoids introducing geometric imperfection and initial stress. [START_REF] Lecieux | Numerical wrinkling prediction of thin hyperelastic structures by direct energy minimization[END_REF] adopted the same shell element as Flores to analyze the wrinkling of membranes by direct energy minimization. In all these models, a large number of degrees of freedom are needed to capture the details of wrinkles that results in heavy computational cost and convergence is difficult to achieve. This kind of shell models will be referred as full shell model.

In this Chapter, we present a new finite element based on a reduced-order model established recently from a multiple scale analysis [START_REF] Damil | New nonlinear multi-scale models for wrinkled membranes[END_REF][START_REF] Damil | Membrane wrinkling revisited from a multiscale point of view[END_REF]. It couples a non-linear 2D membrane model with an envelope equation governing the evolution of wrinkles. This membrane theory is an extension of the Landau-Ginzburg bifurcation equation and it has been deduced by using the method of Fourier series with slowly variable coefficients [START_REF] Damil | Influence of local wrinkling on membrane behaviour: a new approach by the technique of slowly variable Fourier coefficients[END_REF]. The origin of the reduction lies in the double scale approach so that the needed spatial meshes can be considered as macroscopic and are not related with the wrinkling wavelength.

The finite element technology is described in details and this has been implemented in a home-made code and also in the commercial package ABAQUS as a user element. In the first code, the governing equations are solved by the Asymptotic Numerical Method that is an efficient and robust path following technique in the presence of bifurcations. Several numerical simulations will be discussed in order to evaluate the range of validity of the reduced model, to establish the efficiency of the reduction and to analyze the case of a uniaxially stretched plate [START_REF] Friedl | Buckling of stretched strips[END_REF][START_REF] Jacques | On mode localisation in tensile plate buckling[END_REF] that is very sensitive to boundary conditions. This Chapter is structured as follows. Section 2.2 recalls the famous Föppl-Von Karman equations for isotropic plates that will be considered as the reference model. In Section 2.3, the finite element procedure will be described. It includes the reduced-order model of [START_REF] Damil | New nonlinear multi-scale models for wrinkled membranes[END_REF][START_REF] Damil | Membrane wrinkling revisited from a multiscale point of view[END_REF], a 2D finite element discretization and the resolution technique by ANM. In Section 2.4, two numerical examples are investigated to evaluate the range of validity and to bring out the strong influence of boundary conditions on wrinkling. In Section 2.5, a dimensionless parameter K is introduced that is important to predict the appearance of membrane instability. Finally, in Section 2.6, we describe the implementation into ABAQUS as user element (UEL), especially to obtain a reliable evaluation in terms of computation time.

Starting plate model

In this section, the well known Föppl-Von Karman equations for elastic isotropic plates will be considered as the reference model. Sometimes we also use a commercial code, where the shell finite element model remains valid for large rotations, but the difference between these two models is weak in the cases considered in this Chapter. The Föppl-Von Karman equations can be written as:

                           div N = 0, (a)
N = L m • γ, (b) 2γ = ∇u + t ∇u + ∇w ⊗ ∇w, (c) D∆ 2 w -div(N • ∇w) = 0. (d) (2.1)
where u = (u, v) is the in-plane displacement and w is the out-plane displacement, 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

[L m ] = Eh 1 -ν 2     1 ν 0 ν 1 0 0 0 1 -ν 2     . (2.2)
The corresponding internal work W int can be split into a membrane part W m and a bending part W b as follows:

                           W int (u, w) = W b (w) + W m (u, w), 2W b (w) = D ∫ ∫ ω ( (∆w) 2 -2(1 -ν) ( ∂ 2 w ∂x 2 ∂ 2 w ∂y 2 - ( ∂ 2 w ∂x∂y ) 2 )) dω, 2W m (u, w) = ∫ ∫ ω ⟨γ⟩[L m ]{γ}dω = Eh 1 -ν 2 ∫ ∫ ω (γ 2 x + γ 2 y + 2(1 -ν)γ 2 xy + 2νγ x γ y ) dω.
(2.3)

Macroscopic membrane model

In this part, a macroscopic finite element model will be presented that is based on a reduced-order modelling introduced in [START_REF] Damil | New nonlinear multi-scale models for wrinkled membranes[END_REF][START_REF] Damil | Membrane wrinkling revisited from a multiscale point of view[END_REF]. The mathematical model has been deduced from the full shell model by the method of Fourier series with slowly variable coefficients [START_REF] Damil | Influence of local wrinkling on membrane behaviour: a new approach by the technique of slowly variable Fourier coefficients[END_REF]. The principle of this Fourier-related approach is to write the unknown field in the following form:

U (x, y) = +∞ ∑ m=-∞ U m (x, y)e miQx .
(2.4)

where the wavenumber Q is given and the macroscopic unknown field U m (x, y)

slowly varies on a period

[ x, x + 2π Q ]
of the oscillation. In the case of Föppl-Von Karman equations (2.1), the unknown field is U (x, y) = (u(x, y), w(x, y), N(x, y), γ(x, y)).

Of course, in practice only a finite number of Fourier coefficients will be considered.

As pictured in Fig.

(2.1), at least two functions U 0 (x, y) and U 1 (x, y) are necessary to describe nearly periodic patterns: U 0 (x, y) can be identified to the mean value while U 1 (x, y) represents the envelope or the amplitude of the spatial oscillations.

In this work, the unknown fields U (x, y) are expressed only in terms of two harmonics: the mean field U 0 (x, y) and the first order harmonics U 1 (x, y)e iQx and U 1 (x, y)e -iQx . The mean value U 0 (x, y) is real-valued, while the other envelops can be complex-valued. So the envelope of the first harmonic U 1 (x, y) can be written as U 1 (x, y) = r(x, y)e iφ(x) , where r(x, y) represents the amplitude modulation and φ(x) the phase modulation. If the phase linearly varies (φ(x) = qx + φ 0 ), this type of approach is able to describe quasi-periodic responses whose wavenumber Q + q slightly differs from the a priori chosen Q. Hence, the method could account for a change of wavenumber.

Here we want to use the simplest consistent model that is able to couple a membrane stress with the variable amplitude of wrinkling pattern. Hence we limit ourselves to real values of the envelope U 1 (x, y). As shown in [START_REF] Mhada | About macroscopic models of instability pattern formation[END_REF], this phase locking can lead to perturbations near the boundary, but an alternative reduced model would be much more expensive. Moreover as in [START_REF] Damil | Influence of local wrinkling on membrane behaviour: a new approach by the technique of slowly variable Fourier coefficients[END_REF] and [START_REF] Damil | Membrane wrinkling revisited from a multiscale point of view[END_REF], the displacement field is reduced to a wrinkling amplitude, i.e. u 1 = 0, w 0 = 0. Therefore the macroscopic displacement is reduced to three components u 0 = (u 0 , v 0 ) and w 1 that can be rewritten for simplicity as

(u, v) def = (u 0 , v 0 ), w def = w 1 .
It would be not difficult to keep the mean deflection w 0 as in [START_REF] Liu | A new Fourierrelated double scale analysis for instability phenomena in sandwich structures[END_REF], what will be necessary to couple local and global buckling.

Membrane strain energy

The derivation rules in [START_REF] Damil | Influence of local wrinkling on membrane behaviour: a new approach by the technique of slowly variable Fourier coefficients[END_REF] have been extended in this bi-dimensional framework, see [START_REF] Damil | Membrane wrinkling revisited from a multiscale point of view[END_REF]. This leads to the following useful formulae (the indices 1 and 0 refer to the number m of the harmonics as in Eq. (2.4)):

{(∇w) 1 } =            ∂w 1 ∂x + iQw 1 ∂w 1 ∂y            =            ∂w ∂x + iQw ∂w ∂y            , (2.5)
and

{γ} def = {γ 0 } =                          ∂u ∂x + ( ∂w ∂x ) 2 + Q 2 w 2 ∂v ∂y + ( ∂w ∂y ) 2 ∂u ∂y + ∂v ∂x + 2 ∂w ∂x ∂w ∂y                          . (2.6)
The membrane strain formula Eq. (2.6) is quite similar to the strain of the initial Von Karman model. It can be split, first into a linear part ε(u) that is the symmetric part of the displacement gradient and corresponds to the pure membrane strain, second into a nonlinear part γ nl (w) that is more or less equivalent to wrinkling deformation in Roddeman et al. (1987a) and is given by

{γ wr (w)} =                          ( ∂w ∂x ) 2 + Q 2 w 2 ( ∂w ∂y ) 2 2 ∂w ∂x ∂w ∂y                          . (2.7)
The main difference with the classical Von Karman strain (2.1-c) is the extension Q 2 w 2 in the direction of the wrinkles. If the linear strain is compressive, as wrinkling strain Q 2 w 2 is always positive, the wrinkling leads to a decrease of the membrane strain.

One may wonder why derivatives are kept in Eqs. (2.5),(2.6),(2.7). Indeed, the assumption of slowly varying envelopes leads to ∂/∂x ≪ Q, ∂/∂y ≪ Q. This would lead to pure membrane theory with a nonlinear relation between deflection w and membrane stress N, but it has been established [START_REF] Damil | Membrane wrinkling revisited from a multiscale point of view[END_REF] that this could not permit to define the wrinkling wavelength and the wrinkling stress.

By substituting Eq. (2.6) into the third formula of Eq. (2.3), the membrane strain energy is obtained:

2W m (u, w) = Eh 1 -ν 2 ∫ ∫ ω                 ( ∂u ∂x + ( ∂w ∂x ) 2 + Q 2 w 2 ) 2 + ( ∂v ∂y + ( ∂w ∂y ) 2 ) 2 +2(1 -ν) ( 1 2
( ∂u ∂y + ∂v ∂x ) + ∂w ∂x ∂w ∂y ) 2 +2ν ( ∂u ∂x + ( ∂w ∂x ) 2 + Q 2 w 2 ) ( ∂v ∂y + ( ∂w ∂y ) 2 )                 dω.
(2.8)

The variation of the membrane energy (2.8) can be written in a concise form:

δW m = ∫ ∫ ω ⟨δγ⟩ {N} dω, (2.9) in which, {N} def = {N 0 } = [L m ] {γ} .
(2.10)

Bending energy

The bending energy is calculated in the same framework: u 1 = (u 1 , v 1 ) = (0, 0), w 0 =0 and w 1 is real. The computation of the energy can be simplified by keeping only the 0th order in bending energy. Here the derivation process of bending energy is not further described here, more details could be seen in [START_REF] Damil | Membrane wrinkling revisited from a multiscale point of view[END_REF]. So the bending energy could be expressed as:

W b (w) = D ∫ ∫ ω { Q 4 w 2 -2Q 2 w∆w + 4Q 2 ( ∂w ∂x ) 2 + 2(1 -ν)Q 2 [ w ∂ 2 w ∂y 2 + ( ∂w ∂y ) 2 ]} dω.
(2.11)

Hence,

δW b = ∫ ∫ ω ⟨δβ⟩{M} dω, (2.12) in which {β} = t ⟨w ∂w ∂x ∂w ∂y ⟩, (2.13) {M} = [L f ]{β}, [ L f ] = 2D     Q 4 0 0 0 6Q 2 0 0 0 2Q 2     . (2.14)
Partial differential equations are easily deduced from the reduced elastic energies

(2.8),(2.11) and few analytical solutions are available [START_REF] Damil | Membrane wrinkling revisited from a multiscale point of view[END_REF]. In this Chapter, we focus on the numerical aspects.

Weak form of the macroscopic membrane model

The macroscopic membrane model is deduced from the principle of virtual work

δW int = δW ext , (2.15)
where δW int and δW ext are the internal and external virtual work. In this Chapter, as the body forces are neglected, δW ext could be expressed in the following form: (2.16) where δu and δw are virtual in-plane and out-plane displacements, {F m } and F b denote respectively the in-plane and out-plane external forces, and λ is a "load parameter". Hence, Eq. (2.15) is transformed into

δW ext = λ ∫ Sσ {⟨δu⟩{F m } + δwF b } dS,
∫ ∫ ω { ⟨δγ⟩[L m ]{γ} + ⟨δβ⟩[L f ]{β} } dω = λ ∫ Sσ {⟨δu⟩{F m } + δwF b } dS. (2.17)
The unknowns in Eq. (2.17) can be noted by the following vector: (2.18)

{θ} = t ⟨
By using this definition, the strain vector {γ} could be written in the following form:

{γ} = ( [H] + 1 2 [A(θ)] ) {θ} , (2.19)
where

[H] =     1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 1 0 0 0 0     , (2.20) [A(θ)] = 2     0 0 0 0 Q 2 w ∂w ∂x 0 0 0 0 0 0 0 ∂w ∂y 0 0 0 0 0 ∂w ∂y ∂w ∂x     .
(2.21) Hence, the stress vector {N} reads:

{N} = [L m ] ( [H] + 1 2 [A(θ)] ) {θ} . (2.22)
To unify the unknowns, we define a transform matrix [T β ] that establishes the relation between {θ} and {β}:

{β} = [T β ]{θ}, (2.23)
where

[T β ] =     0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1     . (2.24)
Finally, the governing equation of the macroscopic membrane model Eq. (2.17) is in the following form:

∫ ∫ ω ⟨δθ⟩ { ( t [H] + t [A(θ)] ) {N} + t [T β ][L f ][T β ] {θ} } dω = λ ∫ Sσ {⟨δu⟩{F m } + δwF b } dS.
(2.25)

In this Section, this macroscopic membrane model would be referenced as "new reduced model".

Discretization

The 2D-Q8 plate finite element, defined by eight nodes and three degrees of freedom per node (u, v and w), is used. The nodal displacement of this 2D element is described by using the classical serendipity shape functions: (2.26) where {q} e are the nodal unknowns per element and N i are the shape functions.

                 u = ⟨N u ⟩{q} e = ⟨N 1 0 0 • • • N 8 0 0⟩{q} e , v = ⟨N v ⟩{q} e = ⟨0 N 1 0 • • • 0 N 8 0⟩{q} e , w = ⟨N w ⟩{q} e = ⟨0 0 N 1 • • • 0 0 N 8 ⟩{q} e ,
Therefore the elementary unknown vector {θ} e is in the following form: (2.27) where

{θ} e = [G]{q} e ,
[G] =               N u,x 0 0 N u,y 0 0 0 N v,x 0 0 N v,y 0 0 0 N w 0 0 N w,x 0 0 N w,y               .
(2.28)

Using this discretization, the elementary virtual internal work could be expressed in the following semi-discretized form:

δW e int = ⟨δq⟩ e ∫ ∫ ωe t [G] ( t [H] + t [A(θ(q))] ) {N} + t [G] t [T β ][L f ][T β ][G] {q} e dω.
(2.29)

We define

[B(q)] = [H + A(θ(q))] [G],
(2.30) (2.33)

[B l ] = [H] [G], [B nl (q)] = [A(θ(q))] [G]. ( 2
where [N] is the shape function matrix, {F} is the applied external force.

Finally, the discretization form of the governing equation for macroscopic model is expressed by (2.34) and

∑ e ⟨δq⟩ e ∫ ∫ ωe { t [B(q)] {N} + t [G] t [T β ][L f ][T β ][G] {q} e } dω = λ ∑ e ⟨δq⟩ e {f } e ,
{N} = [L m ] ( [B l ] + 1 2 [B nl (q)]
) {q} e .

(2.35)

Implementation of the ANM

The above non-linear equations are solved by Asymptotic Numerical Method (ANM) (see [START_REF] Cochelin | Méthode asymptotique numérique[END_REF][START_REF] Cochelin | Asymptotic-numerical methods and Padé approximants for non-linear elastic structures[END_REF]). The solution at the end of j th step, (q j , N j , λ j ) is assumed to be known so as to determine the solution at the end of (j + 1) th step, (q j+1 , N j+1 , λ j+1 ). Following the perturbation technique, an approached solution path, in the step j + 1, is represented by truncated power series at order N order : (2.36)

N j+1 = N j + N order ∑ k=1 a k N k = N j + aN 1 + a 2 N 2 + . . . ,
q j+1 = q j + N order ∑ k=1
a k q k = q j + aq 1 + a 2 q 2 + . . . , (2.37)

λ j+1 = λ j + N order ∑ k=1 a k λ k = λ j + aλ 1 + a 2 λ 2 + . . . , (2.38)
where a is a path parameter defined as

a = ⟨Q 1 ⟩({Q j+1 } -{Q j }) + λ 1 (λ j+1 -λ j ).
(2.39)

Substituting the approached solution Eqs. (2.36)-(2.39) in the equations Eqs.

(2.34), (2.35), results in a linear system of equations at each order k. The algebraic equations for the order 1 and the order k are described as following:

Order 1                  ∑ e ⟨δq⟩ e ∫∫ ωe ( t [B(q j )] ) {N 1 } + t [B nl (q 1 )]{N j } + t [G] t [T β ][L f ][T β ][G] {q 1 } e dω = ∑ e ⟨δq⟩ e λ 1 {f } e , {N 1 } = [L m ] [ B(q j ) ] {q 1 } e , ⟨Q 1 ⟩ {Q 1 } + λ 1 2 = 1.
(2.40)

The third formula in Eq. (2.40) represents a complementary condition resulting from the definition of the path parameter. The term t [B nl (q 1 )]{N j } can be expressed as: (2.41) where

t [B nl (q 1 )]{N j } = t [G][ Nj ][G]{q 1 } e ,
[ Nj ] = 2               0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Q 2 N j x 0 0 0 0 0 0 0 N j x N j xy 0 0 0 0 0 N j xy N j y               . (2.42)
The use of the matrix [ N] permits to transform the internal work into a bilinear form, and this makes it easy to apply the ANM.

Substituting Eq. (2.41) into Eq. (2.40), we get: (2.43) where

∑ e ⟨δq⟩ e [ k(q j ) ] {q 1 } e = λ 1 ∑ e ⟨δq⟩ e {f } e ,
[ k(q j ) ]
is the elementary stiffness matrix given below:

[ k(q j ) ] = ∫ ∫ ωe t [B(q j )][L m ][B(q j )] + t [G][ Nj ][G] + t [G] t [T β ][L f ][T β ][G] dω.
(2.44)

Suppose that the vectors {Q}, {F} are the assembly of {q} e , {f } e . After assembling the elementary stiffness matrix, a linear equation for the global system at order 1 is given below:

[K(Q j )]{Q 1 } = λ 1 {F}. (2.45) Order k                  ∑ e ⟨δq⟩ e ∫∫ ωe ( t [B(q j )] ) {N k } + t [B nl (q k )]{N j } + t [G] t [T β ][L f ][T β ][G] {q} e dω = ∑ e ⟨δq⟩ e (λ k {f } e + {f nl intk } e ), {N k } = [L m ] [ B(q j ) ] {q k } e + {N nl k }, ⟨Q 1 ⟩ {Q k } + λ 1 λ k = 0.
(2.46)

In contrast to first order, there are two right-hand sides at each order k > 1 caused by nonlinear terms: the first one, {f nl intk } e and the second one, {N nl k }. They depend on the displacements and stresses at orders anterior to the order k: where

{f nl intk } e = - k-1 ∑ m=1 ∫ ∫ ωe t [B nl (q k-m )]{N m } dω, (2.47) {N nl k } = [L m ] k-1 ∑ m=1 1 2 [B nl (q k-m )]{q m } e . ( 2
{f nl k } e = {f nl intk } e - ∫ ∫ ωe t [ B(q j ) ] {N nl k } dω.
(2.50)

After assembling the elementary stiffness matrix, the linear equation for each order k is:

[K(Q j )]{Q k } = λ k {F} + {F nl k }.
(2.51)

Here {F nl k } is the assembly of {f nl k } e . Eqs. (2.45) and (2.51) mean that all orders in one step have the same stiffness matrix, that is to say, only one stiffness matrix will be inversed in one step. After the values of {Q k } are calculated, then the solutions at step (j + 1)th can be obtained by:

{Q j+1 } = {Q j } + a{Q 1 } + a 2 {Q 2 } + a 3 {Q 3 } + . . . (2.52)
The maximum value of path parameter a is adopted with the same principle as in [START_REF] Cochelin | Méthode asymptotique numérique[END_REF],

a = ( ϵ ∥Q 1 ∥ ∥Q k ∥ ) 1 N order -1 , (2.53)
where ϵ is a precision parameter, which is determined by the user.

Numerical results and discussions

The new reduced model will be investigated through two numerical examples: 1) rectangular membrane under biaxial compression-tension load, 2) rectangular membrane under uniaxial tension load. In the first example, there are many wrinkles and this will be an opportunity to check that a macroscopic mesh can describe wrinkles even if the mesh size is larger than the wrinkling wavelength. The second example is a classical case of membrane wrinkling under tensile load [START_REF] Friedl | Buckling of stretched strips[END_REF][START_REF] Jacques | On mode localisation in tensile plate buckling[END_REF]. In this case, the number of waves is much smaller and the onset of instability is very sensitive to loading and boundary conditions.

The results are compared with commercial software ABAQUS to assess advantages and disadvantages of the new reduced model; the shell element S8R5 with reduced integration is adopted in ABAQUS and referred as 'full shell model'. 2.1. The wavenumber Q is chosen with the [START_REF] Cerda | Geometry and Physics of Wrinkling[END_REF]. Remember that there are 10 finite elements in the length for about 30 wrinkles.

Rectangular membrane under biaxial load

formula Q = 4 √ 12π 2 (1 -ν 2 ) 1 hL x 4 √ σ y E , see
This illustrates that the macroscopic mesh is not related to the size of the wrinkles. The wrinkling amplitude and shape are further studied by plotting the central cross section (Y = 100 mm). There are several ways to improve the reduced model. A simple idea would be to increase the number of harmonics. Another possibility is to go back to a complex envelope w 1 , indeed it has been proved [START_REF] Mhada | About macroscopic models of instability pattern formation[END_REF] that this improves the prediction of the boundary behavior by taking phase shift into account. A best possible improvement is a bridging technique as in [START_REF] Hu | A bridging technique to analyze the influence of boundary conditions on instability patterns[END_REF] that permits to go back to the full model in a small region near the boundaries.

Rectangular membrane under uniaxial load

The clamped rectangular membrane is stretched on short edges as shown in Fig.

(2.8), the two long edges are completely free. Geometric and material parameters are as following: L x = 200 mm, L y = 1400 mm, h = 0.05 mm, E = 70 000 MPa, ν = 0.3. The mesh details of each model are as in Table 2.2. The wavenumber Q is chosen from a full shell computation in ABAQUS.

The membrane would be stretched by two kinds of loading approaches: displacement load and stress load. Let us specify the boundary conditions along the . As shown in [START_REF] Friedl | Buckling of stretched strips[END_REF] and [START_REF] Jacques | On mode localisation in tensile plate buckling[END_REF], the buckling level depends on this dimensionless parameter and the buckling stress decrease when the ratio C increase. Because the ratio C in displacement load case is smaller than in stress load, see Table . 2.3, the displacement load case turns out to give a bigger critical load. Besides, the stress ratio C has been observed in the whole process of loading. During pre-buckling period, C remains a constant, are quite similar except slight differences in the wrinkling amplitude and distribution. In fact, one can notice a discrepancy in wavelength and phase as compared with full shell model. As mentioned before, the envelope w 1 is assumed to be real, what provides a constant wavenumber Q through the procedure. Therefore these changes can not be described by the new reduced model.

In conclusion, despite slight differences in boundary conditions, the new reduced model is able to predict accurately and fast the bifurcation point with much less degrees of freedom than full shell model, and correctly describe the post-buckling process.

A dimensionless parameter to govern wrinkling initiation

Wrinkling is characterised by an instability with a short wavelength in a direction,

x in this Chapter, in such a way that the Laplacian ∆ can be approximated by ∂ 2 ∂x 2 . Let us come back to the bending equation (2.1-d) of the Föppl-Von Karman model, with a biaxial prestress N = hσ x (y)e x ⊗ e x + hσ y e y ⊗ e y that is not necessarily constant in the direction of the wrinkles. Thus the Eq. (2.1-d) becomes:

D ∂ 4 w ∂x 4 + h|σ x (y)| ∂ 2 w ∂x 2 -hσ y ∂ 2 w ∂y 2 = 0.
(2.54)

The stress ratio C is a first dimensionless parameter that influences the wrinkling process and that is constant in the pre-buckling phase. Here we shall define another dimensionless parameter from a characteristic length l c defined as in [START_REF] Jacques | On mode localisation in tensile plate buckling[END_REF] that is related to the size of the compressed area. Thus one can define non-dimensional variables x = xQ, ȳ = y/l c , where the wrinkling wavenumber Q is unknown. This leads to the following eigenvalue problem:

DQ 4 ∂ 4 w ∂ x4 + h|σ x (y)|Q 2 ∂ 2 w ∂ x2 - hσ y l c 2 ∂ 2 w ∂ ȳ2 = 0.
(2.55)

Classically [START_REF] Cerda | Geometry and Physics of Wrinkling[END_REF] the wrinkling wavelength is obtained by requiring that the first and the third term are of the same order:

DQ 4 ∼ hσ y l c 2 ⇒ Q 2 ∼ √ σ y E 1 hl c , (2.56)
so that the critical stress |σ crit x | ≈ 2DQ 2 can be related to the quantities E, σ y , h, l c :

|σ crit x | ∼ √ Eσ y h l c . (2.57)
That is why we discuss here the influence of the non-dimensional loading parameter

K = |σ x | √ Eσ y l c h .
(2.58)

The definition of the characteristic length l c can be questionable in the same way as the characteristic length defining the Reynolds number in fluid mechanics.

This is quite clear for the clamped membrane under uniform biaxial load studied in Section 2.4.1, in which case there is an analytical solution [START_REF] Damil | New nonlinear multi-scale models for wrinkled membranes[END_REF] in a double sine shape so that the characteristic length can coincide with the plate width:

l c = L y .
Nevertheless even in this simple case, this definition depends on boundary conditions. In the uniaxial tensile cases of Section 2.4.2, the compression σ x is not uniform and the wrinkling is due to one or two previously mentioned compressive zones, but the most relevant characteristic length is not the size of these compressive zones. Indeed there is a large region where the compressive stress is very small and therefore has little influence on wrinkling. That is why we propose to define l c by the size of the region where |σ x | is greater than |σ max x |/10, see Fig. 2.12. For very short plates, there is only one compressive zone, the characteristic length l c is also defined as before (Fig. 2.14a). For short plates where the two compressive zones overlap (Fig. 2.14b), we propose to define the characteristic length as the half of this area where |σ x | is greater than |σ max x |/10. In fact, this characteristic length can be understood as a longitudinal wrinkling wave length and it is generally unknown.

For instance in [START_REF] Jacques | On mode localisation in tensile plate buckling[END_REF], such a length has been defined from a semi-analytical minimization of the bifurcation load. We tried to define it a priori from the shape of the initial stress field, but of course it is difficult to establish a universal rule.

To check the relevance of the new non-dimensional loading parameter K, several full nonlinear shell computations have been performed, for various geometries and for four loading cases: a clamped plate with biaxial load as in Fig. the sides y = 0, y = L y . Note that in these clamped-free cases, the stress σ x (y) is not uniform. The parameter K in Eq. (2.58) is defined by accounting the maximum of this stress.

Full shell finite element results are presented in

Implementation of new reduced model in ABAQUS

To simplify the simulation procedure and use it for more complicate membrane structures, the new reduced model was implemented in the commercial finite element code ABAQUS as a user element (UEL) subroutine.

In general procedures, ABAQUS uses Newton-Raphson's method as a numerical technique for solving nonlinear equilibrium equations. Path following calculations in the presence of bifurcation require the application of an arc-length technique.

We used the so-called 'modified Riks analysis'.

In this algorithm, three quantities have to be defined: the stiffness matrix [K],

the residual force vector RHS(1) that depends on the right-hand-side force vector {F} and the increment of external load RHS(2). In the case of mechanical equilibrium of a finite element subject to surface tractions {t} and body forces {f } with stress [σ], the RHS force vector is defined as

{F} = ∫ S [N] • {t}dS + ∫ V [N] • {f }dV - ∫ V {β} : [σ]dV, (2.59)
where [N] and {β} are interpolation matrices.

In modified Riks static analysis, the implemented UEL updates the right-handside vectors RHS(1) = {F} (residual force vector)and RHS(2) = △λ(∂{F}/∂{λ}) (increments of external loads on the element) of the overall system of equations. In this Chapter, as the body force is neglected, the residual force vector {F} could be updated as (2.60) in which {f } e represents the total applied loads that must be passed into UEL as the distributed loads by array ADLMAG. The residual increment RHS(2) could be passed in by array DDLMAG that contains the increments of the distributed loads that are currently active.

{F} = {f } e - ∫ ∫ ωe t [B(q)] {N} + t [G] t [T β ][L f ][T β ][G] {q} e dω,
Last elementary stiffness matrix [K] could be derived on the basis of previous work and expressed as

[K] = ∫ ∫ ωe t [B(q)][L m ][B(q)] + t [G][ N][G] + t [G] t [T β ][L f ][T β ][G]dω, (2.61)
which is the same with Eq. (2.44) in Section 2.3.5 and passed into UEL by array AMATRX.

To verify the user element subroutine code, the clamped rectangular membrane under biaxial load of Section 2.4.1 is tested, the parameters remaining the same.

Conclusions

A new reduced-order finite element model has been introduced and implemented in a home-made code and in the commercial package ABAQUS via the user element subroutine UEL. The evolution of the wrinkles is defined in a macroscopic way, by a finite element discretization of some envelopes of the wrinkling patterns. Thus one can use coarse meshes, the size of the elements being related to the variations of the envelope and not to the details of wrinkles. This is confirmed by our numerical experiments and the reduced model needs much less degrees of freedom and computation time than the full shell model. Two numerical examples have been analyzed

and the new reduced model is able to describe correctly and quickly the details of membrane wrinkling, even if the number of wrinkles is small. Especially the classical problem of a rectangular membrane under uniaxial load has been studied in various loading cases and the boundary conditions have a tremendous effect on the stability of the structure. The numerical results show the strong influence of a dimensionless loading parameter on the bifurcation and post-bifurcation behavior.

The present reduced-order model is based on the method of Fourier series with variable coefficients that has been widely used for several structural models: beam on foundation [START_REF] Damil | Influence of local wrinkling on membrane behaviour: a new approach by the technique of slowly variable Fourier coefficients[END_REF][START_REF] Mhada | About macroscopic models of instability pattern formation[END_REF][START_REF] Hu | A bridging technique to analyze the influence of boundary conditions on instability patterns[END_REF], sandwich models [START_REF] Liu | A new Fourierrelated double scale analysis for instability phenomena in sandwich structures[END_REF][START_REF] Yu | Multi-scale techniques to analyze instabilities in sandwich structures[END_REF] or 2D hyperelasticity. Sometimes the behavior near the boundary is not perfectly accounted by the reduced model, but it is possible to re-introduce locally the full model, for instance by using the Arlequin coupling method [START_REF] Hu | A bridging technique to analyze the influence of boundary conditions on instability patterns[END_REF][START_REF] Yu | Multi-scale techniques to analyze instabilities in sandwich structures[END_REF][START_REF] Dhia | Multiscale mechanical problems: the Arlequin method[END_REF], 2006;[START_REF] Hu | Multi-scale nonlinear modelling of sandwich structures using the Arlequin method[END_REF]. Another possible variant is to introduce more harmonics, for instance as in [START_REF] Mhada | About macroscopic models of instability pattern formation[END_REF] where a complex envelope model has been applied, but compromises between speed and accuracy have to be found.

Abstract

In this Chapter, the technique of the slowly varying Fourier coefficient introduced in Chapter 2 is extended for the study of local wrinkling in two-dimensional film/substrate system, in which, the displacement field is transformed into the slowly variable Fourier coefficient, i.e., the macroscopic displacement field that permits to capture the wrinkling evolution in the system with much less degrees of freedom than the full finite element model. The derived macroscopic non-linear system is solved by the Asymptotic Numerical Method that is very efficient and reliable to capture the bifurcation point and the post-buckling path in wrinkling analyses. In particular, the importance of using the first harmonic of Fourier series in approximating the axial displacement in substrate is discussed and a spurious phenomenon related to the hypothesis of the used approximation functions within the Fourier-series approach, i.e., oscillation locking, is pointed out. To overcome this phenomenon, modifications on either the Fourier series or the constitutive equations of the substrate are proposed. The efficiency and accuracy of the proposed macroscopic model are demonstrated by the wrinkling simulations for several kinds of film/substrate systems.

Present Chapter corresponds to the published research paper (Huang et al., 2017) [ Thin-Walled Structures, 114: 144-153, 2017.]. A self-consistent notation is adopted.

Keywords: Fourier series, Film/substrate system, Wrinkling, Asymptotic Numerical Method.

Introduction

Compressed stiff films bonded to a compliant substrate can buckle into a pattern presenting sinusoidal wrinkles with constant wavelength [START_REF] Chen | Herringbone buckling patterns of compressed thin films on compliant substrates[END_REF]; [START_REF] Huang | Nonlinear analyses of wrinkles in a film bonded to a compliant substrate[END_REF]) when compression reaches a critical value. The wrinkles may have an undesirable effect on the system and should be often avoided.

However, the periodic nature of the wrinkles in film/substrate systems has nowadays found some applications such as buckling-based metrology method [START_REF] Stafford | A buckling-based metrology for measuring the elastic moduli of polymeric thin films[END_REF]), optical gratings [START_REF] Harrison | Sinusoidal phase grating created by a tunably buckled surface[END_REF]) and stretchable electronics [START_REF] Rogers | Materials and Mechanics for Stretchable Electronics[END_REF]). Since the critical compressive load and details of the instability pattern are often of interest, an accurate yet efficient characterisation of wrinkling of film/substrate system is necessary.

Wrinkling in film/substrates system is very similar to the symmetrical wrinkling in sandwich panel, in which wrinkling stresses for three modes (single sided face wrinkling, in-phase wrinkling and out-of-phase wrinkling) are expressed by a unified, single expression through approximate structural analysis [START_REF] Allen | Analysis and Design of Structural Sandwich Panels[END_REF]; [START_REF] Niu | Modeling of wrinkling in sandwich panels under compression[END_REF]). Recently, some non-linear analyses have been developed

to comprehend and characterise the wrinkle morphologies. The stiff thin films are usually modelled as a non-linear elastic beam or plate in the works [START_REF] Bowden | Spontaneous formation of ordered structures in thin filmsofmetals supported on an elastomeric polymer[END_REF]; [START_REF] Chen | Herringbone buckling patterns of compressed thin films on compliant substrates[END_REF]; [START_REF] Huang | Nonlinear analyses of wrinkles in a film bonded to a compliant substrate[END_REF]; [START_REF] Song | Buckling of a stiff thin film on a compliant substrate in large deformation[END_REF]). Differences in these works mainly arise from the used kinematics for the substrate and whether the shear deformation at the film/substrate interface is considered or not. [START_REF] Chen | Herringbone buckling patterns of compressed thin films on compliant substrates[END_REF] modelled the substrate as an elastic foundation of springs where the in-plane displacements are ignored. Results show that the herringbone mode constitutes a minimum energy configuration. [START_REF] Huang | Nonlinear analyses of wrinkles in a film bonded to a compliant substrate[END_REF] adopted the three-dimensional elastic field for the substrate and investigated the influence of the Young's modulus, the Poisson's ratio and the thickness of the substrate on the critical strain, amplitude and wavelength of the sinusoidal wrinkles.

In their works, the shear deformation in the film/substrate interface was neglected.

Mei et al. ( 2011) assumed a continuous shear strain across the film/substrate interface. This showed that a significant error of the critical strain and the wrinkling wavelength emerges in [START_REF] Huang | Nonlinear analyses of wrinkles in a film bonded to a compliant substrate[END_REF] when the Poisson's ratio of the substrate decreases. In the above mentioned works, linear elastic constitutive law was assumed in the substrate. [START_REF] Song | Buckling of a stiff thin film on a compliant substrate in large deformation[END_REF] considered finite strains and a non-linear constitutive law in the substrate and showed that the wrinkling wavelength decreases with the increase of the pre-strain rather than remaining constant as in [START_REF] Huang | Nonlinear analyses of wrinkles in a film bonded to a compliant substrate[END_REF]. Furthermore, [START_REF] Brau | Multiple-length-scale elastic instability mimics parametric resonance of nonlinear oscillators[END_REF] pointed out that a quadratic non-linearity of the substrate can trigger the period-doubling instability that cannot be captured by a linear substrate. Hutchinson ( 2013) also studied the influence of the quadratic non-linearity of a neo-Hookean substrate on the mode evolution of the wrinkles in film/substrate systems. [START_REF] Shariyat | 3D energy-based finite element elasticity approach for shear postbuckling analysis of functionally graded plates on elastic foundations[END_REF] recently investigated the instability of rectangular FGM plate on elastic foundation under shear loading and found that the stiffness of the elastic foundation greatly affects the angles of deformation waves. Taking into account both accuracy and computational efficiency, higher order functions can be used to approximate kinematics of the substrate and thus to reduce the number of unknowns, see [START_REF] Reddy | A simple higher-order theory of laminated composite plate[END_REF]; [START_REF] Vidal | A family of sinus finite elements for the analysis of rectangular laminated beams[END_REF]; [START_REF] Ferreira | Analysis of thick isotropic and cross-ply laminated plates by Generalized Differential Quadrature Method and a unified formulation[END_REF]; [START_REF] Tornabene | Radial basis function method applied to doubly-curved laminated composite shells and panels with a general higher-order equivalent single layer theory[END_REF][START_REF] Tornabene | Higher-order structural theories for the static analysis of doubly-curved laminated composite panels reinforced by curvilinear fibers[END_REF]. [START_REF] Yang | A new family of finite elements for wrinkling analysis of thin films on compliant substrates[END_REF] proposed a higher-order kinematics to model substrate based on the Carrera's Unified Formulation (CUF) (Carrera (2003a); [START_REF] He | Multi-scale modelling of sandwich structures using hierarchical kinematics[END_REF]; Cinefra et al. ( 2012)). The finite element method was used to solve the higher-order model and accurate results were obtained with low computational cost.

In this Chapter, stemming from an efficient multi-scale approach established by Damil andPotier-Ferry (2010, 2008) that exploits the periodic nature of the wrinkles, an effective two-dimensional Fourier-based model is developed to study the sinusoidal wrinkling in thin stiff films on compliant substrates. The problem unknowns arising from the assumed kinematics (here addressed as "microscopic model" where Euler-Bernoulli's beam theory is used for the film and a twodimensional plane-strain solid for the substrate) are expanded by Fourier series, which leads to a new problem with the Fourier coefficients as new unknowns exhibiting much slower variation than the original unknowns. This latter problem is called "macroscopic model". The derived macroscopic model has the advantage to require very few degrees of freedom to accurately describe the problem under study.

Compared to the Landau-Ginzburg technique, the Fourier-based method has two advantages: 1) not only the bifurcation point but also the post-buckling path can be captured and 2) the coupled global and local instability patterns can be incorporated and characterised, see [START_REF] Liu | A new Fourierrelated double scale analysis for instability phenomena in sandwich structures[END_REF]. The established nonlinear system shows strong nonlinearity near and after bifurcation point, and it is a difficult task to solve this kind of nonlinear problem. In this Chapter, the Asymptotic Numerical Method (ANM), known as an effective and robust nonlinear solver in tracing bifurcation path in instability problems proposed by [START_REF] Damil | A new method to compute perturbed bifurcations: Application to the buckling of imperfect elastic structures[END_REF]; [START_REF] Cochelin | Méthode asymptotique numérique[END_REF][START_REF] Cochelin | Asymptotic-numerical methods and Padé approximants for non-linear elastic structures[END_REF], is used to solve the established nonlinear equations.

So far, the approach based on Fourier series has been successfully used for the wrinkling analysis of non-linear beams resting on Winkler's foundation [START_REF] Damil | Influence of local wrinkling on membrane behaviour: a new approach by the technique of slowly variable Fourier coefficients[END_REF]; [START_REF] Mhada | About macroscopic models of instability pattern formation[END_REF]), sandwich beams [START_REF] Yu | Multi-scale techniques to analyze instabilities in sandwich structures[END_REF]; Mhada et al. ( 2013)) and thin metal films [START_REF] Damil | New nonlinear multi-scale models for wrinkled membranes[END_REF][START_REF] Damil | Membrane wrinkling revisited from a multiscale point of view[END_REF]; [START_REF] Huang | Macroscopic simulation of membrane wrinkling for various loading cases[END_REF]; [START_REF] Attipou | Thermal wrinkling of thin membranes using a Fourier-related double scale approach[END_REF]). For the above Fourier-related models, the neutral axial displacement is always assumed not fluctuating and the first order harmonics are ignored, which matches the practical kinematics in these beam or plate elements.

However, in the current film/substrate model, the substrate is discretized by continuum elements and shows oscillation at the film/substrate interface. Therefore, the first harmonic of the axial displacement should be introduced to characterise such oscillation. Otherwise, a spurious stiffening effect (named as "oscillation locking") would occur if only the zero order of Fourier series for the axial displacement is used as introduced in previous works. To overcome this difficulty, two manners are proposed: 1) the first one is to enrich the Fourier-series expansion and 2) the second one is to modify the elastic coefficient of the substrate.

This Chapter is structured as follows. The two-dimensional microscopic model for film/substrate systems is introduced in Section 3.2. In Section 3.3, the Fourierrelated macroscopic model and the corresponding finite element are derived. In Section 3.4, numerical tests are performed to assess the established macroscopic model, and discussion on the importance of the first harmonic in continuum element and the oscillation locking is made.

Microscopic model

Kinematics

A two-dimensional elastic stiff film bound to an elastic soft compliant substrate is considered as depicted in Fig. 3.1: x and z are the axial and through-the-thickness coordinates, h f and h s are the thickness of the top film and the substrate, respectively, h t is the total thickness of the structure. The length and the width are denoted respectively by L and b. The thin film is modelled as an elastic Euler-Figure 3.1: Sketch of an elastic thin stiff film on a compliant substrate.

Bernoulli's beam:

U f (x, z) = u f (x) - ( z - h f + 2h s 2 ) w f ,x (x) W f (x, z) = w f (x) z ∈ [h s , h t ] (3.1)
where superscript "f " stands for the film, U (x, z) and W (x, z) are the axial and through-the-thickness components of the displacement field U (x, z) and u f and w f the two unknown displacement functions. A coordinate subscript preceded by comm ", x" stands for a partial derivative. The substrate, denoted by "s", is modelled as a plane-strain elastic solid:

U s (x, z) = u s (x, z) W s (x, z) = w s (x, z) z ∈ [0, h s ] (3.2)
The displacement field at interface of the film and substrate should satisfy the following compatibility conditions:

U f (x, h s ) = U s (x, h s ) W f (x, h s ) = W s (x, h s ) (3.3)
The Constrained Variational Principle (CVP) is used to ensure the congruency of the displacement field at the interface:

Γ = {(x, y, z) : x ∈ [0, L] , y ∈ [-b/2, b/2] , z = h s } (3.4)
The following constrain term L (µ, u) is obtained by introducing the Lagrange multipliers µ = {µ 1 , µ 2 } as fictitious gluing forces:

L (µ, u) = ∫ Γ { µ 1 [ U f (x, h s ) -U s (x, h s ) ] + µ 2 [ W f (x, h s ) -W s (x, h s ) ]} dΓ (3.5)

Geometric equations and constitutive law

The geometric equations and constitute law are supposed to meet the following hypotheses:

1. the material behaviour is linear elastic and described by the Hooke's law, 2. the geometrical non-linearity is considered in the film only and the straindisplacement relationship is described by the Green-Lagrange strain, 3. the substrate undergoes small displacements.

The above assumptions are translated into the following equations:

σ f xx = E f ϵ f xx = E f [ U f ,x + 1 2 ( W f ,x ) 2 ] σ s xx = C s 11 ϵ s xx + C s 13 ϵ s zz = C s 11 U s ,x + C s 13 W s ,z σ s zz = C s 13 ϵ s xx + C s 33 ϵ s zz = C s 13 U s ,x + C s 33 W s ,z σ s xz = 2C s 55 ϵ s xz = C s 55 ( U s ,z + W s ,x ) (3.6)
with ϵ ij and σ ij being the strain and stress components, C ij the material stiffness coefficients

C s 11 = C s 33 = 1 -ν s (1 + ν s ) (1 -2ν s ) E s C s 13 = ν s (1 + ν s ) (1 -2ν s ) E s C s 55 = E s 2 (1 + ν s ) (3.7)
E the Young's modulus and ν the Poisson's ratio.

The weak form of the governing equations is obtained by the Principal of Virtual Displacement accounting for the constraint in Eq. (3.5):

δP int + δL = δP ext (3.8)
where δ stands for a virtual variation, P int is the internal work and P ext the external work. The internal virtual work is expressed as:

δP int = ∫ V f σ f xx δϵ f xx dV + ∫ V s (σ s xx δϵ s xx + σ s zz δϵ s zz + 2σ s xz δϵ s xz ) dV (3.9)
where the V stands for volume. By using Eqs. (3.1) to (3.7) and integrating through the thickness of the film, the virtual internal work in Eq. (3.9) rewrites:

δP int = ∫ L 0 δγ f T S f dx + b ∫ ω s δγ sT S s dω (3.10)
where the superscript "T" stands for transposition and ω s = L × h s ∈ R 2 is the substrate integration domain on the plane xz. S f and S s are the film and substrate stress vectors:

S f = D f [ H f + 1 2 A ( θ f ) ] θ f S s = D s H s θ s (3.11)
in which

D f =   E f h f b 0 0 1 12 E f bh 3 f   H f = [ 1 0 0 0 0 1 ] A ( θ f ) = [ 0 w f ,x 0 0 0 0 ] (3.12) and D s =     C s 11 C s 13 0 C s 13 C s 33 0 0 0 C s 55     H s =     1 0 0 0 0 0 0 1 0 1 1 0     (3.13) θ is the displacement gradient θ f T = { u f 0,x w f ,x w f ,xx } θ sT = { u s ,x u s ,z w s ,x w s ,z } (3.14)
and γ the generalised strain vector.

Macroscopic model

The double-scale approach based on Fourier series with slowly varying coefficients is applied to derive the film/substrate macroscopic model. The characteristic direction and period are described by a wave-number q that is assumed to be given. The unknown displacement field U (x, z) is written in the following form:

U (x, z) = +∞ ∑ j=-∞ U j (x, z) e ijqx (3.15)
where i is the imaginary unit and the Fourier coefficient U j (x, z) denotes the envelope of the jth order harmonic that is conjugated with U -j (x, z). U j varies slowly in the axial direction over a period

[ x, x + 2π q
] of the oscillation. The wave-number q is defined as

q = π L q 0 (3.16)
in which q 0 is a chosen integer wave-parameter: q 0 ∈ N and q 0 ≥ 2. Here a macroscopic model with three envelops U 0 (x, z), U 1 (x, z) and U -1 (x, z), as pictured in Fig. 3.2, is considered. U 0 can be identified as the mean value while U ±1 represent the envelope or the amplitude of the spatial oscillations.

Furthermore, some other simplified assumptions are made. The in-plane displacement of the neutral plane of the film u f (x) and the applied external force F do not fluctuate according to [START_REF] Damil | Influence of local wrinkling on membrane behaviour: a new approach by the technique of slowly variable Fourier coefficients[END_REF]. This means that their Fourier series expansion reduces to a single term

u f = u f 0 F = F 0 (3.17)
In this Chapter, we focus on the local instability of film/substrate system, so the mean field of transverse displacement is nearly unchanged along the x direction and it leads to

W 0 (x, z) = 0 (3.18)
According to the above assumptions, the transverse displacement W can be explicitly expressed as

W (x, z) = 2 Ŵ (x, z) cos (qx + ϕ) (3.19)
where the envelop Ŵ ∈ R and the phase angle ϕ are used in place of the complex conjugate W ±1 . The phase angle ϕ could be variant to describe the phase modulation in the evolution of instability [START_REF] Mhada | About macroscopic models of instability pattern formation[END_REF]. For the sake of simplicity, ϕ is assumed to be constant and equal to ± π 2 (see [START_REF] Liu | A new Fourierrelated double scale analysis for instability phenomena in sandwich structures[END_REF] and [START_REF] Yu | Multi-scale techniques to analyze instabilities in sandwich structures[END_REF] for more details) and the previous equation becomes

W (x, z) = 2 Ŵ (x, z) sin (qx) (3.20)
By substituting Eqs. (3.17) and (3.20) into Eq. (3.1), the displacement fields of the film can be rewritten as

U f (x, z) = u f 0 (x) - ( z - h f + 2h s 2 ) ( 2qw f 1 cos(qx) + 2w f 1,x sin(qx) ) W f (x, z) = 2w f 1 (x)sin (qx) z ∈ [h s , h t ] (3.21)
Finally, some properties are introduced to obtain the macroscopic governing equations. Given two functions a (x) and b (x) expanded by Fourier series with slowly varying coefficients, the following identities hold, see [START_REF] Damil | Influence of local wrinkling on membrane behaviour: a new approach by the technique of slowly variable Fourier coefficients[END_REF]:

∫ L 0 a (x) b (x) dx = ∫ L 0 ∞ ∑ j=-∞ a j (x) b -j (x) dx (3.22) ( da dx ) j = ( a ′ ) j = ( d dx + ijq ) a j = (a j ) ′ + ijqa j (3.23) ( d 2 a dx 2 ) j = (a ′′ ) j = ( d dx + ijq ) 2 a j = (a j ) ′′ + 2ijq (a j ) ′ -j 2 q 2 a j (3.24) (ab) j = +∞ ∑ j 1 =-∞ a j 1 b j-j 1 (3.25)

Internal virtual work of the film

δP f int = ∫ L 0 +1 ∑ j=-1 δγ f T -j S f j dx = ∫ L 0 ( δγ f T 0 S f 0 + 2δγ f T 1R S f 1R + 2δγ f T 1I S f 1I ) dx (3.26)
The constitutive law of the film reads

S f 0 = D f γ f 0 S f 1R = D f γ f 1R S f 1I = D f γ f 1I (3.27)
and the strains are written as

γ f 0 = H f θ f 0 + 1 2 A ( θ f 0 ) θ f 0 + A ( θ f 1R ) θ f 1R + A ( θ f 1I ) θ f 1I γ f 1R = H f θ f 1R + A ( θ f 0 ) θ f 1R γ f 1I = H f θ f 1I + A ( θ f 0 ) θ f 1I (3.28)
The vectors S f , γ f and θ f are the Fourier series expansion of the unknowns and the subscripts "0", "1R" and "1I" stand for the zero-(real value), first-order harmonic (j = ±1, real-value part "R" and complex-value part "I").

For instance, the Fourier series expansion of θ f can be written as

θ f = +1 ∑ j=-1 θ f j e ijqx (3.29)
where, according to the hypotheses given in the previous paragraph

θ f T 0 = { u f 0,x 0 0 } θ f T 1 = { 0 w f 1,x + iqw f 1 w f 1,xx -q 2 w f 1 + i2qw f 1,x } (3.30)
θ f 1 can be also written as

θ f 1 = ℜ ( θ f 1 ) + iℑ ( θ f 1 ) = θ f 1R + iθ f 1I (3.31)
To avoid spurious oscillations (see details in [START_REF] Damil | Influence of local wrinkling on membrane behaviour: a new approach by the technique of slowly variable Fourier coefficients[END_REF]), the second derivative of the envelope w f 1 in Eq. (3.30) is eliminated. After integration by part, the second term in Eq. (3.30) reads

θ f T 1 = { 0 w f 1,x + iqw f 1 -q 2 w f 1 + i √ 6qw f 1,x } (3.32)
By introducing the following generalised stress S f gen , strain Γ f gen and displacement gradient θ f gen vectors:

S f T gen = { S f T 0 S f T 1R S f T 1I } Γ f T gen = { γ f T 0 2γ f T 1R 2γ f T 1I } θ f T gen = { θ f T 0 θ f T 1R θ f T 1I } (3.33)
The generalized macroscopic stress-strain relation can be consequently written as:

S f gen = D f gen Γ f gen Γ f gen = H f gen θ f gen + 1 2 A ( θ f gen ) θ f gen (3.34)
where

D f gen =      D f 0 0 0 1 2 D f 0 0 0 1 2 D f      H f gen =     H f 0 0 0 2H f 0 0 0 2H f     A ( θ f gen ) = 2       1 2 A ( θ f 0 ) A ( θ f 1R ) A ( θ f 1I ) A ( θ f 1R ) A ( θ f 0 ) 0 A ( θ f 1I ) 0 A ( θ f 0 )       (3.35)
The virtual work of the film finally reads:

δP f int = ∫ L 0 δΓ f T gen S f gen dx (3.36)

Internal virtual work of the substrate

The virtual work of the substrate is obtained in a similar manner:

δP s int = b ∫ ω s +1 ∑ j=-1 δγ sT -j S s j dω s = b ∫ ω s δΓ sT gen S s gen dω (3.37)
where Γ s gen and S s gen are the generalised strain and the generalised stress in the substrate, respectively. The following constitutive and geometric equations hold

S s gen = D s gen Γ s gen Γ s gen = H s gen θ s gen (3.38)
where the strain-displacement relationship is assumed to be linear. The matrices D s gen and H s gen are defined in a manner similar to D f gen and H f gen , see Eq. (3.35), and they are not reported for the sake of brevity. The generalised displacement gradient θ s gen is:

θ sT gen = { θ sT 0 θ sT 1R θ sT 1I } (3.39) in which θ sT 0 = { u s 0,x u s 0,z 0 0 } θ sT 1R = { u s 1R,x -qu s 1I u s 1R,z w s 1,x w s 1,z } θ sT 1I = { qu s 1R + u s 1I,x u s 1I,z qw s 1 0 } (3.40)
It should be noted that the first-order harmonic of the axial displacement should be accounted for, otherwise the in-plane displacement U s would vary slowly along the x direction, which will result in a spurious over-stiffening phenomenon named "oscillation locking", see details in Section 3.4.3.

The congruency of the displacement field at film and substrate interface Γ is imposed by means of the Lagrange multiplier µ and the coupling variational term L in Eq. (3.5), which is realised by coupling the envelopes of the displacement field

U f 0 (x, h s ) = U s 0 (x, h s ) U f 1 (x, h s ) = U s 1 (x, h s ) W f 0 (x, h s ) = W s 0 (x, h s ) W f 1 (x, h s ) = W s 1 (x, h s ) (3.41)
Its virtual variation thus reads

δL = L (δµ, u) + L (µ, δu) = ∫ Γ { δµ u 0 [ U f 0 (x, h s ) -U s 0 (x, h s ) ] + δµ u 1R [ U f 1R (x, h s ) -U s 1R (x, h s ) ] + δµ u 1I [ U f 1I (x, h s ) -U s 1I (x, h s ) ] + δµ w 1 [ W f 1 (x, h s ) -W s 1 (x, h s ) ] } dΓ+ ∫ Γ { µ u 0 [ δU f 0 (x, h s ) -δU s 0 (x, h s ) ] + µ u 1R [ δU f 1R (x, h s ) -δU s 1R (x, h s ) ] + µ u 1I [ δU f 1I (x, h s ) -δU s 1I (x, h s ) ] + µ w 1 [ δW f 1 (x, h s ) -δW s 1 (x, h s ) ] } dΓ = b ∫ L 0 { δµ u 0 [ u f 0 -u s 0 ] + δµ u 1R [ h f qw f 1 -2u s 1R ] + δµ u 1I [ h f w f 1,x -(-2u s 1I ) ] + δµ w 1 [ 2w f 1 -2w s 1 ] } dx+ b ∫ L 0 { µ u 0 [ δu f 0 -δu s 0 ] + µ u 1R [ h f qδw f 1 -2δu s 1R ] + µ u 1I [ h f δw f 1,x -(-2δu s 1I ) ] + µ w 1 [ 2δw f 1 -2δw s 1 ] } dx (3.42)
The matrix form of Eq. (3.42) reads:

δL (µ, u) = ∫ L 0 [ δµ T ( D f l θ f l -D s l θ s l ) + ( δθ f T l D f T l -δθ sT l D sT l ) µ ] dx (3.43)
where the subscript "l" stands for the terms related to the Lagrange multiplier, θ f l and θ s l are, respectively, the displacement gradient vectors of the film and the substrate

θ f T l = { u f 0 w f 1 w f 1,x } θ sT l = { u s 0 u s 1R u s 1I w s 1 } (3.44)
D f l and D s l are the corresponding constant matrices that can be derived conveniently from Eq. (3.42), which are not shown here for the sake of brevity.

External work and weak form of the governing equation

Due to an external load F applied on the film lateral surface, see Fig. 3, the external virtual work is supposed to be:

δP ext = λb δθ f T gen F (3.45)
where λ is a "load parameter". The macroscopic governing equation, finally, reads:

∫ L 0 [ δθ f T gen ( H f T gen + A T ( θ f gen )) S f gen ] dx + b ∫ ω s δθ sT gen D s gen θ s gen dω+ ∫ L 0 [ δµ T ( D f l θ f l -D s l θ s l ) + ( δθ f T l D f T l -δθ sT l D sT l ) µ ] dx = λb δθ f T gen F (3.46)

Finite element discretisation

The finite element method is used to solve the non-linear system in Eq. (3.46). A three-node one-dimensional element with two degrees of freedom

v ef T i = { u ef 0i w ef 1i } (3.47)
for each node i is used for the film, whereas two-dimensional eight-node elements with four degrees of freedom

v esT i = { u es 0i u es 1Ri u es 1Ii w es 1i } (3.48)
for each node i is used for the substrate. The Lagrange multiplier is discretised as the film unknowns:

v elT i = { µ e u 0i µ e u 1Ri µ e u 1Ii µ e w 1i } (3.49)
The elemental nodal unknown vectors can be written as:

v ef T = { v ef T I v ef T II v ef T III } v esT = { v esT I v esT II v esT III v esT IV v esT V v esT V I v esT V II v esT V III } v elT = { v elT I v elT II v ef T III } (3.50)
The roman numbers denote the node order in an element. For example, I indicates the first node of the element.

The vector unknowns v ef T and v elT are characterised by three-node quadratic shape functions. The serendipity eight-node shape functions are used for the unknown vector v esT . The shape functions are not presented for the sake of brevity, which can be found in [START_REF] Bathe | Finite element procedures[END_REF].

The element displacement gradient vectors θ ef gen , θ es gen , θ ef l , θ es l and the Lagrange multiplier µ e are expressed as:

θ ef gen = G f v ef θ es gen = G s v es θ ef l = G f l v ef θ es l = G s l v es µ e = G l v el (3.51)
where G f , G s , G f l , G s l and G l contain the shape functions and their derivatives, see [START_REF] Crisfield | Non-Linear Finite Element Analysis of Solids and Structures[END_REF].

Finally, the weak form of the governing equations discretised by the finite element reads

∫ le 0 δv ef T B T ( v ef ) S ef gen dx + ∫ ω s e δv esT K s v es dω+ ∫ le 0 [ δv elT ( C 1 v ef -C 2 v es ) + ( δv ef T C T 1 -δv esT C T 2 ) v el ] dx = λδv ef T f e (3.52)
where the following relations and matrix terms are introduced

S ef gen = D f gen [ B f lin + 1 2 B f nlin ( v ef ) ] v ef B ( v ef ) = B f lin + B nlin ( v ef ) B f lin = H f gen G f B nlin ( v ef ) = A ( v ef ) G f K s = G sT D s gen G s C 1 = G lT D f l G f l C 2 = G lT D s l G s l (3.53)
The resulting non-linear problem is solved by the Asymptotic Numerical Method [START_REF] Cochelin | Asymptotic-numerical methods and Padé approximants for non-linear elastic structures[END_REF]). For each element, the unknowns are determined at the nodes and the stress at the Gauss points. Details on the use of the ANM can be easily found in our papers (Liu et

Numerical results and discussions

The film/substrate system is subjected to compressive concentrated forces of equal modulus F and opposite verse as illustrated in Fig. 3.3. The forces act along the axial direction at the through-the-thickness mid point of the film. The lower Figure 3.3: Sketch of the film/substrate system under compression forces surface of the substrate (z/h t = 0) is fixed. This boundary condition is obtained by setting u s 0 (x, 0), u s 1R (x, 0), u s 1I (x, 0) and w s 1 (x, 0) equal to zero. The vertical displacements at x/L = 0 and 1 are also fixed by imposing a nil value for w s 1 (0, z), w s 1 (L, z), w f 1 (0, z) and w f 1 (L, z). The results obtained by the macroscopic model are compared with a twodimensional non-linear elasticity solution obtained by the finite element software ABAQUS. The latter is referred to as either "2D full model" or "ABAQUS". For this ABAQUS solution, the film is modelled by 3-node beam elements (B22) and the substrate by 8-node bi-quadratic plane-strain elements (CPE8R) with reduced integration [START_REF] Yang | A new family of finite elements for wrinkling analysis of thin films on compliant substrates[END_REF]). To ensure the displacement congruency at the film/substrate interface, the degrees of freedom in the beam elements are bound to the plane-strain elements by using "tie" constraints. Mesh densities are examined to ensure convergent results.

Validation

Numerical examples with different Young's modulus ratios E s /E f and thickness ratios h s /h f are performed to validate the proposed macroscopic model. This first set of analyses is addressed as "problem I". Table 3.1 presents the material properties and geometric parameters. E s /E f ranges from 10 -5 to 10 -3 representing typical Table 3.1: Material and geometric parameters of the film/substrate, problem I.

E f (10 4 MPa) E s /E f ν s L (10 -1 m) b (m) h f /L h s /h f
6.9 10 -5 -10 -3 0.3 5.0 1.0 1/200 10 -10 3 materials for stiff film layer resting on a compliant substrate [START_REF] Huang | Nonlinear analyses of wrinkles in a film bonded to a compliant substrate[END_REF].

The thickness ratio h s /h f varies between 10, for very thin substrates, and 10 3 for very thick substrates.

For the proposed macroscopic model, a mesh with 8 elements along the x direction and 40 elements in z direction in the substrate is used. This implies about 4.33 • 10 3 degrees of freedom (DOFs). Linear eigenvalue buckling analyses using function BUCKLE in the ABAQUS [START_REF] Cao | Buckling and post-buckling of a stiff film resting on an elastic graded substrate[END_REF]) are performed to obtain the critical compressive load λ c and the corresponding half-wave number q 0 of the considered film/substrate systems.

Convergence and computation cost

The model convergence is studied by varying the axial and through-the-thickness discretisation of the substrate. The same film/substrate material and geometrical parameters as in [START_REF] Yang | A new family of finite elements for wrinkling analysis of thin films on compliant substrates[END_REF] are considered and summarised in Table 3.3.

The left and right edges of the film/substrate system are fixed. This set of analyses is addressed as "problem II".

Table 3.3: Material and geometric parameters of the film/substrate, problem II. of the mesh scheme 1 (8 × 31 elements) are an order of magnitude lower than that of ABAQUS solution. Fig. 3.5 shows the bifurcation path of the transverse displacement at the point (x = L/2, z = h s + h f /2) versus the applied load for the macroscopic model I with three mesh schemes. The critical load obtained by all the considered mesh schemes converges to λ c = 0.0498 MPa and agrees well with the reference solution by ABAQUS, which shows that the macroscopic model with mesh 1 is sufficient to accurately predict the buckling load of the considered film/substrate system. Besides, the responses in the post-buckling process also match fairly well.

E f (10 5 MPa) E s /E f ν s L (10 -3 m) b (10 -3 m) h f /L h s /h f 1.8 10 -5 0.3 1.0 1.0 10 -3 10 3
The proposed macroscopic model I is able to correctly describe the film/substrate wrinkling with reduced computational cost.

Here, we further compare the computation efforts of the Fourier-based macroscopic model and the ABAQUS model in the theoretical aspect. A qualitative analysis on the computing cost of the two models is conducted by building the relations between the element number and the geometrical and material parameters in the film/substrate system.

For the ABAQUS model, the element size l f e should be proportional to the wavelength of the wrinkles λ * , i.e., l f e ∼ λ * , which means that a larger wavelengh corresponds to a larger element size. For the film/substrate system with very thick substrate in the plane-strain case [START_REF] Huang | Nonlinear analyses of wrinkles in a film bonded to a compliant substrate[END_REF]; [START_REF] Mei | Concomitant wrinkling and buckle-delamination of elastic thin films on compliant substrates[END_REF], the wavelength can be scaled by

λ * ∼ h f ( Ēf Ēs ) 1/3 (3.54)
where Ē = E/(1 -ν 2 ) is the plane-strain modulus. The element aspect ratio is assumed equal to one [START_REF] Yu | A novel two-dimensional finite element to study the instability phenomena of sandwich plates[END_REF]. Then, the element number (approximate to For the proposed macroscopic model, several elements are generally sufficient to well characterize the wrinkles in the x direction since only the slowly varying amplitude of the wrinkles needs to be simulated. Therefore, the element number in the x direction, n mx e , can be assumed as a constant independent on the geometrical and material parameters [START_REF] Huang | Macroscopic simulation of membrane wrinkling for various loading cases[END_REF]). The most important factor to determine the problem size for the macroscopic model should be the vertical displacement in the substrate, W s , which shows a sharp but monotonic variation in the area near the film. According to our previous study [START_REF] Yang | A new family of finite elements for wrinkling analysis of thin films on compliant substrates[END_REF]), the thickness of such an area with sharp-varying in W s is about λ * . Therefore, the element size in the z direction can be scaled by l mz e ∼ λ * and corresponding element number satisfies n mz e ∼ h s /λ * . Following the above definitions, the total element number follows

n m e = n mx e • n mz e ∼ h s λ * (3.56)
Let n f e and n m e be the approximations of the computation cost for the two modeling methods. Then, the ratio of computation cost should satisfy

r = n f e n m e ∼ L λ * ∼ L h f ( Ēs Ēf ) 1/3 (3.57)
We can easily find that for a larger slenderness ratio L h f or modulus ratio

Ēs

Ēf , the macroscopic model will show better computation performance than the ABAQUS model.

Discussions on Fourier-based approach and oscillation locking

The first harmonic of Fourier series

The Fourier-based approach recently has been efficiently used in analysing instabilities with spatially periodic patterns in thin structures such as the beam on

Winkler's foundation [START_REF] Damil | Influence of local wrinkling on membrane behaviour: a new approach by the technique of slowly variable Fourier coefficients[END_REF]; [START_REF] Mhada | About macroscopic models of instability pattern formation[END_REF]) and the thin metal sheet [START_REF] Damil | New nonlinear multi-scale models for wrinkled membranes[END_REF][START_REF] Damil | Membrane wrinkling revisited from a multiscale point of view[END_REF]; [START_REF] Huang | Macroscopic simulation of membrane wrinkling for various loading cases[END_REF]; Attipou et al.

(

) 2015
). In these structures, the Fourier-based models are usually discretized with structural elements, i.e., beam element and plate element, in which the neutral axial displacements are always assumed not to fluctuate along the x direction within the Fourier series approach, i.e., u f = u f 0 (see Eq. (3.17)). The assumption is valid for these structural elements since the fluctuation in the area such as the bottom or top face of these structures can be described by the derivatives of deflection W f (see Eq. (3.21)). Thus, in the Fourier-series expansion of the neutral axial displacement, the first harmonic representing the fluctuation is overlooked in these Fourier-based models.

However, for the film/substrate system that is composed of two kinds of elements (structural elements in the film and continuum elements in the substrate), it is necessary to take account for the first harmonic of the axial displacement in the substrate. Due to the compatibility of axial displacements at the film/substrate interface, when wrinkles occur in the film/substrate system, the oscillating bottom face of the film induces the occurrence of oscillation in the substrate. However, the zero harmonic U s 0 representing the mean field of the axial displacement can not approximate this oscillation, thus the first harmonic U s 1 should be taken into account, see Fig. If the first harmonic is not considered in the substrate, in another word, the same assumption as made in the previous Fourier-based beam and plate models is used, i.e., U s = U s 0 , the axial displacement continuity at the interface is then not satisfied: the axial displacement of the film's bottom face fluctuates, while the substrate does not. Therefore, the substrate becomes stiffer than actually it is and a higher external load is needed to trigger the buckling. In Fig. 3.7, the critical load predicted by the macroscopic model II is about 12% higher than that obtained by the macroscopic model I or 2D full model. The macroscopic model II is the macroscopic model that considers only the mean field U s 0 of axial displacement and uses planestrain constitutive law in the substrate. This further certify the importance of the first harmonic in approximating the axial displacement in the continuum element.

Oscillation locking

The previously mentioned over-stiffening behavior can be found in locking phenomena such as the membrane locking (Reddy (2004a)) or the Poisson thickness locking (Carrera and Brischetto (2008a,b)). The occurrence of these phenomena are related to the assumed computational approximations. For example, the Poisson thickness locking is caused by the over-simplified kinematic assumption in plate or shell analysis. In this Chapter, the over-stiffening phenomenon in the substrate leading to higher critical load is also caused by an over-simplification of the kinematics. Within the Fourier-series approach, we name this kind of phenomenon as "oscillation locking".

To avoid the oscillation locking, one can use higher order functions, i.e., the first harmonic of the axial displacement in the substrate U s = U s 0 + U s 1 e iqx + U s -1 e -iqx . This is the way to derive the macroscopic model I in Section 3 and the accuracy of this model has been validated in Section 4.1 and 4.2. In addition, we also found an effective way to alleviate the locking problem by modifying the elastic coefficient of the substrate when only the mean field of axial displacement is considered in the substrate. This manner is similar to the way dealing with thickness locking in Carrera and Brischetto (2008a,b). A detailed discussion on the alleviation of oscillation locking is presented below.

In the substrate, according to the assumption U s = U s 0 , the axial displacement shows no fluctuation and varies almost linearly (see the dotted line in Fig. 3.6, thus the axial strain ϵ xx should be almost a constant. In the plane-strain constitutive relation of the substrate, the axial strain ϵ xx has a contribution (ϵ ν xx ) due to the Poisson effect

ϵ ν xx ∝ νϵ zz (3.58)
Rapid fluctuation of the transverse normal strain ϵ zz along the x-direction causes the related ϵ ν xx to rapidly oscillate. This is in contradiction with the kinematic assumption of the macroscopic model, i.e., ϵ xx is constant. Similar to the second remedy in treating the thickness locking in Carrera and Brischetto (2008a), we modify the elastic coefficient by forcing the axial strain ϵ xx to meet the kinematic assumption of the macroscopic model. Thus, the rapid fluctuating axial strain ϵ ν xx should be eliminated, which is achieved by setting a nil Poisson ratio in the substrate. The modified elastic coefficients become

C s 11 = C s 33 = E s C s 13 = 0 C s 55 = G (3.59)
where G is the shear modulus: G = E s /2(1 + ν s ). It leads to a reduction in the elastic coefficients and thus reduces the stiffness of the substrate. The macroscopic model with only U s 0 for the in-plane axial displacement and the modified elastic coefficient is referred as "macroscopic model III" in the following.

The critical load and the half-wave number obtained by the macroscopic model III are presented in Table 3.5 for problem I. Good result agreement is observed com-pared with the reference solutions in Table 3.2. The bifurcation curves for different Table 3.5: Half-wave number q 0 and critical load parameter λ c for problem I via the macroscopic model III with elastic coefficient correction. Reference ABAQUS results are presented in Table 3.2. two-dimensional macroscopic model III is fairly close to the reference. Compared with the macroscopic model II, the oscillation locking is obviously alleviated in the macroscopic model III, showing that the reduction of elastic coefficient is also an effective approach. Besides, the degrees of freedom of this solution (1662 DOFs) are almost half of those in the two-dimensional macroscopic model I with a first harmonic term in the axial displacement of the substrate (3290 DOFs). Compared to the first remedy (using higher order approximation functions), this method can save considerable computational cost (the model has less unknowns to be determined).

h s h f 10 10 2 10 3 E s E f q 0 10 -4 λ c (MPa) q 0 10 -4 λ c (MPa) q 0 10 -4 λ c (MPa)

Conclusion

A two-dimensional Fourier-based finite element model has been introduced to study the wrinkling in thin stiff films on compliant substrates. The proposed model is able to accurately and quickly characterize the wrinkling phenomena in film/substrate systems with only several elements in the wrinkle direction. Unlike the 2D full model, the size of the elements in the wrinkle direction is independent on the wrinkling wavelength. Numerical examples have been performed to validate and assess the proposed model by determining the critical wrinkling loads and half-wave number. Results prove that the Fourier-series model yields solutions very close to the reference ABAQUS ones at a reduced computational cost. Introducing the first harmonic of the Fourier series expansion of axial displacement in continuum element has been proved to be important in approximating the actual kinematical in substrate. Furthermore, a spurious over-stiffening phenomenon named as "oscillation locking" is pointed out within the Fourier series approach, which is caused by an over-simplification of actual kinematics, i.e., the elimination of the first harmonic.

Two methods are proposed to avoid the locking effect. The first one stems from a physical consideration and consists in introducing the first harmonic of in-plane displacement in the substrate. The second one is to modify the elastic coefficient of the substrate so as to eliminate the oscillating axial strain originating from Poisson effect and match with the kinematic assumption in Fourier-series approach.

The latter is recommended as far as the computational costs are concerned. As a final remark, it is well known that, as any reduced model, the proposed model may be not accurate near the boundary. A possible solution could be the use of a multi-scale model where a full model is considered close to the boundary and the Fourier-series model in the bulk, see [START_REF] Hu | A bridging technique to analyze the influence of boundary conditions on instability patterns[END_REF]. The two models are bridged by the Arlequin method [START_REF] Hu | A bridging technique to analyze the influence of boundary conditions on instability patterns[END_REF][START_REF] Hu | Multi-scale nonlinear modelling of sandwich structures using the Arlequin method[END_REF]; Ben [START_REF] Ben Dhia | The Arlequin method as a flexible engineering design tool[END_REF]; Hu et al. (2009b)). This approach proves to be accurate yet computationally efficient to study the instabilities in sandwich structures, see [START_REF] Yu | Multi-scale techniques to analyze instabilities in sandwich structures[END_REF].

analyses.

Present Chapter corresponds to the published research paper (Huang et al., 2017) [Composite Structures, 160: 613-624, 2017.]. A self-consistent notation is adopted.

Keywords: Fourier series, CUF, Film/substrate, Wrinkling, Asymptotic Numerical Method

Introduction

System of a stiff layer resting on soft substrate exists widely in nature (human skin [START_REF] Cerda | Geometry and Physics of Wrinkling[END_REF])) and engineering fields (intelligent materials, biomedical field [START_REF] Li | Surface wrinkling of mucosa induced by volumetric growth: Theory, simulation and experiment[END_REF])). When subjected to in-plane compression or exposed to the reduction of temperature field [START_REF] Bowden | Spontaneous formation of ordered structures in thin filmsofmetals supported on an elastomeric polymer[END_REF]), the mismatch of mechanical properties or coefficients of thermal expansion of layers may lead to wrinkles. The wrinkles may pose a limit on the performance of materials or structures [START_REF] Allen | Analysis and Design of Structural Sandwich Panels[END_REF]) and are always thought to be avoided, but nowadays may find some applications such as assembly of materials [START_REF] Li | Mechanics of morphological instabilities and surface wrinkling in soft materials: a review[END_REF]), measuring the mechanical properties of materials in modern metrology [START_REF] Stafford | A buckling-based metrology for measuring the elastic moduli of polymeric thin films[END_REF]).

For these reasons, it is quite necessary to characterize the wrinkling of the system in an accurate and efficient way.

As early as forty years ago, the stability analysis of multi-layered materials was investigated by [START_REF] Allen | Analysis and Design of Structural Sandwich Panels[END_REF] in the framework of sandwich panel designs in airplanes. From then on, much of theoretical work on substrate-bonded films has been proposed on the basis of linear perturbation analysis (see [START_REF] Niu | Modeling of wrinkling in sandwich panels under compression[END_REF]), which focused on determining the critical membrane force and wavelength of wrinkles. Recently, some nonlinear analyses were performed to further investigate and comprehend wrinkle characteristics in film/substrate systems. Using energy calculation, [START_REF] Chen | Herringbone buckling patterns of compressed thin films on compliant substrates[END_REF] showed that the herringbone mode constitutes a minimum energy configuration among computed modes for a film under equal biaxial compression. [START_REF] Huang | Evolution of wrinkles in hard films on soft substrates[END_REF] modeled the substrate as an array of spring dashpots and developed a spectral method to inspect the evolution of wrinkling configurations, whose results showed that the anisotropy of the film forces affects the evolution of wrinkling patterns. As an extension work, [START_REF] Huang | Nonlinear analyses of wrinkles in a film bonded to a compliant substrate[END_REF] extended their original spring model to a model representing the substrate by three-dimensional elastic field and investigated the influence of the Young's modulus, the thickness of substrate on the amplitude and wavelength of the sinusoidal wrinkles. Through experiments and analytical solutions, [START_REF] Jiang | Finite width effect of thin-films buckling on compliant substrate: Experimental and theoretical studies[END_REF] found the width effect on wrinkles of film/substrate that the amplitude and wavelength of wrinkles increase with the film width. Unlike the results obtained by the small deformation theory [START_REF] Chen | Herringbone buckling patterns of compressed thin films on compliant substrates[END_REF]; [START_REF] Huang | Nonlinear analyses of wrinkles in a film bonded to a compliant substrate[END_REF][START_REF] Huang | Evolution of wrinkles in hard films on soft substrates[END_REF]; Jiang et al. ( 2008)), [START_REF] Song | Buckling of a stiff thin film on a compliant substrate in large deformation[END_REF] considered finite strain and nonlinear constitutive law in the substrate, and found that the buckling wavelength does not keep constant as in [START_REF] Chen | Herringbone buckling patterns of compressed thin films on compliant substrates[END_REF]; [START_REF] Huang | Nonlinear analyses of wrinkles in a film bonded to a compliant substrate[END_REF][START_REF] Huang | Evolution of wrinkles in hard films on soft substrates[END_REF]) but relates to the large prestrain. To model the film/substrate system more efficiently, [START_REF] Yang | A new family of finite elements for wrinkling analysis of thin films on compliant substrates[END_REF] proposed a high order model based on the concept of the Carrera's Unified Formulation (CUF) (Carrera (2003a); He et al. ( 2011)), and results verified that the CUF model yielded accurate results with low computational cost. Similar methods based on such enriched kinematical functions have been widely used in the mechanical analysis of composite plates [START_REF] Reddy | A simple higher-order theory of laminated composite plate[END_REF]; [START_REF] Vidal | A family of sinus finite elements for the analysis of rectangular laminated beams[END_REF]; [START_REF] Polit | High-order triangular sandwich plate finite element for linear and non-linear analyses[END_REF]; [START_REF] Ferreira | Radial basis functions collocation and a unified formulation for bending, vibration and buckling analysis of laminated plates, according to a variation of Murakami's zig-zag theory[END_REF]).

In the above film/substrate systems, instability patterns are usually nearly periodic in spatial. Based on an efficient multi-scale approach established by Damil andPotier-Ferry (2010, 2008) that exploits such periodic nature and the effective CUF model of [START_REF] Yang | A new family of finite elements for wrinkling analysis of thin films on compliant substrates[END_REF], a Fourier-related one-dimensional model is developed to study the wrinkling phenomena in thin stiff films on compliant substrates. In the framework of the multi-scale approach (often referred to as "Fourier series with slowly varying coefficients"), the unknowns to be solved in [START_REF] Yang | A new family of finite elements for wrinkling analysis of thin films on compliant substrates[END_REF] are expanded by the Fourier series, which leads to the Fourier coefficients as the new unknowns varying much more slowly than the former unknowns (we address the former unknowns as "microscopic" and the latter "macroscopic" in this Chapter).

Usually using rather coarse spatial meshes is sufficient to describe the slowly varying envelops, therefore computational efficiency can be improved significantly [START_REF] Damil | Influence of local wrinkling on membrane behaviour: a new approach by the technique of slowly variable Fourier coefficients[END_REF]; [START_REF] Liu | A new Fourierrelated double scale analysis for instability phenomena in sandwich structures[END_REF]; [START_REF] Huang | Macroscopic simulation of membrane wrinkling for various loading cases[END_REF]). This approach is firstly applied to study the instability phenomena in beams on Winkler foundation (Damil andPotier-Ferry (2010, 2008)) with a most simplified macroscopic model, in which only the zero harmonic (mean field) and the first order harmonics (envelops) are taken into account. Results showed that the macroscopic model has two advantages compared to the Landau-Ginzburg technique: 1) not only the bifurcation point but also the post-buckling path can be captured and 2) the coupled global and local instability patterns can be incorporated and characterised, see [START_REF] Liu | A new Fourierrelated double scale analysis for instability phenomena in sandwich structures[END_REF]. In this Chapter, a similar most simplified macroscopic model is firstly presented to investigate the instability of the film/substrate structures. As like any reduced models, the Fourier models may loss some accuracy in the vicinity of boundary. To handle this issue, [START_REF] Mhada | About macroscopic models of instability pattern formation[END_REF] took into account phase change in macroscopic model, and [START_REF] Hu | A bridging technique to analyze the influence of boundary conditions on instability patterns[END_REF] proposed to establish multi-scale model containing microscopic model near the boundary and macroscopic model in the bulk, the two models were bridged by the Arlequin method (Ben [START_REF] Ben Dhia | The Arlequin method as a flexible engineering design tool[END_REF]). Results showed that the multi-scale model treats boundary conditions better. For the above Fourier models, only two harmonics (the first order envelops and the mean field) were taken into account, thus the wrinkling wavelength or wavenumber should be a prior set variable, which may lead to some inconvenience by trying different wavelengths or wavenumbers to find the lowest critical buckling load. In this Chapter, we propose another strategy to trace the lowest critical load by deriving a macroscopic model considering several harmonics of the Fourier series.

Results show that, independent on the preset wavelength, the macroscopic model can automatically trace the lowest wrinkling critical load. The established nonlinear system for the instability problems of the film/substrate structures shows strong nonlinearity with bifurcation branches, therefore it is difficult to trace the so-called equilibrium path. Two kinds of widely used methods to solve these nonlinear problems can be the classical predictor-corrector methods (e.g. the Newton-Raphson method and the arc-length method) and the perturbation methods (H. S. Shen ( 2013)). By coupling the perturbation method with numerical method, [START_REF] Damil | A new method to compute perturbed bifurcations: Application to the buckling of imperfect elastic structures[END_REF]; [START_REF] Cochelin | Méthode asymptotique numérique[END_REF][START_REF] Cochelin | Asymptotic-numerical methods and Padé approximants for non-linear elastic structures[END_REF] proposed an attractive nonlinear solver known as the Asymptotic Numerical Method (ANM), which is effective and robust especially in tracing bifurcation path in instability problems. In this Chapter, the ANM is used to solve the established nonlinear equations.

This Chapter is structured as follows. The 1D microscopic model for film/substrate system is reviewed in Section 4.2. In Section 4.3, the macroscopic model is deduced by the Fourier series, and the nonlinear equations are solved by the Asymptotic Numerical Method. In Section 4.4, numerical tests have been carefully investigated to valid the established macroscopic model.

Microscopic model

Kinematics

We consider a two-dimensional elastic stiff film resting on an elastic soft compliant substrate, as depicted in Fig. 4.1. Let x and z be the longitudinal and the transverse coordinates, L and b are the length and the width of the structure, respectively. The thickness of the top film and substrate are, respectively, h f and h s . The substrate is ideally divided into a core and a bottom layer of thickness h c and h b as in [START_REF] Yang | A new family of finite elements for wrinkling analysis of thin films on compliant substrates[END_REF]. This division allows to describe the rapid displacement variation along the z direction close to the thin film by high-order kinematic theories and slowly varying kinematics far away from the top film by low-order polynomials, so that less variables are needed to accurately describe the displacement field in the film/substrate structure. As the same kinematic proposed in [START_REF] Yang | A new family of finite elements for wrinkling analysis of thin films on compliant substrates[END_REF], the thin film is modeled as the Euler-Bernoulli's beam:

Film: U f (x, z) = u f 0 (x) -(z - h f + 2h s 2 )W f ,x (x, z) W f (x, z) = w f (x) z ∈ [h s , h s + h f ] (4.1)
where the superscript "f " stands for the film, the notation ", x" stands for ∂ ∂x , U (x, z) and W (x, z) represent the displacements along the x-and z-axis, respectively. The two unknown functions u f 0 and w f are the components of the displacements of the mid-plane of the film. The substrate is modeled as a plane-strain elastic solid: Core:

U c (x, z) = F τ (z)u c τ (x) W c (x, z) = F τ (z)w c τ (x) z ∈ [h b , h s ] τ ∈ [0, n c ] (4.2) Bottom: U b (x, z) = F τ (z)u b τ (x) W b (x, z) = F τ (z)w b τ (x) z ∈ [0, h b ] τ ∈ [0, n b ] (4.3)
where superscripts "c" and "b" stand for the core and bottom layer of substrate, re-spectively. The notation F τ (z) represents the through-the-thickness approximating function within CUF framework, which generally can be an element of a generic approximation base. Within this work, the Mac Laurin's polynomials z n are adopted as expansion function. The notations u τ and w τ are the unknown displacement functions along the beam axis. According to the Einstein rule, repeated indexes denote summation, thus U c (x, z), as an example, can be explicitly written as:

U c (x, z) = u c 0 + zu c 1 + z 2 u c 2 + • • • + z n u c n n ∈ N (4.4)
As a result, a family of one-dimensional refined beam models can be systematically obtained considering the expansion function order as a free parameter of the formulation, that is, n can be an arbitrary value. In this film/substrate system, the variable n is valued as n c and n b for the core and bottom layer of substrate, respectively.

The continuity of the displacement field between the layers film/core and core/bottom is ensured by

U f (x, h s ) = U c (x, h s ) W f (x, h s ) = W c (x, h s ) U c (x, h b ) = U b (x, h b ) W c (x, h b ) = W b (x, h b ) (4.5)
For the interface of the film and the core substrate, the displacement continuity is ensured by substituting Eqs. (4.1) and (4.2) into Eq. (4.5):

u c 0 = u f 0 + 1 2 h f w f ,x -Fi u c i w c 0 = w f -Fi w c i i ∈ [1, n c ] (4.6)
where Fi = h i s . Then Eq. (4.2) can be rewritten as:

U c = u f 0 + 1 2 h f w f ,x + (F i -Fi )u c i W c = w f + (F i -Fi )w c i z ∈ [h b , h s ] i ∈ [1, n c ] (4.7)
The congruency of the displacement field at the interface

Γ = {(x, y, z) : x ∈ [0, L], y ∈ [-b/2, b/2], z = h b } of
L (µ, u) = ∫ Γ µ 1 [ U c (x, h b ) -U b (x, h b ) ] + µ 2 [ W c (x, h b ) -W b (x, h b ) ] dΓ (4.8)

Geometric equation and constitutive law

The geometric equation and constitute law meet the following assumptions:

1. the material behavior is linearly elastic and meets the Hooke's law, 2. the Green-Lagrange strain is used to describe the relationships between strain and displacement in the film, 3. the small displacement with moderate rotation is assumed for the film and the substrate small displacement.

The above assumptions are given by the following equations:

σ f xx = E f ϵ f xx = E f (U f ,x + 1 2 (W f ,x ) 2 ) σ s xx = C s 11 ϵ s xx + C s 13 ϵ s zz = C s 11 U s ,x + C s 13 W s ,z σ s zz = C s 13 ϵ s xx + C s 33 ϵ s zz = C s 13 U s ,x + C s 33 W s ,z σ s xz = 2C s 55 ϵ s xz = C s 55 ( U s ,z + W s ,x ) (4.9) with C s 11 = C s 33 = 1 -ν s (1 + ν s )(1 -2ν s ) E s C s 13 = ν s (1 + ν s )(1 -2ν s ) E s C s 55 = E s 2(1 + ν s ) (4.10)
where E f is the Young's modulus of the film, E s and ν s are the Young's modulus and Poisson's ratio of the substrate. The subscript "s" represents "c" or "b" for the core or bottom layer of substrate for brevity.

Considering the constrained variational problem, the weak form of the governing equations are obtained by using the Principal of Virtual Displacement as follows:

δP int + δL = δP ext (4.11)
where δP int is the virtual internal work and δP ext the virtual external work. The internal virtual work of the film/substrate system is expressed as:

δP int = ∫ V f σ f xx δϵ f xx dV f + ∫ V c (σ c xx δϵ c xx + σ c zz δϵ c zz + 2σ c xz δϵ c xz )dV c + ∫ V b (σ b xx δϵ b xx + σ b zz δϵ b zz + 2σ b xz δϵ b xz )dV b (4.12)
where V f , V c and V b express the integration volumes of the film, core and bottom substrates, respectively. By using Eqs. (4.1-4.3) and substituting Eqs. (4.6-4.10) into Eq. (4.12), a 1D microscopical model could be deduced, see more details in [START_REF] Yang | A new family of finite elements for wrinkling analysis of thin films on compliant substrates[END_REF]. In the next part, the technique of slowly varying Fourier coefficient is applied to the above equations to derive the one-dimensional macroscopical model.

Macroscopic model

In this part, the macroscopic film/substrate model will be presented based on the concept of Fourier series with slowly varying coefficients. The principle of the Fourier-related approach is to write the unknown field in the following form:

U (x, z) = +∞ ∑ j=-∞ U j (x, z)e ijqx (4.13)
where the Fourier coefficient U j (x, z) denotes the envelope for the jth order harmonic, which is conjugated with U -j (x, z). For simplicity, the macroscopic unknown field U j (x, z) is assumed to slowly vary in a single direction over a period

[ x, x + 2π q
] of the oscillation. The wavenumber q is defined as

q = π L q 0 (4.14)
in which, L is the length of film/substrate and q 0 is a chosen integer wave-parameter: q 0 ∈ N and q 0 ≥ 2. In practice, the number of Fourier coefficients considered is finite. But at least two functions U 0 (x, z) and U 1 (x, z) are needed to describe nearly periodic patterns: U 0 (x, z) can be identified to the mean value while U 1 (x, z)

represents the envelope or the amplitude of the spatial oscillations, see Fig.

General macroscopic formulation

To establish the macroscopic equations satisfied by the envelopes U j , we firstly introduce some rules proposed by [START_REF] Damil | Influence of local wrinkling on membrane behaviour: a new approach by the technique of slowly variable Fourier coefficients[END_REF]. If a(x) and b(x)

are Fourier series with slowly varying coefficients as in Eq. (4.13), the following identities hold:

∫ l 0 a(x)b(x)dx = ∫ l 0 ∞ ∑ j=-∞ a j (x)b -j (x)dx (4.15) ( da dx ) j = (a ′ ) j = ( d dx + ijq ) a j = (a j ) ′ + ijqa j (4.16) ( d 2 a dx 2 ) j = (a ′′ ) j = ( d dx + ijq
) 2 a j = (a j ) ′′ + 2ijq(a j ) ′ -j 2 q 2 a j (4.17)

(ab) j = +∞ ∑ j 1 =-∞ a j 1 b j-j 1 (4.18)
Using the above identities, the macroscopical version of the displacement fields in Eqs. (4.1-4.3) are expressed as

U f j (x, z) = (u f 0 ) j -(z - h f + 2h s 2 ) ( d dx + ijq ) W f j W f j (x, z) = w f j z ∈ [h s , h s + h f ] (4.19) U c j (x, z) = F τ (u c τ ) j W c j (x, z) = F τ (w c τ ) j z ∈ [h b , h s ] τ ∈ [0, n c ] (4.20) U b j (x, z) = F τ ( u b τ ) j W b j (x, z) = F τ ( w b τ ) j z ∈ [0, h b ] τ ∈ [0, n b ] (4.21)
and the macroscopic form of Eq. (4.7) could be written as

U c j = (u f 0 ) j + 1 2 h f ( d dx + ijq)w f j + (F i -Fi )(u c i ) j W c j = w f j + (F i -Fi )(w c i ) j (4.22)
The macroscopic constitutive laws are thus as follows:

( σ f xx ) j = E f ( ϵ f xx ) j = E f   ( d dx + ijq)U f j + 1 2 +∞ ∑ j 1 =-∞ ( d dx + ij 1 q)( d dx + i(j -j 1 )q)W f j 1 W f j-j 1   (σ s xx ) j = (C 11 ϵ s xx + C 13 ϵ s zz ) j = C 11 ( d dx + ijq)U s j + C 13 W s j,z (σ s zz ) j = (C 13 ϵ s xx + C 33 ϵ s zz ) j = C 13 ( d dx + ijq)U s j + C 33 W s j,z (σ s xz ) j = 2C 55 (ϵ s xz ) j = C 55 ( U s j,z + ( d dx + ijq)W s j ) (4.23)
After substituting Eq. (4.23) into Eq. (4.12), the macroscopic version of the internal virtual work is:

δP int = ∫ V f +n j ∑ j=-n j (σ f xx ) j δ(ϵ f xx ) -j dV f ∫ V c +n j ∑ j=-n j ( (σ c xx ) j δ(ϵ c xx ) -j + (σ c zz ) j δ(ϵ c zz ) -j + 2(σ c xz ) j δ(ϵ c xz ) -j ) dV c ∫ V b +n j ∑ j=-n j ( (σ b xx ) j δ(ϵ b xx ) -j + (σ b zz ) j δ(ϵ b zz ) -j + 2(σ b xz ) j δ(ϵ b xz ) -j ) dV b (4.24)
where n j is the number of Fourier coefficients. In the following, before the finite element computational scheme is established, some assumptions are introduced according to [START_REF] Damil | Influence of local wrinkling on membrane behaviour: a new approach by the technique of slowly variable Fourier coefficients[END_REF]:

u f 0 , u c τ , u b τ do not fluctuate, i.e u f 0 def = (u f 0 ) 0 , u c τ def = (u c τ ) 0 , u b τ def = (u b τ ) 0 .

Internal virtual work of the film

As the thin film is modeled as the nonlinear beam, the virtual work of the film could be expressed as:

δP f int = ∫ V f σ f xx δϵ f xx dV f = ∫ L 0 N f δγ + M f δw ,xx dx (4.25) in which γ = u f 0,x + 1 2 w f ,x 2 N f = E f bh f γ M f = E f I f w f ,xx I f = 1/12bh 3 f (4.26)
The macroscopic version of internal virtual work of the film in Eq. (4.25) could be rewritten as:

δP f int = ∫ V f +n j ∑ j=-n j (σ f xx ) j δ(ϵ f xx ) -j dV f = ∫ L 0 +n j ∑ j=-n j (N f j δγ -j + M f j δ(w f ,xx ) -j ) dx (4.27)
With the assumption that the axial force N f (x) does not fluctuate [START_REF] Damil | Influence of local wrinkling on membrane behaviour: a new approach by the technique of slowly variable Fourier coefficients[END_REF]:

N f (x) = N f 0 (x), γ(x) = γ 0 (x), Eq.
(5.42) can be simplified as:

δP f int = ∫ L 0 N f 0 δγ 0 + +n j ∑ j=-n j M f j δ(w f ,xx ) -j dx (4.28)
where

N f 0 = E f bh f ( u f 0,x + 1 2 +n j ∑ j=-n j ( d dx + ijq)( d dx -ijq)w f j w f -j ) M f j = E f I f ( ( d 2 dx 2 + 2ijq d dx -j 2 q 2 )w f j ) (4.29)
The generalized strain {γ f } of the film could be defined as:

{γ f } = ([H] + 1 2 [A(θ f )]){θ f } (4.30)
in which

[H] =                     1 0 0 0 0 0 • • • 0 0 0 0 0 1 0 0 0 • • • 0 0 0 0 0 0 1 0 0 • • • 0 0 0 0 0 0 0 1 0 • • • 0 0 0 0 0 0 0 0 1 • • • 0 0 0 . . . . . . . . . . . . . . . . . . . . . 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1                     [ A ( θ f )] = 2                     0 1 2 w f 0,x 0 q 2 w f 1 w f 1,x 0 • • • j 2 q 2 w f n j w f n j ,x 0 0 0 0 0 0 0 • • • 0 0 0 0 0 0 0 0 0 • • • 0 0 0 0 0 0 0 0 0 • • • 0 0 0 0 0 0 0 0 0 • • • 0 0 0 . . . . . . . . . . . . . . . . . . . . . 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0                     (4.31) ⟨ θ f ⟩ = ⟨ u f 0,x w f 0,x w f 0,xx w f 1 w f 1,x w f 1,xx • • • w f n j w f n j ,x w f n j ,xx ⟩ (4.32)
where

[H] ∈ R (2+3n j )×(3+3n j ) , [ A ( θ f )] ∈ R (2+3n j )×(3+3n j ) and {θ f } ∈ R (3+3n j )×1 . Note that [ A ( θ f )]
and {θ f } are linear functions of u f 0 , w f 0 and w f j , thus the internal virtual work of the film Eq. (4.28) is a quadratic expression in terms of the displacement, and the generalized stress { S f } is defined as:

{ S f } = [D f ]([H] + 1 2 [A(θ f )]) { θ f } (4.33)
in which:

[ D f ] = diag ( E f bh f 2E f bh f E f I f 2q 4 E f I f 12q 2 E f I f 2E f I f • • • 2n 4 j q 4 E f I f 12n 2 j q 2 E f I f 2E f I f ) (4.34)
To unify the unknowns, the unknown vector {θ f } of the film is expressed as a function of the whole unknown vector {θ} ∈ R (

6+4n j +(4+2n j )(nc+n b +1) ) ×1 by introducing a transform matrix [R f ]: { θ f } = [R f ] {θ} (4.35)
where

⟨θ⟩ = ⟨ u f 0 u f 0,x w f 0 w f 0,x w f 0,xx • • • w f n j w f n j ,x w f n j ,xx u c 1 u c 1,x • • • u c nc u c nc,x (w c 1 ) 0 (w c 1,x ) 0 • • • (w c nc ) 0 (w c nc,x ) 0 • • • (w c 1 ) n j (w c 1,x ) n j • • • (w c nc ) n j (w c nc,x ) n j u b 0 u b 0,x • • • u b n b u b n b ,x (w b 0 ) 0 (w b 0,x ) 0 • • • (w b n b ) 0 (w b n b ,x ) 0 • • • (w b 0 ) n j (w b 0,x ) n j • • • (w b n b ) n j (w b n b ,x ) n j µ 1 (µ 2 ) 0 • • • (µ 2 ) n j ⟩ (4.36)
for the sake of brevity, the details of the transform matrix are not presented. Thus, the virtual work of the film could be written as:

δP f int = ∫ L 0 ⟨δγ f ⟩{S f }dx = ∫ L 0 ⟨δθ⟩ T [R f ] T ([H] + [A(θ f )]){S f }dx (4.37)
where T [2] represents the transpose of a matrix.

Internal virtual work of the core substrate

By applying the macroscopic displacement fields Eqs. (4.20), (4.22) and the macroscopic stress Eq. (4.23) of core substrate into Eq. (4.24), the virtual work of the core substrate becomes: model. The mechanism of the phenomenon is carefully explained and two proper remedies are proposed to treat it in Huang et al. (2017a). In this Chapter, one of the two remedies, similar to the way of treating the Poisson locking in plate theories (Carrera and Brischetto (2008a)), is used by modifying the constitutive relations in the substrate. As a result, the modified elastic coefficients are

δP c int = ∫ V c +n j ∑ j=-n j ( (σ c xx ) j δ(ϵ c xx ) -j + (σ c zz ) j δ(ϵ c zz ) -j + 2(σ c xz ) j δ(ϵ c xz ) -j ) dV c = ∫ L 0 [ N c xx0 δu f 0,x + M c 1xz0 δw f 0,x + M c 2xx0 δw f 0,xx + M c xzxxj δw f -j + M c 1xzxxj δw f -j,x + M c 2xxj δw f -j,xx + N c xzi δu c i + N c xxi δu c i,x + M c zzi0 δ(w c i ) 0 + M c 1xzi0 δ(w c i,x ) 0 + M c zzxzij δ(w c i ) -j + M c 1xzij δ(w c i,x ) -j ] (4.38) where i ∈ [1, n c ], j ∈ [1, n j ]. Because U s (x,
C 11 = C 33 = E s , C 13 = 0, C 55 = E s /2(1 + ν s ).
The core stress resultants are:

N c xx0 b = J 11c 00 (u f 0,x + h f 2 w f 0,xx ) + J 11c 0i u c i,x + J 13c 0i,z (w c i ) 0 M c 1xz0 b = J 55c 00 w f 0,x + J 55c 0i,z u c i + J 55c 0i w c 0,x M c 2xx0 b = J 11c 00 h f 2 (u f 0,x + h f 2 w f 0,xx ) + J 11c 0i h f 2 u c i,x + J 13c 0i h f 2 (w c i ) 0 M c xzxxj b = (j 2 q 2 J 11c jj + J 55c jj )w f j -J 11c jj w f j,xx + J 55c ij (w c i ) j -J 13c 0i,z h f j 2 q 2 (w c i ) j M c 1xzxxj b = J 1155c 1jj w f j,x + J 55c 1ij (w c i,x ) j M c 2xxj b = -J 11c jj w j + J 11c 00 h 2 f 2 w f j,xx + J 13c 0i,z h f (w c i ) j N c xzi b = J 55c i,z0 w f 0,x + J 55c i,zm,z u c m + J 55c i,zm (w c m,x ) 0 N c xxi b = J 11c i0 (u f 0,x + h f 2 w f 0,xx ) + J 11c im u c m,x + J 13c im,z (w c m ) 0 M c zzi0 b = J 33c i,zm,z (w c m ) 0 + J 13c i,z0 u f 0,x + J 13c i,z0 h f 2 j 2 q 2 w f 0,xx + J 13c i,zm u c i,x M c 1xzi0 b = J 55c i0 w f 0,x + J 55c im,z u c m + J 55c im (w c m,x ) 0 M c zzxzij b = 2j 2 q 2 J 55c i0 w f j + (2J 33c i,zm,z + 2j 2 q 2 J 55c im )(w c m ) j +2J 13c i,z0 ( h f 2 w f j,xx -j 2 q 2 w f j ) M c 1xzij b = 2J 55c 0i w f j,x + 2J 55c im (w c m,x ) j (4.39)
where J 2c 2 are the constant terms related to the elastic coefficients and the thickness of the core substrate

J xxc im = ∫ hs h b C c xx (F i -Fi )(F m -Fm )dz J xxc i,zm,z = ∫ hs h b C c xx (F i -Fi ) ,z (F m -Fm ) ,z dz J xxc i,zm = ∫ hs h b C c xx (F i -Fi ) ,z (F m -Fm )dz J xxc im,z = ∫ hs h b C c xx (F i -Fi )(F m -Fm ) ,z dz J xxc 00 = C c xx h c J xxc i0 = J xxc 0i = ∫ hs h b C c xx (F i -Fi )dz J xxc i,z0 = J xxc 0i,z = ∫ hs h b C c xx F i,z dz J 11c jj = C c 11 h 2 f 2 j 2 q 2 h c , J 55c jj = 2C c 55 j 2 q 2 h c , J 55c ij = 2j 2 q 2 J 55c 0i J 1155c 1jj = 4J 11c jj + 2J 55c 00 , J 55c 1ij = 2J 55c 0i (4.40)
The superscript "xx" represents 11, 13, 33 or 55, see details in [START_REF] Yang | A new family of finite elements for wrinkling analysis of thin films on compliant substrates[END_REF].

As a consequence, Eqs. (4.38) and (4.39) can be rewritten in the matrix form in terms of displacements:

δP c int = ∫ L 0 ⟨δθ c ⟩ [D c ] {θ c }dx (4.41)
where [D c ] is the generalised stiffness matrix of core layer that is not here explicitly reported for the sake of brevity, {θ c } ∈ R ( 3+3n j +4nc+2ncn j ) ×1 is a vector whose components are the core generalised displacement unknown functions and their derivatives versus the axial coordinate. It reads as:

⟨θ c ⟩ = ⟨u f 0,x w f 0,x w f 0,xx w f 1 w f 1,x w f 1,xx • • • w f n j w f n j ,x w f n j ,xx (u c 1 ) 0 • • • (u c nc ) 0 (u c 1,x ) 0 • • • (u c nc,x ) 0 (w c 1 ) 0 • • • (w c nc ) 0 (w c 1,x ) 0 • • • (w c nc,x ) 0 • • • (w c 1 ) n j • • • (w c nc ) n j (w c 1,x ) n j • • • (w c nc,x ) n j ⟩ (4.42)
By unifying the unknowns, the internal virtual work of the core substrate is formulated as:

δP c int = ∫ L 0 ⟨δθ c ⟩ [D c ] {θ c }dx = ∫ L 0 ⟨δθ⟩ T [R c ] [D c ] [R c ] {θ}dx = ∫ L 0 ⟨δθ⟩ [K c ] {θ}dx (4.43)
where [R c ] is the transition matrix between the total unknown vector{θ} and the unknown vector of the core {θ c }.

Internal virtual work of the bottom substrate

The derivation process of internal virtual work of the bottom substrate is more or less the same as that of the core layer by replacing Eq. ( 4.21) and integrating over through-the-thickness and -width coordinates, the internal virtual work of the bottom layer is given as following:

δP b int = ∫ V b +n j ∑ j=-n j ( (σ b xx ) j δ(ϵ b xx ) -j + (σ b zz ) j δ(ϵ b zz ) -j + 2(σ b xz ) j δ(ϵ b xz ) -j ) dV b = ∫ L 0 [ N b xzn δu b n + N b xxn δu b n,x + M b zzn0 δ(w b n ) 0 + M b 1xzn0 δ(w b n,x ) 0 + M b xznj δ(w b n ) -j + M b 1xznj δ(w b n,x ) -j ] dx (4.44) where n ∈ [0, n b ], j ∈ [1, n j ].
The bottom stress resultants are:

N b xzn b = J 55b n,zp,z u b p + J 55b n,zp (w b p,x ) 0 N b xxn b = J 11b np u b p,x + J 13b np,z (w b p ) 0 M b zzn0 b = J 33b n,zp,z (w b p ) 0 + J 13b n,zp u b p,x M b 1xzn0 b = J 55b np,z u b p + J 55b np (w b p,x ) 0 M b xznj b = (2J 33b n,zp,z + 2j 2 q 2 J 55b np )(w b p ) j M b 1xznj b = 2J 55b np (w b p,x ) j (4.45)
where

J xxb np = ∫ h b 0 C b xx F n F p dz J xxb n,zp,z = ∫ h b 0 C b xx F n,z F p,z dz J xxb n,zp = ∫ h b 0 C b xx F n,z F p dz J xxb np,z = ∫ h b 0 C b xx F n F p,z dz (4.46)
By introducing a transform matrix, the internal virtual work of the bottom substrate could be expressed as:

δP b int = ∫ L 0 ⟨δθ b ⟩ [ D b ] {θ b }dx = ∫ L 0 ⟨δθ⟩ T [ R b ] [ D b ] [ R b ] {θ}dx = ∫ L 0 ⟨δθ⟩ [ K b ] {θ}dx (4.47)
in which

⟨θ b ⟩ = ⟨(u b 0 ) 0 • • • (u b n b ) 0 (u b 0,x ) 0 • • • (u b n b ,x ) 0 (w b 0 ) 0 • • • (w b n b ) 0 (w b 0,x ) 0 • • • (w b n b ,x ) 0 • • • (w b 0 ) n j • • • (w b n b ) n j (w b 0,x ) n j • • • (w b n b ,x ) n j ⟩ (4.48)
where

{θ b } ∈ R 2(n b +1)(n j +2)×1 , [R b ] ∈ R 2(n b +1)(n j +2)× ( 6+4n j +(4+2n j )(nc+n b +1)
) .

The congruency of the displacement field at core and bottom layers' interface Γ is realised by means of the Lagrange multiplier µ through δL , whose variation formulation is as following:

δL =L (δµ, u) + L (µ, δu) = ∫ Γ δµ 1 [ U c (x, h b ) -U b (x, h b ) ] + δ(µ 2 ) 0 [ W c 0 (x, h b ) -W b 0 (x, h b ) ] + δ(µ 2 ) j [ W c j (x, h b ) -W b j (x, h b ) ] dΓ + ∫ Γ µ 1 [ δU c (x, h b ) -δU b (x, h b ) ] + (µ 2 ) 0 [ δW c 0 (x, h b ) -δW b 0 (x, h b ) ] + (µ 2 ) j [ δW c j (x, h b ) -δW b j (x, h b )
] dΓ

δL = b ∫ L 0 δµ 1 [ u f 0 + 1 2 h f n j ∑ j=-n j ( d dx + ijq)w j + ( F i -Fi )u c i -F n u b n ] + δ(µ 2 ) 0 [ w f 0 + ( F i -Fi )(w c i ) 0 -F n (w b n ) 0 ] + δ(µ 2 ) j [ w f j + ( F i -Fi )(w c i ) j -F n (w b n ) j ] dx + b ∫ L 0 µ 1 [ δu f 0 + 1 2 h f (δw f 0,x + n j ∑ j=-n j ( d dx + ijq)δw j ) + ( F i -Fi )δu c i -F n δu b n ] + (µ 2 ) 0 [ δw f 0 + ( F i -Fi )δ(w c i ) 0 -F n δ(w b n ) 0 ] + (µ 2 ) j [ δw f j + ( F i -Fi )δ(w c i ) j -F n δ(w b n ) j ] dx (4.50) in which F τ = h τ b .
After introducing a transform matrix [R cp ], Eq. (4.50) can be rewritten into the matrix form:

δL (µ, u) = ∫ L 0 ⟨δθ cp ⟩[D cp ]{θ cp }dx = ∫ L 0 ⟨δθ⟩ T [R cp ][D cp ][R cp ]{θ}dx = ∫ L 0 ⟨δθ⟩[K cp ]{θ}dx (4.51)

Governing equation

In this Chapter, the body forces are neglected, then the external virtual work could be

δP ext = λ ∫ Sσ ⟨δθ⟩{F}dS (4.52)
where {F} is the external applied force, and λ is a "load parameter". Hence, the macroscopic governing equation is expressed as:

∫ L 0 ⟨δθ⟩ T [R f ] T ([H] + [A(θ f )]){S f }dx + ∫ L 0 ⟨δθ⟩([K c ] + [ K b ] + [K cp ]){θ}dx = λ ∫ Sσ ⟨δθ⟩{F}dS { S f } = [D f ]([H] + 1 2 [A(θ f )])[R f ] {θ} (4.53)

Discretization

The finite element method is used to solve the nonlinear system of Eq. (4.53). A two-node 1D element with (3 + 2n j ) • (n c + n b + 3) degrees of freedom for each node is used, which results in the following nodal unknown vector:

q IT = { u f 0 w f 0 w f 0,x w f 1 w f 1,x • • • w f n j w f n j ,x u c 1 • • • u c nc (w c 1 ) 0 • • • (w c nc ) 0 (w c 1,x ) 0 • • • (w c nc,x ) 0 • • • (w c 1 ) n j • • • (w c nc ) n j (w c 1,x ) n j • • • (w c nc,x ) n j u b 0 • • • u b n b (w b 0 ) 0 • • • (w b n b ) 0 (w b 0,x ) 0 • • • (w b n b ,x ) 0 • • • (w b 0 ) n j • • • (w b n b ) n j (w b 0,x ) n j • • • (w b n b ,x ) n j µ 1 (µ 2 ) 0 (µ 2,x ) 0 • • • (µ 2 ) n j (µ 2,x ) n j } (4.54) {v} = { v I v II } (4.55)
where I and II denote the first and the second node of the element, {v} is the elementary vector of unknowns. The longitudinal unknown displacement functions u f 0 , u c τ , u b τ and the Lagrange multiplier µ 1 are discretised by the Lagrange linear shape functions, while the unknowns transverse displacement functions w f j , (w c τ ) j , (w b τ ) j and the Lagrange multipliers (µ 2 ) j by the Hermite shape functions. The unified discretization of the unknowns can be obtained by:

V L = Γ(1)V I L + Γ(2)V II L = ⟨N L ⟩ {v} V H = Γ(3)V I H + Γ(4)V I H + Γ(5)V II H + Γ(6)V II H = ⟨N H ⟩ {v} (4.56)
in which, V L denotes the linear functions and V H denotes the Hermite functions, ⟨N L ⟩ and ⟨N H ⟩ are corresponding interpolation function matrices whose nonzero elements are determined by:

Γ = ⟨ 1 -ξ 2 1 + ξ 2 (1 -ξ) 2 (2 + ξ) 4 l e (1 -ξ 2 )(1 -ξ) 8 (1 + ξ) 2 (2 -ξ) 4 l e (-1 + ξ 2 )(1 + ξ) 8 ⟩ (4.57)
The parameter ξ is the elementary local coordinate and l e is the length of the element. Here, we take u c i and (w c i ) j for example:

u c i = Γ(1)u cI i + Γ(2)u cII i = ⟨ N c u i ⟩ {v} (4.58) (w c i ) j = Γ(3)(w c i ) I j + Γ(4)(w c i,x ) I j + Γ(5)(w c i ) II j + Γ(6)(w c i,x ) II j = ⟨ N c w ij ⟩ {v} (4.59) where ⟨ N c u i ⟩ and ⟨ N c w ij
⟩ are, respectively, the interpolation function matrix for u c i and (w c i ) j . Other unknowns can be discretised as similar to these definitions. Thus, the unknown vector {θ} is expressed as:

{θ} = [G]{v} (4.60)
where [G] is assembled by ⟨N L ⟩ , ⟨N H ⟩ and their derivatives. By substituting Eq.

(4.60) in Eq. ( 4.53), the discretization form of the governing equation is written as:

∑ e ⟨δv⟩ ∫ le 0 ( T [B f (v)]{S f } + ( [ K c ] + [ K b ] + [ K cp ] ){v} ) dx = λ ∑ e ⟨δv⟩{f } e {S f } = [D f ]([B f l ] + 1 2 [B f nl (v)]){v} (4.61)
where

[B f (v)] = [B f l ] + [B f nl (v)] [B f l ] = [H][R f ][G] [B f nl (v)] = [A(v)][R f ][G] [ K c ] = T [G] [K c ] [G] [ K b ] = T [G] [ K b ] [G] [ K cp ] = T [G] [K cp ] [G] {f} e = ∫ Sσ T [G]{F}dS (4.62)
The Asymptotic Numerical Method (ANM) [START_REF] Cochelin | Asymptotic-numerical methods and Padé approximants for non-linear elastic structures[END_REF]) is used to solve the set of Eq. (4.61). For each element the unknown vector {v} is determined at the nodes and the stress vector {S f } is determined at the Gauss points. Details on the use of the ANM for the resolution of non linear problems are given in Appendix A 4.5. Besides, the various boundary conditions of macroscopic model have been deduced in [START_REF] Liu | A new Fourierrelated double scale analysis for instability phenomena in sandwich structures[END_REF], which will not be illustrated any more in this Chapter.

Numerical results

In this part, the film/substrate system is subjected to compression as illustrated in Fig. 4 The results obtained by the macroscopic model are compared with the numerical results obtained by a two-dimensional nonlinear elasticity solution with the finite element code ABAQUS, which is referred to as "2D full model" in this paper. The film is modeled by 3-node beam elements (B22) and the substrate by 8-node planestrain quadrilateral reduced integration elements (CPE8R). The rotational degree of freedom in the beam elements is constrained when connected with the planestrain elements. To ensure accuracy, very fine meshes are adopted to discrete the film and substrate.

Validation

Table 4.1 presents the material properties and geometric parameters. The Young's modulus ratio E s /E f ranges from 10 -5 to 10 -3 . The thickness ratio h s /h f changes between 10, for very thin substrates, and 10 3 , for very thick substrates. This case is referred to as "problem I" to validate the 1D macroscopic model. The chosen expansion orders n c and n b have been proved to be sufficient enough for the considered cases, see [START_REF] Yang | A new family of finite elements for wrinkling analysis of thin films on compliant substrates[END_REF]. For the present, 10 elements are used along the axis (mesh convergence will be investigated later), which results in about 4.29 • 10 2 DOFs. It is a half less than the minimum DOFs 10 3 by the 1D microscopic model in [START_REF] Yang | A new family of finite elements for wrinkling analysis of thin films on compliant substrates[END_REF]. Besides, the 1D microscopic model's DOFs increase rapidly as the wavelength becomes shorter, whereas the mesh size of the 1D macroscopic model is independent on the wavelength change, which makes the 1D macroscopic model even more attractive considering that accurate results are obtained.

E f (10 4 MPa) E s /E f ν s L(10 -1 m) b(m) h f /L h s /h f
Table 4.2: Half-wave number q 0 and critical load parameter λ c . The 1D macroscopic solution via a n c = 7, n b = 3 and n j = 1 model with 10 axial elements. "2D" is short for the 2D full model by ABAQUS. Problem I.

h s /h f 10 10 2 10 3 q 0 10 -4 λ c q 0 10 -4 λ c q 0 10 -4 λ c E s /E f 1D 2D 1D 2D 1D 2D

Mesh convergence

In this section, the loading and boundary conditions of the film/substrate remain the same as shown in Fig. 4.3, and the material properties and the geometrical data are presented in Table 4.3, which has been studied by the 1D microscopic model in [START_REF] Yang | A new family of finite elements for wrinkling analysis of thin films on compliant substrates[END_REF]. This case is referred to as "problem II".

Here, the thickness of core layer h c = 2ℓ w (ℓ w is the wrinkling wavelength defined as ℓ w = L/q 0 ) proposed in Yang et al. ( 2015) is adopted in the 1D macroscopic model. The expansion orders for the core layer and bottom layer of the substrate are Furthermore, the above critical load corresponds to the wrinkling pattern with the lowest critical load for the considered geometry and material configuration in problem II. The computation expenses will increase for buckling modes with more wavenumbers and smaller wavelengths for the 2D full model or the 1D microscopic model [START_REF] Yang | A new family of finite elements for wrinkling analysis of thin films on compliant substrates[END_REF]), since a larger number of elements are needed to accurately simulate this phenomenon. Sometimes special imperfection should be applied to get the higher mode, which may lead to a great trouble in controlling the nonlinear solution path, see Hu et al. (2009a). However, the 1D macroscopic model easily overcomes the above shortcomings by just changing the wavenumber q 0 , and the corresponding critical load can be detected accurately and fast. whereas the 1D macroscopic model only uses 2 elements (99 DOFs) to describe the higher buckling mode by applying a transverse perturbation force with a value of F/10 5 on the center of the system. The bifurcation curve in Fig. 4.7 shows that the critical load of the macroscopic model is quite accurate with such a limited degrees of freedom. Therefore, the macroscopic model easily overcomes the disadvantages of the 1D microscopic model or the 2D full model and significantly improves the computational efficiency.

n c =
Table 4.3: Material and geometric parameters of the film/substrate. Problem II.

E f (10 5 MPa) E s /E f ν s L(10 -3 m) b(10 -3 m) h f /L h s /h f 1.8 10 -5 0.3 1.0 1.0 10 -3 10 3

Fourier coefficient study

In this section, several Fourier coefficients U j (x, y) (referred specially to the Fourier coefficients of transverse displacement W in the paper) are used to simulate the instability phenomenon. Usually, only the zero harmonic U 0 (mean field) and the first harmonic U 1 (envelop of the oscillation) are adopted, see Damil and Potier-Ferry In this way, we can easily distinguish the harmonic that contributes the most to the bifurcation curves and conveniently predict the lowest wrinkling pattern with a random preset wavenumber. For simplicity, we take a shorter film/substrate system with L = 0.8 • 10 -3 m, and the other parameters remain the same as problem II. This case is referred to as 'problem III'.

For this considered case, the lowest wrinkling pattern has 8 wavenumbers as depicted in Fig. 4.8. The Fourier coefficients U 1 , U 2 , U 3 , U 4 and U 5 are used in the 1D macroscopic model. These coefficients are assumed to be real, that is to say, the amplitude modulation is taken into account and the evolution of the phase is disregarded. The mean field U 0 is quite small and negligible in the considered film/substrate instability phenomenon. Small transverse perturbation forces F/10 5 are applied at the center of the structure. Results have been computed via a model with n c = 5 and n b = 3 and 2 elements along the x axis. First, the wavenumber q 0 is chosen as q 0 = 2. Fig. 4.9 presents the bifurcation curves obtained by the specified 1D macroscopic model. The notation U represents the microscopic transverse displacement, whereas U j is the Fourier coefficient describing the jth

= 3, n j = 1, h c = 2ℓ w ). Problem II.
envelop of the transverse displacement. One can see that the coefficient of the 4th harmonic dominates the bifurcation path and leads to almost the same critical load predicted by the 2D full model. It seems that the 1D macroscopic model could automatically choose the proper Fourier coefficient to ensure the system in a state with minimum energy. Then the wavenumber q is changed to q 0 = 4, and Fig. 4.10 shows that the 1D macroscopic model choose the 2th harmonic as the main fluctuation mode. Besides, the out-of-plane deformation of the film near bifurcation point for the cases q 0 = 2 and q 0 = 4 has been plotted as shown in Fig. 4.11, which indicates that the corresponding main fluctuation mode describes the instability pattern correctly. There in fact exists fluctuation in the transversal displacement for other Fourier coefficients, but they are too small and makes little contribution to the buckling of the system. The 1D macroscopic model accounting for several Fourier coefficients could conveniently and correctly predict the lowest critical load in the film/substrate instability phenomenon. Compared with the Landau-Ginzburg equation deduced from the asymptotic approach, the proposed model has an advantage to couple the evolution of several harmonics or envelops of the microscopic displacement field.

Conclusion

Based on a high order reduced-dimension model proposed by [START_REF] Yang | A new family of finite elements for wrinkling analysis of thin films on compliant substrates[END_REF], where the Euler-Bernoulli's beam was used for the thin film and the Carrera's Unified Formulation(CUF) for the substrate, a new macroscopic model is developed by the technique of slowly varying Fourier coefficient. The computational efficiency is significantly improved compared with the former model, especially when dealing with wrinkling instability phenomena with vast wavenumbers. Besides, several harmonics in Fourier series have been taken into account, and numerical experiments show that the proposed macroscopic model automatically predicts the lowest buckling mode. The present model is based on the method of Fourier series with variable coefficients that has been widely used for several structural models, such as the beam on Winkler foundation [START_REF] Damil | Influence of local wrinkling on membrane behaviour: a new approach by the technique of slowly variable Fourier coefficients[END_REF]; [START_REF] Hu | A bridging technique to analyze the influence of boundary conditions on instability patterns[END_REF]), the sandwich models [START_REF] Liu | A new Fourierrelated double scale analysis for instability phenomena in sandwich structures[END_REF]; [START_REF] Yu | Multi-scale techniques to analyze instabilities in sandwich structures[END_REF]) and the membrane model [START_REF] Damil | Membrane wrinkling revisited from a multiscale point of view[END_REF], [START_REF] Huang | Macroscopic simulation of membrane wrinkling for various loading cases[END_REF]). However, sometimes the results near the boundary show poor accuracy, which can be overcome by introducing the Arlequin method [START_REF] Dhia | Global-local approaches: the Arlequin framework[END_REF], [START_REF] Hu | A bridging technique to analyze the influence of boundary conditions on instability patterns[END_REF][START_REF] Hu | Multi-scale nonlinear modelling of sandwich structures using the Arlequin method[END_REF], [START_REF] Yu | Multi-scale techniques to analyze instabilities in sandwich structures[END_REF]) coupling the microscopic model on the boundary and macroscopic model in the bulk.

Appendix A. Asymptotic numerical method

The solution to many physics problems can be achieved through resolution of nonlinear problems depending on real parameter Θ. The corresponding nonlinear system of equations can be written as follow:

R (Θ; ϑ) = 0 (4.63)
Assuming that a solution of Eq. 4.63 at point j is known, the solution (Θ; ϑ) j+1 at point j + 1 is obtained from this latter following a perturbation technique. An approached solution path at step j + 1 is represented by a n-order Taylor series expansion of Θ and ϑ:

Θ = Θ j + aΘ 1 + a 2 Θ 2 + . . . + a n Θ n = Θ j + a p Θ p ϑ = ϑ j + aϑ 1 + a 2 ϑ 2 + . . . + a n ϑ n = ϑ j + a p ϑ p (4.64)
where a ∈ R is a solution path parameter for the non-linear problem. According to Einstein's notation, p is a dummy index that implicitly stands for summation over the range [1, n]. By replacing Eqs. 4.64 within the non-linear problem, Eq. 4.63 reads:

R ( Θ j + a p Θ p ; ϑ j + a p ϑ p ) = 0 (4.65)
Its n-order Taylor series expansion about the point (Θ; ϑ) j is, then:

R ( Θ j ; ϑ j ) + ∂R ∂Θ j a p Θ p + ∂R ∂ϑ j a p ϑ p + + 1 2 [ ∂ 2 R ∂Θ 2 j a p Θ p a s Θ s + ∂ 2 R ∂ϑ 2 j a p ϑ p a s ϑ s + ∂ 2 R ∂Θ∂ϑ j a p Θ p a s ϑ s ] + + 1 3! [ ∂ 3 R ∂Θ 3 j a p Θ p a s Θ s a q Θ q + ∂ 3 R ∂ϑ 3 j a p ϑ p a s ϑ s a q ϑ q + + ∂ 3 R ∂Θ 2 ∂ϑ j a p Θ p a s Θ s a q ϑ q + ∂ 3 R ∂Θ∂ϑ 2 j a p Θ p a s ϑ s a q ϑ q ] + . . . = 0 (4.66)
where R ( Θ j ; ϑ j ) = 0 since (Θ; ϑ) j belongs to the solution domain. Eq. 4.66 can be re-arranged in the following manner:

R (Θ; ϑ) = a [ ∂R ∂Θ j Θ 1 + ∂R ∂ϑ j ϑ 1 ] + a 2 { ∂R ∂Θ j Θ 2 + ∂R ∂ϑ j ϑ 2 + 1 2 [ ∂ 2 R ∂Θ 2 j Θ 2 1 + ∂ 2 R ∂ϑ 2 j ϑ 2 1 + ∂ 2 R ∂Θ∂ϑ j Θ 1 ϑ 1 ]} + a 3 [ ∂R ∂Θ j Θ 3 + ∂R ∂ϑ j ϑ 3 -F 3 (Θ 1 , Θ 2 ; ϑ 1 , ϑ 2 ) ] + . . . + +a n [ ∂R ∂Θ j Θ n + ∂R ∂ϑ j ϑ n -F n (Θ 1 , Θ 2 , . . . , Θ n-1 ; ϑ 1 , ϑ 2 , . . . , ϑ n-1 ) ] = 0 (4.67)
where, for a give order p of a, F p is a function of the unknowns Θ k and ϑ k up to the order p -1. Eq. 4.67 must be valid for a generic value of a, therefore, the following linear system of N equations in N + 1 unknowns:

∂R ∂Θ j Θ p + ∂R ∂ϑ j ϑ p = F p (Θ 1 , Θ 2 , . . . , Θ p-1 ; ϑ 1 , ϑ 2 , . . . , ϑ p-1 ) (4.68)
is obtained from each term a p . F p is known once the linear system derived form the orders lower than the p-th one have been solved. The problem in Eq. 4.68 is not well posed and a complementary condition for each expansion term p should be provided. This complementary condition is obtained from the definition of the path parameter a as a pseudo-arc-length by projecting the solution increment (Θ; ϑ) j+1 = ( Θ -Θ j ; ϑ -ϑ j ) over the tangent direction (Θ 1 ; ϑ 1 ):

a = ( Θ -Θ j ) T Θ 1 + ( ϑ -ϑ j ) ϑ 1 (4.69)
The following equations are, then, obtained by replacing Eqs. 4.64 within Eq. 4.69:

∥Θ 1 ∥ 2 + ϑ 2 1 = 1 Θ T p Θ 1 + ϑ p ϑ 1 = 0 p ≥ 2 (4.70)
being ∥□∥ the Euclidean norm. As a concluding remark, it should be noted that at each step j, the matrix ∂R ∂Θ j is the same for each expansion term p which allows to save a relevant amount of calculation time. Once the values Θ p and ϑ p have been all computed, the path solution at step j + 1 is obtained by Eq. 4.64. Finally, the maximum value of path parameter a is adopted with the same principle as in [START_REF] Cochelin | Méthode asymptotique numérique[END_REF]:

a = ( ϵ ∥Θ 1 ∥ ∥Θ p ∥ ) 1 p-1 , (4.71)
where ϵ is a precision parameter, which is determined by the user.

Chapter 5 Present Chapter corresponds to the published research paper (Huang et al., 2017) [Computer Methods in Applied Mechanics and Engineering, 318: 270-295, 2017.].

Fourier model for sandwich instability

A self-consistent notation is adopted.

Keywords: Sandwich plate, Fourier series, Instability, Wrinkling, Asymptotic Numerical Method.

Introduction

During recent decades, sandwich structures are widely used in aerospace, automotive and civil engineering fields. For the typical sandwich structures, the skins are usually made of high-strength materials such as lightweight metal alloys so that the skins can carry almost all the axial loads, whereas the core can be made of low-density materials such as carbon, balsa or plastic materials to transmit the transverse normal and shear loads. Thus very attractive properties in terms of light weight and high flexural stiffness can be obtained within this principle of construction. When designing this kind of structures, one of the major concerns belongs to the buckling phenomena (global and local instabilities). Hence, this Chapter aims to accurately and efficiently predict the critical load and describe the structural responses in these phenomena.

Instability phenomena in sandwich structures can be generally divided into three main classes: global buckling, antisymmetrical wrinkling and symmetrical wrinkling, see [START_REF] Allen | Analysis and Design of Structural Sandwich Panels[END_REF]. Over the last years, many works have been devoted to investigate and characterise the global and local instabilities of sandwich structures Léotoing et al. (2002a); [START_REF] Aiello | Local buckling loads of sandwich panels made with laminated faces[END_REF]; Ji andWaas (2007, 2010); Douville and Le Grognec (2013); Sad [START_REF] Saoud | A unified formulation for the biaxial local and global buckling analysis of sandwich panels[END_REF]; [START_REF] Noor | Computational Models for Sandwich Panels and Shells[END_REF]; Reddy (2004b); [START_REF] Dawe | Overall and local buckling of sandwich plates with laminated faceplates[END_REF]; [START_REF] Yuan | Overall and local buckling of sandwich plates with laminated faceplates, part II: Applications[END_REF]; Léotoing et al. (2002b). Generally, the skins are modelled by nonlinear beams or plates due to their thinness. For the core, there are many different ways to model it depending on the type of the core material and its deformability. The simplest model is the socalled one-parameter foundation that only considers the transversal core stiffness. Dawe (2001); Léotoing et al. (2002b), were systematically investigated and assessed and benchmark solutions were also presented for global and local buckling in sandwich structures. They showed that the quadratic approximation as in Léotoing et al. (2002a,b) for the transverse displacement in the core is the minimum requirement to well characterise the local buckling instability and thus leads to few variables. As an extension work of Léotoing et al. (2002a,b), Hu et al. (2009a) proposed a novel one-dimensional finite element considering the axial stiffness of the core, which yielded accurate results with a considerable reduction of computational cost for the global and local buckling of sandwich beams. In a following paper, [START_REF] Hu | Review and assessment of various theories for modelling sandwich composites[END_REF] also reviewed several theories for sandwich beams. As an extension work of Hu et al. (2009a), [START_REF] Yu | A novel two-dimensional finite element to study the instability phenomena of sandwich plates[END_REF] adopted the similar assumptions in Léotoing et al. (2002a,b) to develop a novel 2D finite element for instability analyses of sandwich plates, which has only nine variables and leads to fairly accurate and efficient solutions compared with the numerical results via the 3D finite element method (FEM). However, even for this 2D finite element, intensive computation efforts are still needed to accurately characterise the buckling responses of sandwich plates when the wrinkling wavelength is very small and the wavenumber is extremely large. Furthermore, the nonlinear calculation path is difficult to be controlled because too many solution paths exist around the useful one and the nonlinear system may be sensitive to the perturbation. Thus, it is necessary to develop highly efficient models to accurately and quickly describe instability phenomena in sandwich plates.

Recently, an effective double scale approach has been proposed based on the Fourier series with slowly varying coefficients (Damil andPotier-Ferry (2010, 2008)).

In this approach, exploiting the periodic nature of the instability pattern, the kinematic unknowns (addressed as "microscopic model") are expanded into the form of the Fourier series, which leads to a new problem with Fourier coefficients as new unknowns (addressed as "macroscopic model") that vary much more slowly than the original unknowns. This approach is similar to the Landau-Ginzburg technique but more robust than it, see [START_REF] Damil | Influence of local wrinkling on membrane behaviour: a new approach by the technique of slowly variable Fourier coefficients[END_REF]. So far, the Fourier-related approach has been successfully applied for instability phenomena in different kinds of structures such as the elastic beams on nonlinear Winkler foundation [START_REF] Damil | Influence of local wrinkling on membrane behaviour: a new approach by the technique of slowly variable Fourier coefficients[END_REF]; [START_REF] Mhada | About macroscopic models of instability pattern formation[END_REF]), sandwich beams [START_REF] Liu | A new Fourierrelated double scale analysis for instability phenomena in sandwich structures[END_REF]; [START_REF] Mhada | A 2D Fourier double scale analysis of global-local instability interaction in sandwich structures[END_REF]; [START_REF] Liu | The interactive bending wrinkling behaviour of inflated beams[END_REF]), thin metal sheets [START_REF] Damil | New nonlinear multi-scale models for wrinkled membranes[END_REF][START_REF] Damil | Membrane wrinkling revisited from a multiscale point of view[END_REF]; [START_REF] Huang | Macroscopic simulation of membrane wrinkling for various loading cases[END_REF]; [START_REF] Attipou | Thermal wrinkling of thin membranes using a Fourier-related double scale approach[END_REF]) and film-substrate systems (Huang et al. (2017b)). Numerical results verified that the approach has great advantages in reducing computational source and simplifying nonlinear calculation. As a reduced order model, the Fourier-based model becomes inaccurate near the boundary. It can be overcome by using the Arlequin method (Ben Dhia (1998); Ben [START_REF] Ben Dhia | The Arlequin method as a flexible engineering design tool[END_REF]), which couples the full model used in the region of boundary conditions and the macroscopic model in the bulk, see [START_REF] Hu | A bridging technique to analyze the influence of boundary conditions on instability patterns[END_REF]; [START_REF] Yu | Multi-scale techniques to analyze instabilities in sandwich structures[END_REF]. In this Chapter, based on the Fourier-related approach and the 2D microscopic finite element proposed by [START_REF] Yu | A novel two-dimensional finite element to study the instability phenomena of sandwich plates[END_REF], an effective 2D Fourier-based finite element is developed to study the symmetrical and antisymmetrical wrinkling phenomena in sandwich plates. The proposed 2D finite element model permits to capture the wrinkling evolution with much less degrees of freedom (DOFs) compared with the 2D microscopic model of [START_REF] Yu | A novel two-dimensional finite element to study the instability phenomena of sandwich plates[END_REF] and the 3D FEM solutions by commercial code ABAQUS. The obtained nonlinear system is solved by the Asymptotic Numerical Method (ANM), see [START_REF] Damil | A new method to compute perturbed bifurcations: Application to the buckling of imperfect elastic structures[END_REF] and [START_REF] Cochelin | Asymptotic-numerical methods and Padé approximants for non-linear elastic structures[END_REF]. The ANM is an efficient and robust path following technique in the pres-

The microscopical model

Kinematics

In this Chapter, a sandwich plate is considered as a three-layered panel with a soft core and two stiff skins, as shown in Fig. 5.1: x and y are the in-plane coordinates, z the transverse coordinate, L x the in-plane length and L y the in-plane width. The thickness of the top and bottom skin h s is assumed to be the same for simplicity.

The core and total thickness are h c and h t , respectively.

The same kinematics proposed by [START_REF] Yu | A novel two-dimensional finite element to study the instability phenomena of sandwich plates[END_REF] are used for the considered sandwich plate, where a high-order kinematic is used in the core and the classical plate theory (CPT) for the skins:

Top skin (5.2)

         U t (x, y, z) = u t (x, y) - ( z -hs+hc 2 ) ∂W t (x,y,z) ∂x , V t (x, y, z) = v t (x, y) - ( z -hs+hc 2 ) ∂W t (x,y,z) ∂y , hc 2 < z ⩽ ht 2 W t (x, y, z) = w t (x, y).
Bottom skin As the same assumption in Hu et al. (2009a); [START_REF] Liu | A new Fourierrelated double scale analysis for instability phenomena in sandwich structures[END_REF]; [START_REF] Yu | A novel two-dimensional finite element to study the instability phenomena of sandwich plates[END_REF]; Léotoing et al. (2002b), a linear shear stress in the thickness direction is considered in the core layer, which can well approximate the actual shear stress distribution in both antisymmetrical and symmetrical instability modes. Thus the total number of unknowns can firstly be reduced with the following relations:

         U b (x, y, z) = u b (x, y) - ( z + hs+hc 2 ) ∂W b (x,z) ∂x , V b (x, y, z) = v b (x, y) - ( z + hs+hc 2 ) ∂W b (x,y,z) ∂y , -ht 2 < z ⩽ -hc 2 W b (x, y, z) = w b (x, y), ( 5
   α * = -1 3 ψ ,x , β * = -1 3 ψ ,y ,
(5.4) where the ", x" stands for ∂ ∂x . As the skins and the core are assumed to be perfectly bounded to each other, the displacements at the interface between the skins (top or bottom) and the core should satisfy the continuity conditions

                           U t (x, y, hc 2 ) = U c (x, y, hc 2 ), U b (x, y, -hc 2 ) = U c (x, y, -hc 2 ), V t (x, y, hc 2 ) = V c (x, y, hc 2 ), V b (x, y, -hc 2 ) = V c (x, y, -hc 2 ), W t (x, y, hc 2 ) = W c (x, y, hc 2 ), W b (x, y, -hc 2 ) = W c (x, y, -hc 2 ).
(5.5) By substituting Eqs. (5.4) and (5.5) into Eqs. (5.1)-( 5.3), one gets the following relations:

                           w c = 1 2 (w t + w b ) -h 2 c 4 ψ, φ = 1 hc (w t -w b ), α = hs 2hc (w t ,x + w b ,x ) + 1 hc (u t -u b ) + h 2 c 12 ψ ,x , β = hs 2hc (w t ,y + w b ,y ) + 1 hc (v t -v b ) + h 2 c 12 ψ ,y , u c = hs 4 (w t ,x -w b ,x ) + 1 2 (u t + u b ) -h 2 c 4 u * , v c = hs 4 (w t ,y -w b ,y ) + 1 2 (v t + v b ) -h 2 c 4 v * .
(5.6)

Finally, the unknowns are reduced from seventeen to nine:

u t , u b , v t , v b , u * , v * , w t ,
w b , ψ.

Strain-displacement relation and constitutive laws

The constitutive and geometric equations are determined by the following hypotheses:

• the material behavior is isotropic in each layer and meets the Hooke's law;

• the geometrical non-linearity is considered in the skins only with moderate rotation and the strain-displacement relationship is described by the Von Karman strain;

• the core undergoes small displacements.

Thus, the strain-displacement relations for the skins and the core can be expressed as:

                                               ϵ s xx = U s ,x + 1 2 (W s ,x ) 2 , ϵ s yy = V s ,y + 1 2 (W s ,y ) 2 , 2ϵ s xy = U s ,y + V s ,x + W s ,x W s ,y , ϵ c xx = U c ,x , ϵ c yy = V c ,y , ϵ c zz = W c ,z , 2ϵ c xy = U c ,y + V c ,x , 2ϵ c yz = V c ,z + W c ,y , 2ϵ c xz = U c ,z + W c ,x ,
(5.7) in which the superscript s denotes either the top or bottom skin for brevity. For the skin layers, the stress-strain relationships are

       σ s xx σ s yy σ s xy        = E s 1 -ν 2 s     1 ν s 0 ν s 1 0 0 0 (1-νs) 2            ϵ s xx ϵ s yy 2ϵ s xy        , (5.8)
while the stress-strain relations for the core can be expressed as: (5.10) and

                     σ c xx σ c yy σ c zz σ c xy σ c yz σ c xz                      =            C 11 C 12 C 13 0 0 0 C 12 C 22 C 23 0 0 0 C 13 C 23 C 33 0 0 0 0 0 0 C 44 0 0 0 0 0 0 C 55 0 0 0 0 0 0 C 66                                 ϵ c xx ϵ c yy ϵ c zz 2ϵ c xy 2ϵ c yz 2ϵ c xz                      , (5.9) in which C 11 = C 22 = C 33 = λ c + 2G c , C 12 = C 13 = C 23 = λ c , C 44 = C 55 = C 66 = G c ,
λ c = E c ν c (1 + ν c )(1 -2ν c ) , G c = E c 2(1 + ν c ) .
(5.11)

In Eqs. (5.8-5.11), E and ν are respectively the Young's modulus and the Poisson's ratio, λ the Lame's first parameter and G the shear modulus. Eqs. (5.8) and (5.9) can be rewritten in the following compact form:

σ s = C s ϵ s , σ c = C c ϵ c .
(5.12)

Weak form of the microscopical model

The weak form of the governing equations is obtained by the Principal of Virtual Work (PVW):

δP int = δP ext (5.13)
where δ stands for a virtual variation, P int is the internal work and P ext the external work. The internal virtual work is expressed as: (5.14) where the superscript T represents transposition of a matrix. The volumes

δP int = δP c int +δP t int +δP b int = ∫ Vc δϵ T c σ c dV c + ∫ Vt δϵ T t σ t dV t + ∫ V b δϵ T b σ b dV b ,
V t = H t × Ω, V c = H c × Ω and V b = H b ×
Ω stand for top skin, core and bottom skin integration domains where

H t = { z : h c 2 ≤ z ≤ h t 2 } , H c = { z : - h c 2 ≤ z ≤ h c 2 } , H b = { z : - h t 2 ≤ z ≤ - h c 2 } and Ω = [0, L x ] × [ - L y 2 , L y 2 ]
. By taking into account of the strain-displacement and stress-strain relationships, one can easily get the unknown vector of the microscopic model: where the term □ f is related to the microscopic model. In the following, we will present the details to separately derive the virtual work of the core and the skins of the microscopic model.

⟨q f ⟩ = ⟨ u t u t ,x u t ,y v t v t ,x v t ,y u b u b ,x u b ,y v b v b ,x v b ,y u * u * ,x u * ,y v * v * ,x v

Virtual work of the core

By using the core strains in Eq. (5.7), the stresses in Eq. (5.9), the core kinematics in Eq. (5.2), the relations in Eqs. (5.4) and (5.6) and after integrating the first term of the right hand side of Eq. (5.14) through the thickness of the core, the internal virtual work of the core is expressed as: (5.16) in which {q c f } ∈ R 30×1 is the unknown vector of the core layer

δP c int = ∫ Ω ⟨ δq c f ⟩ [C c ] { q c f } dΩ,
⟨ q c f ⟩ = ⟨ u c ,x u c ,y v c ,x v c ,y α α ,x α ,y β β ,x β ,y u * u * ,x u * ,y v * v * ,x v * ,y α * α * ,x α * ,y β * β * ,x β * ,y w c ,x w c ,y φ φ ,x φ ,y ψ ψ ,x ψ ,y ⟩ ,
(5.17)

and [C c ] ∈ R 30×30 is the elasticity matrix of the core presented in Table 5.1. To relate the core layer's unknown {q c f } to the total unknown of the microscopic model {q f }, a transformation matrix [R c f ] ∈ R 30×36 is introduced, whose non-zero elements are listed in Table 5.2, then one gets

{q c f } = [R c f ]{q f }.
(5.18)

Thus the generalized elasticity matrix of the core

[D c f ] ∈ R 36×36 can be easily obtained by [D c f ] = [R c f ] T [C c ][R c f ].
Finally, the virtual work of the core layer Eq. (5.16) can be rewritten as

δP c int = ∫ Ω ⟨δq f ⟩ [ D c f ] {q f }dΩ.
(5.19)

Virtual work of the skins

By substituting Eq. (5.1) or (5.3) into Eq. (5.7), the strains of the skins (as the kinematic assumption for the top and bottom skins is similar, here we take the kinematic field of the top skin as example) finally read (5.20)

[R c f ] ∈ R 30×36 . (1,2)= 1 2 (1,8)= 1 2 (1,14)=-
             ϵ s xx = u s ,x -(z - h s + h c 2 )w s ,xx + 1 2 (w s ,x ) 2 , ϵ s yy = v s ,y -(z - h s + h c 2 )w s ,
Following the derivation in Section 5.2.3.1 and using Eqs. (5.8) and (5.20), the internal virtual work of the skins is: (5.21) where {γ s f } and {S s f } are, respectively, the generalized strain and stress vector of the skins, and {γ s f } is expressed by the unknown vector of the skins {q s f }: (5.22) in which

δP s int = ∫ Ω ⟨ δγ s f ⟩ { S s f } dΩ,
{ γ s f } = ( [H s f ] + 1 2 [A f (q s f )] ) { q s f } ,
[H s f ] =           
1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

           , [A f (q s f )] =           
0 0 0 0 w s ,x 0 0 0 0 0 0 0 0 0 w s ,y 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 w s ,y 0 w s ,x 0 0 0 0 0 0 0 0 0

           , ( 5
[R s f ] ∈ R 9×36 , one gets { q s f } = [R s f ] {q f } .
(5.25)

The details of the non-zero elements of [R The generalized stress vector {S s f } is expressed as: (5.26) in which the generalized elasticity matrix in the skin

{ S s f } = [D s f ]{γ s f },
[D s f ] reads [D s f ] = E s 1 -ν 2 s            h s ν s h s 0 0 0 0 ν s h s h s 0 0 0 0 0 0 h 3 s 12 ν s h 3 s 12 0 0 0 0 ν s h 3 s 12 h 3 s 12 0 0 0 0 0 0 (1-νs)hs 2 0 0 0 0 0 0 (1-νs)h 3 s 6            .
(5.27)

The internal virtual work of a skin layer finally reads

δP s int = ∫ Ω ⟨δq f ⟩ [R s f ] T ( [H s f ] T + [A f (q f )] T ) { S s f } dΩ.
(5.28)

Governing equations of microscopical model

The external virtual work is expressed as:

δP ext = λ f ∫ Ω ⟨δq f ⟩{F f }dΩ, (5.29)
where {F f } denotes the external load and λ f is the load proportional parameter.

Finally the governing equation of the microscopical model reads:

                         ∫ Ω ⟨ δq f ⟩ [ R t f ] T ( [ H t f ] T + [ A f ( q t f )] T ) { S t f } + ⟨ δq f ⟩ [ R b f ] T ( [ H b f ] T + [ A f ( q b f )] T ) { S b f } + ⟨ δq f ⟩ [ D c f ] {q f }dΩ = λ f ∫ Ω ⟨δq f ⟩ {F f }dΩ, { S t f } = [ D t f ] ( [ H t f ] + 1 2 [ A f ( q t f )] ) [ R t f ] {q f } , { S b f } = [ D b f ] ( [ H b f ] + 1 2 [ A f ( q b f )] ) [ R b f ] {q f } .
(

The macroscopical model

The general methodology

In this part, a macroscopic finite element model will be derived from the microscopical model by using the technique of the Fourier series with slowly varying Fourier coefficients. The main idea is to expand the unknowns of the microscopic model into the Fourier series, whose slowly varying coefficients are introduced as the new unknowns, i.e., the Fourier-based macroscopic model. In this study, we limit the expansion related to the rapid oscillations only in the x direction as in [START_REF] Damil | New nonlinear multi-scale models for wrinkled membranes[END_REF][START_REF] Damil | Membrane wrinkling revisited from a multiscale point of view[END_REF]; [START_REF] Huang | Macroscopic simulation of membrane wrinkling for various loading cases[END_REF]: (5.31) where U (x, y) = {u t u b v t v b u * v * w t w b ψ} represents the rapid oscillating unknowns of the microscopic model. The Fourier coefficient U j (x, y), the envelope for j th order, is the new macroscopical unknown. The U j (x, y) varies slowly over a period [x, x + 2π q ] of the oscillation. The wavenumber q is defined as (5.32) in which L x is the length of the sandwich plate and q 0 is the number of half waves: q 0 ∈ N and q 0 ≥ 2 that is a given number. By selecting the wavenumber q 0 , the instability mode could be controlled. In this Chapter, we consider a macroscopic model with three envelopes U 0 (x, y), U 1 (x, y) and U -1 (x, y), as pictured in Fig.

U (x, y) = +∞ ∑ j=-∞ U j (x, y)e ijqx ,
q = π L x q 0 ,
⟨q r ⟩ = ⟨ u t 0 u t 0,x u t 0,y v t 0 v t 0,x v t 0,y u b 0 u b 0,x u b 0,y v b 0 v b 0,x v b 0,y u * 0 u * 0,x u * 0,y v * 0 v * 0,x v * 0,y u * 1R u * 1R,x u * 1R,y u * 1I u * 1I,x u * 1I,y v * 1R v * 1R,x v * 1R,y v * 1I v * 1I,x v * 1I,y w t 1 w t 1,x w t 1,xx w t 1,y w t 1,yy w t 1,xy w b 1 w b 1,x w b 1,xx w b 1,y w b 1,yy w b 1,xy ψ 1 ψ 1,x ψ 1,xx ψ 1,y ψ 1,yy ψ 1,xy ⟩ . (5.39)
Then, the internal virtual work of the core in Eq. (5.35) can be rewritten as

δP c int = ∫ Ω ⟨δq f ⟩ 0 [ D c f ] {q f } 0 + ⟨δq f ⟩ 1 [ D c f ] {q f } -1 + ⟨δq f ⟩ -1 [ D c f ] {q f } 1 dΩ (5.40) = ∫ Ω ⟨δq f ⟩ 0 [ D c f ] {q f } 0 + 2⟨δq f ⟩ 1R [ D c f ] {q f } 1R + 2⟨δq f ⟩ 1I [ D c f ] {q f } 1I dΩ = ∫ Ω ⟨δq r ⟩ [D c r ] {q r }dΩ, in which the generalized macroscopic elastic matrix [D c r ] ∈ R 48×48 is [D c r ] = [T 0 ] T [ D c f ] [T 0 ] + 2[T R ] T [ D c f ] [T R ] + 2[T I ] T [ D c f ] [T I ] .
(5.41) in which, the matrices [D s r ], [H s r ], [A r (q s r )] can be easily assembled by the corresponding terms of Eqs. (5.23) and (5.27) in the microscopic model

[D s r ] =      D s f 0 0 0 1 2 D s f 0 0 0 1 2 D s f      [H s r ] =     H s f 0 0 0 2H s f 0 0 0 2H s f     [A r (q s r )] = 2      1 2 A f ( {q s f } 0 ) A f ({q s f } 1R ) A f ({q s f } 1I ) A f ({q s f } 1R ) A f ({q s f } 0 ) 0 A f ({q s f } 1I ) 0 A f ({q s f } 0 )      .
(5.47)

The internal virtual work of the skins in Eq. (5.42) then reads:

δP s int = ∫ Ω ⟨δγ s r ⟩{S s r }dΩ = ∫ Ω ⟨δq s r ⟩ ( [H s r ] T + [A r (q s r )] T ) {S s r }dΩ.
(5.48)

The skins macroscopic unknown vector {q s r } can be related to the total unknown vector of the macroscopic model {q r } by (5.49) where the transformation matrix

{q s r } = [R s r ]{q r }
[R s r ] ∈ R 27×48 is [R s r ] =     R s f T 0 R s f T R R s f T I     .
(5.50)

Finally, the macroscopic version of the virtual internal work of the skins can be

δP s int = ∫ Ω ⟨δq r ⟩[R s r ] T ( [H s r ] T + [A r (q r )] T ) {S s r }dΩ.
(5.51)

Governing equations of macroscopical model

Considering the load proportional to a parameter λ r and neglecting the body forces, the external virtual work is defined as

δP ext = λ r ∫ Ω ⟨δq r ⟩{F r }dΩ.
(5.52)

Similar to the parameters in the microscopical model, {F r } ∈ R 48×1 denotes the external load and λ r is the load proportional parameter. Finally the governing equation of the macroscopical model reads

                         ∫ Ω ⟨δq r ⟩ [ R t r ] T ( [ H t r ] T + [ A r ( q t r )] T ) { S t r } + ⟨δq r ⟩ [ R b r ] T ( [ H b r ] T + [ A r ( q b r )] T ) { S b r } dΩ + ∫ Ω ⟨δq r ⟩ [D c r ] {q r }dΩ = λ r ∫ Ω ⟨δq r ⟩ {F r }dΩ, { S t r } = [D s r ] ( [ H t r ] + 1 2 [ A r ( q t r )] ) [ R t r ] {q r } , { S b r } = [D s r ] ( [ H b r ] + 1 2 [ A r ( q b r )] ) [ R b r ] {q r } .
(5.53)

The unknowns in Eq. (5.53) denote to the slowly varying envelopes of the displacement field, which are completely different from those rapidly fluctuating displacement field in Eq. (5.30) or in [START_REF] Yu | A novel two-dimensional finite element to study the instability phenomena of sandwich plates[END_REF]. However, one can notice that the governing equations for the microscopic and macroscopic models are almost the same in the form, which brings a lot of conveniences to implement the new macroscopic model based on our previous work [START_REF] Yu | A novel two-dimensional finite element to study the instability phenomena of sandwich plates[END_REF]).

Discretization

In this part, the Finite Element Method (FEM) is used to solve Eq. (5.53). The in-plane displacements related to u, v, u * , v * are discretized by the 4-node Lagrange bi-linear shape functions while the transversal displacements w, ψ by the 4-node non-conforming rectangular functions, see Reddy (2004a); [START_REF] Zienkiewicz | The finite element method for analysis of elastic isotropic and orthotropic slabs[END_REF]; [START_REF] Melosh | Basis of derivation of matrices for the direct stiffness method[END_REF]. The non-conforming element can lead to accurate results with less degrees of freedom compared to the conforming element, and the membrane locking and shear locking effect is not obvious in this element, see [START_REF] Yu | A novel two-dimensional finite element to study the instability phenomena of sandwich plates[END_REF]; Reddy (2004a). The node numbering is shown in Fig. 5.3. The standard is small enough to ensure the accuracy of the results. Unless other statements, the following 3D FEM solution procedure is the same as mentioned above.

Figure 5.4: The configuration of the rectangular sandwich plate. Its four edges are simply supported and uni-axial loads are applied.

Table 5.5: Geometric and material parameters of the sandwich plate.

L x [m] L y [m] E s [Pa] ν c(s) E c /E s h c /h s L x /h t 0.5
0.25 6.9 × 10 10 0.3 10 -4 10 50

Validation and mesh convergence study

To validate the proposed 2D macroscopic model, numerical test with the same rectangular sandwich plate under uni-axial compression in [START_REF] Yu | A novel two-dimensional finite element to study the instability phenomena of sandwich plates[END_REF] is conducted. Meanwhile, the convergence with respect to the mesh is studied for the 2D macroscopic model. The obtained wrinkling pattern and critical load are mainly compared with the 3D FEM solution. As shown in Fig. 5.4, the rectangular sandwich plate is simply supported around four edges, which is realised by constraining the derivative of the envelope of the transversal displacement w s 1,x for the two edges of length L y and w s 1 for the other two edges in the 2D macroscopic model (similar way to apply boundary conditions can be found in the works [START_REF] Liu | A new Fourierrelated double scale analysis for instability phenomena in sandwich structures[END_REF]; [START_REF] Mhada | About macroscopic models of instability pattern formation[END_REF])), and the x-directional in-plane compressive loads are applied at the two ends of the skins that are perpendicular to the x axis. The geometric and material parameters presented in Table 5.5 are applied to all the numerical examples in this Chapter. Fig. 5.5 qualitatively presents the color map of the transversal displacement obtained by the 2D macroscopic model and the 3D FEM model, which is the lowest buckling mode for the considered sandwich configuration. We can see that the two models display the same instability pattern with five half waves in the x direction and one half wave in the y direction, which demonstrates that the 2D macroscopic model is capable to well characterise the wrinkling pattern. Table 5.6 presents 5.6) and are all in good agreement with that obtained by the 3D FEM model and the 2D microscopic model in [START_REF] Yu | A novel two-dimensional finite element to study the instability phenomena of sandwich plates[END_REF]. One notes that for the proposed 2D macroscopic model only 2 elements in the wrinkling direction can be sufficient to accurately predict the local wrinkling phenomenon. It leads to just 285 DOFs, which corresponds to only about 0.06% DOFs of the 3D FEM model and even saves about 82% DOFs compared with the 2D microscopic model of [START_REF] Yu | A novel two-dimensional finite element to study the instability phenomena of sandwich plates[END_REF]. This numerical test indicates the accuracy and efficiency of the 2D macroscopic sandwich model in describing the instability phenomenon of the sandwich plate. Besides, we would like to mention that the mesh size of the proposed model is independent on the wrinkling wavelength, which makes this model even more attractive in efficiency than the other two models when considering instability phenomena with very small wavelength or extremely large wavenumbers (see details in Section 5.5.3).

Antisymmetrical wrinkling

In this section, the antisymmetrical wrinkling with larger wavenumbers for the same sandwich plate (L x /L y = 2) in Fig. 5.4 will be investigated based on the proposed 2D macroscopic model. The material properties, geometrical data and boundary conditions remain the same as in Section 5.5.1. Fig. 5.7 shows the antisymmetrical wrinkling patterns with seven half waves in the x direction and two half waves in the y direction, which are obtained by the 2D macroscopic model and the 3D FEM model. The two wrinkling morphologies match fairly well. For the 3D FEM model, the mesh density is the same as in Section 5.5.1. Generally, at least 4 elements should be used in a period of the wrinkle to get reasonable results, and thus 100 elements in the wrinkling direction is accurate enough to characterise seven or more half waves in the high modes The proposed 2D macroscopic model only uses 8 elements (285 DOFs), which saves about 99.94% computational resource compared with the 3D FEM solution. Fig. 5.9 also presents the critical loads obtained by the two models for different high antisymmetrical modes, in which the half wave numbers are respectively 1 and 2 in the y direction and changing from 5 to 13 in the x direction. All the critical loads agree very well with each other, in which the maximum difference is only 1.9%. Here, we should emphasize that the total element number remains the same as 8 for all the considered buckling modes and only the wavenumber q 0 is changed in the 2D macroscopic model. This is a great advantage compared to the general microscopic models (such as the 3D FEM model or the 2D microscopic sandwich plate model of [START_REF] Yu | A novel two-dimensional finite element to study the instability phenomena of sandwich plates[END_REF]), especially when dealing with instability phenomena with very small wavelength. For microscopic models, the number of the elements is related to the wavenumber or the wavelength, and thus the increase in wavenumber or the decrease in wavelength will lead to a larger computational cost. Besides, special perturbation forces should be carefully introduced in the microscopic models in order to obtain a certain buckling mode, which causes difficulties in controlling the nonlinear calculation. The perturbation sensitivity of the buckling mode was observed in [START_REF] Liu | A new Fourierrelated double scale analysis for instability phenomena in sandwich structures[END_REF] and [START_REF] Léotoing | Using new closed-form solutions to set up design rules and numerical investigations for global and local buckling of sandwich beams[END_REF] where very tiny perturbations can shift global buckling into local wrinkling for a certain sandwich beam. Therefore, the proposed 2D macroscopic sandwich finite element model has high efficiency, robustness and accuracy to describe the antisymmetrical wrinkling in sandwich plates. Furthermore, these features provide a different way, from the numerical aspect, to determine the lowest buckling mode and the corresponding critical load of sandwich plates by calculating the critical load with different wavenumber q 0 and then choosing the lowest one.

Symmetrical wrinkling

The symmetrical wrinkling patterns are investigated for the same sandwich plate as in Fig. 5.4. Two different kinds of loading cases are considered in this part: the uniaxial compressive load (the same as Section 5.5.1) and the equi-biaxial compressive load. For the latter, the compressive loads are applied on both the bottom and top skins in the x and y directions, and the four edges of the sandwich plate are still simply supported. Due to the symmetry of this kind of wrinkling, the transversal displacements of the mid-plane of the core layer could be constrained to zero. This is realised by setting the enrichment function of the transversal displacement ψ equal to 0 in the 2D macroscopic model.

Symmetrical wrinkling for uni-axial compressive load case

The sandwich plate is subjected to the uni-axial compressive load in the x direction.

Here, the one-dimensional wrinkling pattern is investigated (see Fig. 5.10). For the 2D macroscopic model and the 3D FEM model, the mesh schemes keep the same as above (see Section 5.5.1). To trigger the post-buckling, the vertical perturbation forces F • 10 -5 are applied on the top skin along the middle line x = Lx 2 and the opposite ones with equal value are applied on the bottom skin. Fig. 5.10 presents the transversal displacement contour obtained by the two models, which has 25 half waves in the x direction. The wrinkling morphologies are quite similar. The transversal displacement near the bifurcation point for the middle line y = 0 of the sandwich plate is further depicted in Fig. 5.11, which is quite similar to the symmetrical wrinkling in sandwich beams analysed with efficient one-dimensional finite elements (see Hu et al. (2009a) and [START_REF] Liu | A new Fourierrelated double scale analysis for instability phenomena in sandwich structures[END_REF]). The symmetrical wrinkling in the sandwich plate with the thickness ratio h c /h s = 10 can be treated as a thin film resting on a thin substrate. According to [START_REF] Huang | Nonlinear analyses of wrinkles in a film bonded to a compliant substrate[END_REF], the analytical wavelength for the one-dimensional wrinkling is about 0.047 m in this kind of structure, which is close to the wavelength 0.04 m numerically predicted by the 2D macroscopic model. Fig. 5.12 presents the curve of the transversal 5.5 for geometric and material details) is under uni-axial load. 2%. Thus, if such wrinkling patterns are simulated using the microscopic model of [START_REF] Yu | A novel two-dimensional finite element to study the instability phenomena of sandwich plates[END_REF], the value and location of perturbation forces should be carefully applied to detect the useful buckling mode, which leads to significant difficulties in nonlinear calculation. Nevertheless, it can be quite easy for the 2D macroscopic model to distinguish and detect the useful buckling mode by just changing the wavenumber q 0 .

Symmetrical wrinkling for equi-biaxial compressive load case

In this part, the skins of the sandwich plate are subjected to the equi-biaxial compressive load in the x and y directions. The symmetrical wrinkling with checkerboard patterns occurs, which alleviates the in-plane compressive stresses in all directions and leads to lower energy than that of the one-dimensional wrinkling, see [START_REF] Chen | Herringbone buckling patterns of compressed thin films on compliant substrates[END_REF]. Fig. 5.14 presents the color map of the transversal displacement for the sandwich plate (one quarter of the plate near the upper left boundary is zoomed for more clarity), which has 17 half waves in the x direction and 9 half waves in the y direction. For the 2D macroscopic model, the mesh scheme N ex × N ey = 2 × 18 (1083 DOFs) is used. One notes that as the displacement field is expanded in the x direction, only 2 elements are used in this direction. The mesh scheme for the 3D FEM model remains unchanged (474939 DOFs). 5.5 under equi-biaxial compressive load with vast wavenumbers by the 2D macroscopic model.

Conclusion

In this Chapter, a Fourier-related double-scale analysis for the instability phenomena of sandwich plates was performed. Based on the 2D sandwich model proposed by [START_REF] Yu | A novel two-dimensional finite element to study the instability phenomena of sandwich plates[END_REF], in which the classical plate theory is used in the skins and high-order kinematics in the core, an effective 2D sandwich plate finite element model has been developed by using the technique of slowly variable Fourier coefficients. Compared to the former model, it reduces dramatically the number of degrees of freedom, thus saving significant computational cost, and makes nonlinear solution path controllable, especially when dealing with the wrinkling instability phenomena with vast wavenumbers. Antisymmetrical and symmetrical wrinkling in sandwich plates exposed to uni-axial and equi-biaxial compressive loads with four edges simply supported have been investigated. The results have been verified by the 3D FEM solution as well as the 2D microscopic results of [START_REF] Yu | A novel two-dimensional finite element to study the instability phenomena of sandwich plates[END_REF], which demonstrates that the proposed Fourier-related sandwich plate could accurately and fast predict the critical load and describe the post-bifurcation response for the considered instability phenomena.

The proposed model can be improved from the following aspects. One of the promising direction is to adopt the microscopic model on the boundaries and couple the proposed Fourier-based model and the microscopic model by the Arlequin method (Ben [START_REF] Dhia | Multiscale mechanical problems: the Arlequin method[END_REF][START_REF] Dhia | Global-local approaches: the Arlequin framework[END_REF]). This will improve the accuracy of the proposed model (especially for more complicated boundary conditions) without significantly enhancing the computation cost, see [START_REF] Hu | A bridging technique to analyze the influence of boundary conditions on instability patterns[END_REF]; [START_REF] Yu | Multi-scale techniques to analyze instabilities in sandwich structures[END_REF]; [START_REF] Hu | Multi-scale nonlinear modelling of sandwich structures using the Arlequin method[END_REF]. In addition, though we need to prescribe the wavenumber as an input, we could introduce several Fourier harmonics to develop a multi-harmonic model that is able to capture the most suitable buckling mode, which has been validated in our recent work (Huang et al. (2017b)). Finally, we propose to extend our model by taking account the mean field of the transversal displacement to model the coupled global and local buckling of sandwich plates, similar to the work of [START_REF] Liu | A new Fourierrelated double scale analysis for instability phenomena in sandwich structures[END_REF] for modeling the coupled instability in sandwich beams.

Chapter 6

Conclusion and future perspectives

In this thesis, the instability phenomena in three common engineering structures, i.e., membrane, film/substrate and sandwich structures are investigated, in which several kinds of accurate yet efficient multi-scale models are developed by using the Technique of Slowly Variable Fourier Coefficients (TSVFC). The established nonlinear systems are mainly solved by the Asymptotic Numerical Method (ANM) that is an efficient and robust path following technique in the presence of bifurcations.

So far, some specific results are achieved as follows:

• Based on the Von Karman thin plate equations, the two-dimensional doublescale macroscopic model for the membrane structure has been constructed in Chapter 2 by expanding the displacement of the Von Karman equations into Fourier series and taking the Fourier coefficients as the new unknowns of the macroscopic model. Due to the property of the slow variation of the new unknowns, the needed spatial meshes can be very coarse and are not related with the wrinkling length. The established model has also been implemented into the commercial package ABAQUS as user element, which believes to lay a foundation for the future study of the instability phenomena in complex membrane structures. Numerical examples have been performed to demonstrate that the macroscopic model is able to accurately and quickly describe the details of membrane wrinkling, even if the number of the wrinkling is small. By especially analyzing the classical problem of a rectangular membrane under uniaxial load, we found that 1) the membrane instability is very sensitive to boundary conditions, for instance, the difference in the critical values due to uniaxial tensile stress and uniaxial displacement stretch is as much as 6 times;

2) there exists a dimensionless parameter that is almost constant at the presence of the wrinkles independent of the boundary condition, the loading case and the geometry size and shape of the membrane, of which the establishment has important guiding significance for rational design of membrane structure, rapid prediction and effective prevention of membrane instability.

• The two-dimensional Fourier double-scale macroscopic model has been developed in Chapter 3, which has been proved having high computational effi-ciency and good numerical stability by studying the instability of a stiff film bounded to a soft substrate system that is under uniaxial compressive stress.

By virtue of the deformation property of the film/substrate (see [START_REF] Yang | A new family of finite elements for wrinkling analysis of thin films on compliant substrates[END_REF]), we further combine the Carrera's Unified Formulation (CUF) with the TSVFC to deduce a one-dimensional Fourier double-scale macroscopic model in Chapter 4. This numerical model demonstrates impressive computational efficiency under the condition of ensuring accuracy, which makes it an admirable tool for large-scale instability problems such as in [START_REF] Mei | Principles and applications of micro and nanoscale wrinkles[END_REF].

• In the framework of the TSVFC, different Fourier coefficients have been used

to deduce the Fourier model for film/substrate structures. It is pointed out that the in-plane displacement should be approximated by at least the zeroorder and first-order Fourier coefficients, otherwise the oscillation locking phenomenon will occur. In addition, the Fourier double-scale model, which simultaneously accounts for a number of higher order Fourier coefficients, quickly predicts the minimum critical load of the structure, and the result is independent of the set wavelength. This could be an alternative option to trace the lowest critical buckling load in the framework of the TSVFC.

• Based on the high-order laminate theory, the two-dimensional Fourier doublescale macroscopic model of sandwich plate has been constructed in Chapter 5. Antisymmetric and symmetric wrinkling in sandwich plates exposed to uni-axial and equi-biaxial compressive loads with four edges simply supported have been investigated. The results demonstrate the Fourier-related sandwich model could accurately and fast predict the critical load and describe the postbifurcation response for the instability phenomena in sandwich structures.

The present thesis leaves some topics that would be the possible research directions to address in the future:

• The Technique of Slowly Variable Fourier Coefficients is computationally efficient and stable. Subsequent work could be done to further study instabilities in the polar coordinate systems such as circular membranes, cylindrical shells, curved film/substrate structures.

• When the local boundary solution is concerned, the Fourier model, similar to all other reduced models, may lack enough accuracy near these areas due to the difficulty for accurately imposing the boundary conditions. To overcome

Figure 1

 1 Figure 1.1: Instability in thin metal sheet during cold rolling (Damil and Potier-Ferry (2010)) Figure 1.2: Instability in metal film under uni-axial tension

Figure 2 . 1 :

 21 Figure 2.1: At least two macroscopic fields are necessary to describe a nearly periodic response: the mean field and the amplitude of the fluctuation.

Fig. ( 2

 2 Fig. (2.2) presents a clamped rectangular membrane under biaxial tension-compression load. A uniform tensile load N y = 10 N/mm is firstly preloaded, then an increas-

Figure 2 . 2 :

 22 Figure 2.2: Rectangular membrane under biaxial load. Local wrinkles occurs in the membrane when the applied compressive stress increases. Fig. (2.3) presents the wrinkling patterns just after the bifurcation. It appears that the new reduced model can correctly predict the same shape and number of wrinkles as the full shell model does. A more quantitive analysis is shown in Fig. (2.4). One can obviously conclude that the new reduced model agrees very well with the full shell model near the bifurcation point where the average compressive load is N x ≈ -0.0900 N/mm.

  Fig. (2.5) shows that the new reduced model

Figure 2 . 3 :

 23 Figure 2.3: Wrinkling patterns near the bifurcation of a rectangular membrane under biaxial load for new reduced model (a), full shell model (b).

  Figure 2.4: Bifurcation path of a rectangular membrane under biaxial load with tensile load N y = 10 N/mm.

Figure 2 . 5 :

 25 Figure 2.5: Central cross section (Y=100 mm) of a rectangular membrane under biaxial load near bifurcation point, N y = 10 N/mm, N x = -0.0905 N/mm.

Figure 2 . 6 :

 26 Figure 2.6: Central cross section (Y=100 mm) of a rectangular membrane under biaxial load further after bifurcation point, N y = 10 N/mm, N x = -0.1024 N/mm.

Figure 2 . 7 :

 27 Figure 2.7: Details of the out-of-plane displacement near the boundary in Fig. 2.6.

Figure 2 . 8 :

 28 Figure 2.8: Rectangular membrane under uniaxial load.

Figure 2 . 9 :Figure 2

 292 Figure 2.9: Wrinkling patterns near the bifurcation of a rectangular membrane under uniaxial displacement load for full shell model (a) and new reduced model (b), N y = 39.2 N/mm.

  |σ x | is proportional to σ y , while after membrane buckling, |σ x | increases rapidly, so C raises rapidly too. Thus membrane instability could be observed through the change of the value of C.

Figure 2

 2 Figure 2.12: Compressive stress distribution near bifurcation with new reduced model for two load cases: (a) displacement load; (b) stress load.

Fig. ( 2

 2 Fig. (2.13) plots the out-of-plane displacement along the cross section in Y = 168 mm farther away from bifurcation point. As in previous numerical tests, they

Figure 2

 2 Figure 2.13: Central cross section (Y=168 mm) of rectangular membrane under uniaxial stress load (a), displacement load (b).

Figure 2

 2 Figure 2.14: Compressive stress distribution near bifurcation with full shell model for rectangular membrane under uniaxial displacement load: (a) 200 × 400; (b) 200 × 600.

Figure 2 .

 2 Figure 2.15: The evolution of K in a rectangular membrane under uniaxial displacement load with full shell model.

Fig. ( 2

 2 Fig. (2.16) proves that the new reduced model implemented into UEL gives the same bifurcation point as the full shell model. Besides, in the same computational condition, the CPU cost time is 422.9 seconds with implemented elements, much less than 40508 seconds needed by the full shell model. This means that the new macroscopic model yields very strong time reductions, especially when there are many wrinkles.

  model . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.2.1 Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.2.2 Geometric equations and constitutive law . . . . . . . . . . . 40 3.3 Macroscopic model . . . . . . . . . . . . . . . . . . . . . . . . 42 3.3.1 Internal virtual work of the film . . . . . . . . . . . . . . . . 44 3.3.2 Internal virtual work of the substrate . . . . . . . . . . . . . 46 3.3.3 External work and weak form of the governing equation . . . 48 3.3.4 Finite element discretisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.4.2 Convergence and computation cost . . . . . . . . . . . . . . . 52 3.4.3 Discussions on Fourier-based approach and oscillation locking 56 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

Figure 3

 3 Figure 3.2: At least two macroscopic envelopes are necessary to describe a nearly periodic response: the mean field and the amplitude of the fluctuation.

  By applying Eq. (3.15) and by virtue of the properties in Eqs. (3.22)-(3.25), the virtual internal work of the film becomes:

  al. (2012); Yu et al. (2015); Yang et al. (2015)) and thus are not presented in this work. The obtained macroscopic model is referred as "macroscopic model I".

Fig. 3 .Figure 3

 33 Fig. 3.4 presents the instability pattern in the neighbourhood of the bifurcation point as a colour map. The macroscopic model correctly predicts the shape and the number of wrinkles. For the ABAQUS model, a mesh N ex × N ez = 80 × 81 with 3.95 • 10 4 DOFs is used to model the film/substrate structure. N ex and N ez represent the element number along the x-and z-directions, respectively. The nonlinear post-buckling solution in the ABAQUS is performed by the Riks algorithm, where the eigenmode of interest scaled by a very small factor is introduced as

Figure 3

 3 Figure 3.5: Bifurcation curves for ABAQUS solution and two-dimensional macroscopic model I with different meshes, problem II.

Figure 3

 3 Figure 3.6: Axial displacement distribution along x at z = 0.95 mm after bifurcation point in the substrate. Problem II.

Figure 3 . 7 :

 37 Figure 3.7: Bifurcation curves, problem II, three different kinds of macroscopic model and mesh 1.

Figure 4 . 1 :

 41 Figure 4.1: Sketch of an elastic thin stiff film on a compliant substrate

  the core and bottom layers is ensured via the Constrained Variational Principle (CVP). By introducing the Lagrange multipliers µ = ⟨µ 1 , µ 2 ⟩ as fictitious gluing forces, it results in the constrain term L (µ, u):

  4.2.These two harmonics of the Fourier series have been used to describe instability phenomena in beams on Winkler's foundation[START_REF] Damil | Influence of local wrinkling on membrane behaviour: a new approach by the technique of slowly variable Fourier coefficients[END_REF]), sandwich beams[START_REF] Yu | Multi-scale techniques to analyze instabilities in sandwich structures[END_REF];[START_REF] Mhada | A 2D Fourier double scale analysis of global-local instability interaction in sandwich structures[END_REF]) and thin metal films[START_REF] Damil | New nonlinear multi-scale models for wrinkled membranes[END_REF][START_REF] Damil | Membrane wrinkling revisited from a multiscale point of view[END_REF];[START_REF] Huang | Macroscopic simulation of membrane wrinkling for various loading cases[END_REF]), and rather good estimations in wrinkling pattern and critical load have been observed.

Figure 4 . 2 :

 42 Figure 4.2: At least two macroscopic fields are necessary to describe a nearly periodic response: the mean field and the amplitude of the fluctuation.

  z) does not fluctuate along the x direction according to Section 3.1, i.e., U s (x, z) = U s 0 (x, z), the in-plane strain ϵ s xx and stress σ s xx show no oscillation: ϵ s xx = (ϵ s xx ) 0 , σ s xx = (σ s xx ) 0 . This would cause a spurious phenomenon called "oscillation locking" in the macroscopic film/substrate

  j ∈ [1, n j ]. Then substituting the Eqs. (4.20) and (4.21) into Eq. (4.49) gets:

  .3. The bottommost surface of the substrate (z = 0) is clamped. The kinematic boundary condition is obtained straightforwardly by imposing the nodal DOFs corresponding to u b 0 , (w b 0 ) j and (w b 0,x ) j equal to zero. The vertical displacements W is simply supported and equal to zero at x/L = 0 and 1. Two concentrated forces of equal modulus F and opposite verse are applied along the axial direction at the through-the-thickness mid point of the film.

Figure 4 . 3 :

 43 Figure 4.3: Sketch of the film/substrate system under compression forces

  the critical load parameter λ c and the corresponding half-wave number q 0 for problem I. The critical load and wavenumber decrease as the substrate becomes softer and thicker. For a given E s /E f , the critical load and wavenumber change quite little when h s /h f exceeds 10 2 where the substrate becomes quite thick. This shows that only limited part of the substrate near film affects on the instability phenomenon while the remaining part is unstrained and shows almost no influence on such phenomenon[START_REF] Huang | Nonlinear analyses of wrinkles in a film bonded to a compliant substrate[END_REF];[START_REF] Yang | A new family of finite elements for wrinkling analysis of thin films on compliant substrates[END_REF]). The 1D macroscopic model yields accurate results compared with the 2D full model. The results are obtained by a 1D macroscopic model with n c = 7, n b = 3 and n j = 1 (only zero and first order Fourier coefficients considered) and the mean field W 0 being zero.

  5, n b = 3 and n j = 1. In ABAQUS, convergence has been carefully examined, and 50 elements along both the axial and transversal directions are enough to simulate the considered instability phenomenon. The DOFs of this corresponding mesh is about 1.57 • 10 4 . The instability pattern for the problem II is presented in Fig. 4.4, where the 1D macroscopic model can correctly predict the wrinkling pattern compared with the 2D full model. Besides, an analysis of mesh convergence is carried out to ensure the accuracy of the proposed 1D macroscopic model. Fig. 4.5 presents the bifurcation curves for four mesh schemes (ne = 2, 4, 8, 16) of the 1D macroscopic model along the axial direction. The results converge to λ c = 0.0504 with only 2 macroscopic elements (99 DOFs), which is very close to the critical load λ c = 0.0495 obtained by the 2D full model (15705 DOFs). The 1D macroscopic model saves about 99% degrees of freedom and significantly decreases computation cost compared to the 2D full model.

  Fig. 4.6 presents the higher mode of problem II with 33 wavenumbers. The 2D full model needs 200 elements along axis direction (242805 DOFs) to get the convergent solution,

(

  2010);[START_REF] Liu | A new Fourierrelated double scale analysis for instability phenomena in sandwich structures[END_REF];[START_REF] Mhada | About macroscopic models of instability pattern formation[END_REF];[START_REF] Hu | A bridging technique to analyze the influence of boundary conditions on instability patterns[END_REF];[START_REF] Yu | Multi-scale techniques to analyze instabilities in sandwich structures[END_REF]. To determine the lowest critical load for a given structure, one should simulate several wrinkling patterns corresponding to different wavenumbers q 0 , and then choose the lowest one among them. Here, we introduce a new strategy to avoid the intensive simulation efforts by accounting several Fourier harmonics with different wavenumbers simultaneously in a single simulation.

Figure 4 . 4 :

 44 Figure 4.4: Transverse displacement 40 • w (10 -3 m) color map: (a) the 2D full model and (b) the 1D macroscopic model n c = 5, n b = 3 and n j = 1. Problem II.

Figure 4 . 5 :

 45 Figure 4.5: Bifurcation curves for four mesh schemes ("ne" is short for the number of element) of the 1D macroscopic model (n c = 5, n b = 3, n j = 1, h c = 2ℓ w ). Problem II.

Figure 4 . 6 :

 46 Figure 4.6: Transverse displacement 2•w (10 -3 m) color map with 33 wavenumbers: (a) 2D full model and (b) 1D macroscopic model n c = 5, n b = 3 and n j = 1. Problem II.

Figure 4 . 7 :

 47 Figure 4.7: Bifurcation curves by the 1D macroscopic model with 33 wavenumbers (n c = 5, n b = 3, h c = 2ℓ w ). Problem II.

Figure 4 . 8 :Figure 4

 484 Figure 4.8: Transverse displacement 40 • w (10 -3 m) color map: (a) 2D full model and (b) 1D macroscopic model n c = 5, n b = 3 and n j = 5. Problem III.

Figure 4

 4 Figure 4.10: Bifurcation curves by the 1D macroscopic model with q 0 = 4. Problem III.

Figure 4 .

 4 Figure 4.11: Transverse displacement of film by the 1D macroscopic model near bifurcation point: (a) q 0 = 2 and (b) q 0 = 4. Problem III.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.2.2 Strain-displacement relation and constitutive laws . . . . . . 102 5.2.3 Weak form of the microscopical model . . . . . . . . . . . . . 103 5.2.4 Governing equations of microscopical model . . . . . . . . . . 108 5.3 The macroscopical model . . . . . . . . . . . . . . . . . . . . . 108 5.3.1 The general methodology . . . . . . . . . . . . . . . . . . . . 108 5.3.2 Internal work of the core . . . . . . . . . . . . . . . . . . . . . 110 5.3.3 Internal work of the skins . . . . . . . . . . . . . . . . . . . . 113 5.3.4 Governing equations of macroscopical model . . . . . . . and mesh convergence study . . . . . . . . . . . . 118 5.5.2 Antisymmetrical wrinkling . . . . . . . . . . . . . . . . . . . 120 5.5.3 Symmetrical wrinkling . . . . . . . . . . . . . . . . . . . . . . 123 5.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 Abstract This Chapter presents a Fourier-related double scale analysis to study the instability phenomena of sandwich plates. By expanding the displacement field into Fourier series, the sandwich plate model proposed by Yu et al. (2015), using the classical plate theory in the skins and high-order kinematics in the core, is transformed into a new Fourier-based reduced two-dimensional sandwich plate model with the slowly varying Fourier coefficients as macroscopic unknowns. The resulting nonlinear equations are solved by the Asymptotic Numerical Method (ANM), which is very efficient and reliable to capture the bifurcation point and the post-buckling path in wrinkling analyses. Both antisymmetrical and symmetrical wrinkling for sandwich plates under uni-axial and equi-biaxial compressive loads are studied and the numerical results demonstrate that the Fourier-based finite element model can accurately yet efficiently predict wrinkling patterns and critical loads, especially when dealing with wrinkling phenomena with extremely large wavenumbers.

  ence of bifurcation compared to the classical nonlinear solution methods such as the Newton-Raphson or the arc-length (Riks) method. Numerical examples are performed to assess the 2D Fourier-based finite element model and very accurate results are obtained with significant reduction in computational cost.The Chapter is structured as follows. In Section 5.2, the same kinematic assumptions as[START_REF] Yu | A novel two-dimensional finite element to study the instability phenomena of sandwich plates[END_REF] are used to derive a 2D microscopic model for sandwich plates. Then, on the basis of the 2D microscopic model and using the Fourierrelated approach, the weak form of governing equations of the Fourier-based 2D macroscopic model is derived in Section 5.3. In Section 5.4, the finite element method is used to discretize the proposed 2D Fourier-based model. In Section 5.5, numerical examples relating to the symmetrical and antisymmetrical wrinkling are performed to examine and assess the validation and efficiency of the 2D Fourierbased finite element model.

Figure 5 . 1 :

 51 Figure 5.1: Sketch of the sandwich plate.

  x, y, z) = u c (x, y) + zα(x, y) + z 2 u * (x, y) + z 3 α * (x, y),V c (x, y, z) = v c (x, y) + zβ(x, y) + z 2 v * (x, y) + z 3 β * (x, y), -hc 2 < z ⩽ hc 2W c (x, y, z) = w c (x, y) + zφ(x, y) + z 2 ψ(x, y).

  .3) in which, the superscripts t, b and c represent the top, bottom and core layers, respectively. The symbols U and V denote the in-plane displacements in the x and y directions and W denotes the transverse displacement in the z direction, and u, v and w, respectively, refer to the in-plane and transversal displacements of the neutral plane in each layer. The unknown functions α, u * , α * , β, v * , β * , φ and ψ are enrichment functions with respect to the classical zig-zag theory. One notes that both symmetrical and antisymmetrical instability modes can be described by this kinematic assumption, because the displacement field is separately assumed for each layer and a quadratic through-the-thickness function is used for the transverse displacement of the core layer. For the present, there are seventeen unknowns in this kinematic assumption.

  .30) This model permits to accurately and efficiently describe the global and local wrinkling of sandwich plates Yu et al. (2015). It will be referred to as 'microscopical model' in this Chapter, because the displacement unknowns in this model vary rapidly compared with the unknowns (the envelopes of the displacement) of the forthcoming Fourier-based model (macroscopic model), and thus the mesh sizes of the microscopic model should be much smaller than the wavelength to accurately describe local oscillations. Based on this microscopical model, the technique of slowly varying Fourier coefficients will be used to derive the macroscopical model in the next part.

5. 2 .Figure 5 . 2 :

 252 Figure 5.2: At least two macroscopical fields are necessary to describe a nearly periodic response: the mean field and the amplitude of fluctuation.

Figure 5

 5 Figure 5.5: The wrinkling patterns obtained by (a) the 2D macroscopic sandwich plate model and (b) the 3D FEM model.

Figure 5 . 6 :

 56 Figure 5.6: Bifurcation analysis for validation and mesh convergence of the 2D macroscopic sandwich model.

Figure 5

 5 Figure 5.7: The antisymmetrical wrinkling patterns obtained by (a) the 2D macroscopic sandwich plate model and (b) the 3D FEM model.

Figure 5

 5 Figure 5.8: Bifurcation path for the antisymmetric wrinkling with 7 half waves in the x direction and 2 half waves in the y direction.

Figure 5 . 9 :

 59 Figure 5.9: Critical loads for different antisymmetric modes with the 3D FEM model and the 2D macroscopic sandwich model. "y-wave" is short for half-wave numbers in the y direction.

Figure 5

 5 Figure 5.10: The symmetrical wrinkling patterns obtained by (a) the 2D macroscopic sandwich plate model and (b) the 3D FEM model. The sandwich plate (see Table5.5 for geometric and material details) is under uni-axial load.

Figure 5 .

 5 Figure 5.11: The transversal displacement along the line y = 0 for the symmetrical wrinkling by the 2D macroscopic model. The sandwich plate is under uni-axial load.

Fig. 5 .Figure 5

 55 Fig.5.13 displays the bifurcation curves for the symmetrical wrinkling with different half waves in the x direction. One notices that the difference of the critical loads for the five cases is quite small, where the maximum discrepancy is less than

Figure 5

 5 scheme for the 3D FEM model remains unchanged (474939 DOFs). Fig.5.15 shows the curve of the transversal displacement at the middle point x = L x /2, y = 0

Figure 5 .

 5 Figure 5.14: The symmetrical checkerboard wrinkling pattern for the sandwich plate of Table 5.5 under equi-biaxial compressive load by the 2D macroscopic model. One quarter of the plate near the upper left boundary is zoomed.

Figure 5 .

 5 Figure 5.15: The bifurcation path for the symmetrical checkerboard wrinkling pattern of the sandwich plate under equi-biaxial compressive load.

Figure 5 .

 5 Figure 5.16: The symmetrical checkerboard wrinkling pattern for the sandwich plate of Table 5.5 under equi-biaxial compressive load with vast wavenumbers by the 2D macroscopic model.

Table 2

 2

	.1: Mesh details for the rectangular membrane under biaxial load.
		X	Y Total element DOF
	Full shell model (ABAQUS) 200 20	4000	98646
	New reduced model	10	16	160	1599

Table 2

 2

	.2: Mesh details for rectangular membrane under uniaxial load.
		X	Y Total element DOF
	Full shell model(ABAQUS) 20 50	1000	24846
	New reduced model	8	30	240	2391

Table 2

 2

	.3. In these numerical

Table 3 .

 3 Table3.2 presents the results of the two parameters obtained by the macroscopic model I and 2D full model. Good results agreement is observed for all examples. This confirms the validation of the proposed macroscopic FE model in predicting the critical load and wavenumber in the film/substrate systems. From Table3.2, one notices that the critical load and the wavenumber decrease as the substrate becomes thicker and softer. For a fixed E 2: Half-wave number q 0 and critical load parameter λ c . "Present" stands for Fourier-series macroscopic model I and "ABAQUS" is 2D full model by ABAQUS.

	Problem I.
	h s /h f	10			10 2			10 3
		q 0	10 -4 λ c (MPa)	q 0	10 -4 λ c (MPa)	q 0	10 -4 λ c (MPa)
	E s /E f		Present ABAQUS		Present ABAQUS		Present ABAQUS
	10 -3	12	143	143	9	108	109	9	111	109
	10 -4	8	39.5	38.5	5	23.6	23.4	5	23.7	23.4
	10 -5	4	12.8	12.8	2	5.91	5.85	2	5.82	5.76

s /E f ratio, the results change very little when h s /h f varies from 10 2 to 10 3 , because only a limited part of the substrate near the interface contributes to the instability phenomenon while the remaining part far from the interface is unstrained. This is also evident from the wrinkling color map of a thin film resting on a very thick substrate in Fig.

3

.4. These prominent features in the wrinkling phenomenon of a stiff film on an elastic compliant substrate are similar to that revealed by

[START_REF] Huang | Nonlinear analyses of wrinkles in a film bonded to a compliant substrate[END_REF]

through analytical solutions and

[START_REF] Yang | A new family of finite elements for wrinkling analysis of thin films on compliant substrates[END_REF]

through numerical results with finite element method.

Table 3

 3

	.4. The DOFs

Table 3

 3

.4: Mesh details for problem II.

Table 4 .

 4 1: Material and geometric parameters of the film/substrate. Problem I.

Table 5 .

 5 2: Non-zero elements of the transformation matrix

Table 5 .

 5 4: Relevant transformation matrices for the core layer, in which the non-zero elements are listed.

	Matrix Size
	[T R ]	36 × 48 (13,19)=1, (14,20)=1, (14,22)=-q, (15,21)=1, (16,25)=1, (17,26)=1,
		(17,28)=-q, (18,27)=1, (19,31)=1, (20,32)=1, (21,33)=1, (21,31)=-q 2 ,
		(22,34)=1, (23,35)=1, (24,36)=1, (25,37)=1, (26,38)=1, (27,39)=1,
		(27,37)=-q 2 , (28,40)=1, (29,41)=1, (30,42)=1, (31,43)=1, (32,44)=1,
		(33,45)=1, (33,43)=-q 2 , (34,46)=1, (35,47)=1, (36,48)=1.
	[T I ]	36 × 48 (13,22)=1, (14,19)=q, (14,23)=1, (15,24)=1, (16,28)=1, (17,25)=q,
		(17,29)=1, (18,30)=1, (20,31)=q, (21,32)=2q, (24,34)=q, (26,37)=q,
		(27,38)=2q, (30,40)=q, (32,43)=q, (33,44)=2q, (36,46)=q.
	[T0]	36 × 48 (1,1)=1, (2,2)=1, (3,3)=1, (4,4)=1, (5,5)=1, (6,6)=1, (7,7)=1, (8,8)=1,
		(9,9)=1, (10,10)=1, (11,11)=1, (12,12)=1, (13,13)=1, (14,14)=1, (15,15)=1,
		(16,16)=1, (17,17)=1, (18,18)=1.

  N ex N ey N ez Total elements DOFs Critical load (10 4 N/m)

	Shell/Volume/Shell	100 50	6	30000	474939	2.02
	2D Microscopic Model	20	4	-	80	1575	2.00
	2D Macroscopic Model	2	4	-	8	285	2.02
	2D Macroscopic Model	4	4	-	16	475	2.02
	2D Macroscopic Model	8	4	-	32	855	2.02
	Table 5.6: The mesh scheme for the 3D FEM model, the 2D microscopic model and
	the 2D macroscopic model under the uni-axial compressive load.
		2.5	-4 x 10 2D macroscopic sandwich plate (8 elements, 285 DOFs)
			2D macroscopic sandwich plate (16 elements, 475 DOFs)
			2D macroscopic sandwich plate (32 elements, 855 DOFs)
	Maximum transversal displacement (m)	0.5 1 1.5 2	Shell/Volume/Shell (30000 elements, 474939 DOFs) 2D microscopic sandwich plate (80 elements, 1575 DOFs)
		0	0	0.5		1	1.5	2
						Load F (N/m)		x 10 4

Acknowledgment

This thesis has been written during my position as a joint-

Table 5.1: Non-zero elements of the elasticity matrix C c ∈ R 30×30 . (1,1)=h1E ′ (1,4)=h1λ (1,12)=h3E ′ (1,16)=h3λ (1,25)=h1λ (2,2)=h1Gc (2,3)=h1Gc (2,13)=h3Gc (2,15)=h3Gc (3,2)=h1Gc (3,3)=h1Gc (3,13)=h3Gc (3,15)=h3Gc (4,1)=h1λ (4,4)=h1E ′ (4,12)=h3λ (4,16)=h3λ (4,25)=h1λ (5,5)=h1Gc (5,17)=3h3Gc (5,23)=h1Gc (5,29)=h3Gc (6,6)=h3E ′ (6,10)=h3λ (6,18)=h5E ′ (6,22)=h5λ (6,28)=2h3λ (7,7)=h3Gc (7,9)=h3Gc (7,19)=h5Gc (7,21)=h5Gc (8,8)=h1Gc (8,20)=3h3Gc (8,24)=h1Gc (8,30)=h3Gc (9,7)=h3Gc (9,9)=h3Gc (9,19)=h5Gc (9,21)=h5Gc (10,6)=h3λ (10,10)=h3E ′ (10,18)=h5λ (10,22)=h5E ′ (10,28)=2h3λ (11,11)=4h3Gc (11,26)=2h3Gc (12,1)=h3E ′ (12,4)=h3λ (12,12)=h5E ′ (12,16)=h5λ (12,25)=h3λ (13,2)=h3Gc (13,3)=h3Gc (13, 13) = h5Gc (13, 15) = h5Gc (14,14)=4h3Gc (14,27)=2h3Gc (15,2)=h3Gc (15,3)=h3Gc (15, 13) = h5Gc (15, 15) = h5Gc (16,1)=h3λ (16,4)=h3E ′ (16,12)=h5λ (16,16)=h5E ′ (16,25)=h3λ (17,5)=3h3Gc (17,17)=9h5Gc (17,23)=3h3Gc (17,29)=3h5Gc (18,6)=h5E ′ (18,10)=h5λ (18,18)=h7E ′ (18,22)=h7λ (18,28)=2h5λ (19,7)=h5Gc (19,9)=h5Gc (19,19)=h7Gc (19,21)=h7Gc (20,8)=3h3Gc (20,20)=9h5Gc (20,24)=3h3Gc (20,30)=3h5Gc (21,7)=h5Gc (21,9)=h5Gc (21,19)=h7Gc (21,21)=h7Gc (22,6)=h5λ (22,10)=h5E ′ (22,18)=h7λ (22,22)=h7E ′ (22,28)=2h5λ (23,5)=h1Gc (23,17)=3h3Gc (23,23)=h1Gc (23,29)=h3Gc (24,8)=h1Gc (24,20)=3h3Gc (24,24)=h1Gc (24,30)=h3Gc (25,1)=h1λ (25,4)=h1λ (25,12)=h3λ (25,16)=h3λ (25,25)=h1E ′ (26,11)=2h3Gc (26,26)=h3Gc (27,14)=2h3Gc (27,27)=h3Gc (28,6)=2h3λ (28,10)=2h3λ (28,18)=2h5λ (28,22)=2h5λ (28,28)=4h3E ′ (29,5)=h3Gc (29,17)=3h5Gc (29,23)=h3Gc (29,29)=h5Gc (30,8)=h3Gc (30,20)=3h5Gc (30,24)=h3Gc (30,30)=h5Gc

Abstract

In this Chapter, a new Fourier-related double scale approach is presented to study the wrinkling of thin films on compliant substrates. By using the method of Fourier series with slowly variable coefficients, the 1D microscopic model proposed by [START_REF] Yang | A new family of finite elements for wrinkling analysis of thin films on compliant substrates[END_REF] is transformed into a 1D macroscopic film/substrate model whose mesh size is independent on the wrinkling wavelength. Numerical tests prove that the new model improves computational efficiency significantly with accurate results, especially when dealing with wrinkling phenomena with vast wavenumbers. Besides, we propose a strategy to efficiently trace the wrinkling pattern corresponding to the lowest critical load by accounting for several harmonics of Fourier series in this new model. The established nonlinear system is solved by the Asymptotic Numerical Method (ANM), which has advantages of efficiency and reliability for stability Then improvement is made in the two-parameter model by accounting for both the compressibility of the core and the shear effect between the skin and the core [START_REF] Aiello | Local buckling loads of sandwich panels made with laminated faces[END_REF]. These models can be only practical on the condition that the core thickness is relatively thick compared to the skins thickness, i.e., the interaction between the two skins is negligible. To adequately describe the kinematic behavior in the core, the core can be modeled as a three-dimensional (3D) continuous solid.

On the basis of the elastic bifurcation theory in the 3D framework, the closed-form analytical expressions of the critical loads and the related buckling modes (global and local buckling) are recently presented for sandwich plates and beams, see Douville and Le Grognec (2013); Sad [START_REF] Saoud | A unified formulation for the biaxial local and global buckling analysis of sandwich panels[END_REF]. However, from the numerical aspects, when the 3D solid is decretized by hexahedral elements [START_REF] Madenci | The finite element method and applications in engineering using ANSYS ®[END_REF], very fine meshes are needed to capture the response due to the small-scale wrinkles and the relatively large thickness, which leads to large computational cost. Thus, to accurately characterise sandwich structures' mechanical behavior, including buckling phenomena, in a computation-cost-effective way, many works are resorted to the two-dimensional (2D) modelling framework, where the 3D problem reduces to the 2D problem by making suitable assumptions with respect to the displacement or stress field in the thickness direction. Generally these works can be classified into two main approaches: the equivalent single layer (ESL) approach and the layer-wise (LW) one, see Reddy (2004b). For both approaches, various additional functions are usually introduced (such as linear function [START_REF] Kreja | Finite elements based on a first-order shear deformation moderate rotation shell theory with applications to the analysis of composite structures[END_REF]; Zhang et al. (2016a,b)), cubic function [START_REF] Reddy | A simple higher-order theory of laminated composite plate[END_REF]; [START_REF] Ferreira | Analysis of composite plates using a layerwise theory and multiquadrics discretization[END_REF]), sinusoidal function [START_REF] Touratier | An efficient standard plate theory[END_REF]; [START_REF] Vidal | Assessment of the refined sinus model for the non-linear analysis of composite beams[END_REF]), exponential function [START_REF] Karama | Mechanical behaviour of laminated composite beam by the new multi-layered laminated composite structures model with transverse shear stress continuity[END_REF]) and Legendre's polynomials (D' Ottavio and Polit (2015))

etc.) to cooperate for multilayer structures with different geometrical, material properties under various loading cases over the years. Comprehensive assessments and comparisons on various approaches and models used for modelling multilayered plates as well as sandwich plates can be found in Reddy (2004b) and[START_REF] Carrera | Historical review of Zig-Zag theories for multilayered plates and shells[END_REF]. From the aspect of numerical simulation, it is necessary to balance accuracy and computation cost for a certain model when handling buckling problems in multilayered plates or sandwich plates.

Based on the Carrera Unified Formulation (CUF), [START_REF] Carrera | Selection of appropriate multilayered plate theories by using a genetic like algorithm[END_REF] proposed the so-called Best Theory Diagram to select the appropriate unknown variables upon application of genetic algorithms. Recently, D 'Ottavio et al. (2016) proposed an equally clever high order LW model named sublaminate generalized unified formulation (S-GUF), in framework of which several kinematic models, e.g. high order models in Léotoing et al. (2002a); [START_REF] Dawe | Overall and local buckling of sandwich plates with laminated faceplates[END_REF]; Yuan and which means that their expansion reduces into a single term: y), and F (x, y) = F 0 (x, y). As the in-plane displacements of the core in fact show fluctuation in the wrinkling direction when the sandwich plates buckle, three Fourier terms should be retained for the enrichment functions of the core's in-plane displacement kinematics:

(5.33)

Secondly, as we are interested in the local instability of the sandwich plates, the mean field of the transverse displacement W 0 (x, y, z) for both the skins and the core is eliminated, i.e., W 0 (x, y, z) = 0. Thus, the transverse displacement W (x, y, z) can be:

where the envelope W 1 is real and equals to W -1 . The phase ϕ could be variant so as to describe the phase modulation in the evolution of instability, see [START_REF] Mhada | About macroscopic models of instability pattern formation[END_REF]. To simplify our work, ϕ is assumed to be constant and equals to ± π 2 , which is simple yet accurate for the general cases (see [START_REF] Liu | A new Fourierrelated double scale analysis for instability phenomena in sandwich structures[END_REF]; [START_REF] Yu | Multi-scale techniques to analyze instabilities in sandwich structures[END_REF]; [START_REF] Huang | Macroscopic simulation of membrane wrinkling for various loading cases[END_REF]) but may lead to a loss of accuracy near the boundary (see [START_REF] Mhada | About macroscopic models of instability pattern formation[END_REF]; [START_REF] Huang | Macroscopic simulation of membrane wrinkling for various loading cases[END_REF]). For the model exactly treating the boundary condition, one can resort to a multi-scale model similar to [START_REF] Hu | A bridging technique to analyze the influence of boundary conditions on instability patterns[END_REF].

Finally, one can get the envelopes of the displacement unknowns of the micro-

One notes that if the zero harmonic of the transverse displacement is taken into account, the coupling of global and local buckling of sandwich plates (studied by a similar high order Fourier model in [START_REF] Liu | A new Fourierrelated double scale analysis for instability phenomena in sandwich structures[END_REF] for sandwich beams) can also be characterised.

Internal work of the core

By using the derivation rules for Fourier series with slowly varying coefficients in [START_REF] Damil | Influence of local wrinkling on membrane behaviour: a new approach by the technique of slowly variable Fourier coefficients[END_REF], the internal virtual work of the core in Eq. (5.19) is expanded into the Fourier form:

(5.35)

The unknown vector {q f } j is the envelope corresponding to the j th harmonic of the microscopic unknowns {q f } and their relation reads

(5.36)

As an example, we write

) 0 and

) 0 in the following form:

(5.37)

According to Eq. (5.37), the envelopes {q f } ±1 can be written as: (5.38) in which {q f } ±1 are divided into two parts: the real one {q f } 1R and the complex one i {q f } 1I . Then one can get the following relation: and their non-zero elements are listed in Table 5.4. According to the assumptions in Section 5.3.1, the microscopic unknowns {q f } in Eq. (5.15) is transformed into

Internal work of the skins

Similar to Section 5.3.2, the internal virtual work of the skins in Eq. ( 5.21) is also expressed in the following Fourier form: (5.42) in which the envelopes of the generalized stress vector {S s f } in Eq. ( 5.26) are (5.43) and the envelopes of the generalized strain vector {γ s f } in Eq. ( 5.22) are written as

(5.44)

By introducing the following generalized macroscopic stress {S s r }, strain {γ s r } ∈ R 18×1 and displacement gradient {q s r } ∈ R 27×1 in skins: , (5.45) the generalized macroscopic stress-strain and strain-displacement relations of the skins can be consequently rewritten into the similar form as that in the microscopic model (see Eqs. (5.22) and (5.26)):

) {q s r },

(5.46) bi-linear isotropic shape functions read as

where the parameters ξ and η are the elementary local coordinates.

Figure 5.3: The element natural reference system and nodes coordinates.

The non-conforming element consists of three degrees of freedom (w, ∂w ∂x , ∂w ∂y ) at each vertex node. The interpolation functions are

where (5.56) in which i = 1, 2, 3, 4, ξ 0 = ξ i ξ, η 0 = η i η, ξ i and η i are the corresponding node coordinates in the natural reference system. Then, the macroscopical discretized unknowns could be written as (5.57) where [G r ] is the interpolation matrix that is assembled by N L , N H and their derivation in terms of the local coordinates. The vector

is the nodal unknowns of an element with 19 degrees of freedom per node:

(5.58)

Finally, by substituting Eq. ( 5.57) into Eq. (5.53) the discretized form of the governing equation for the macroscopical model is expressed as:

(5.59)

The resulting non-linear problem Eq. (5.59) is solved by the Asymptotic Numerical Method (ANM) [START_REF] Cochelin | Asymptotic-numerical methods and Padé approximants for non-linear elastic structures[END_REF]. For each element, the macroscopic unknowns are determined at the nodes and the stress at the Gauss points. A short brief introduction on the use of the ANM for the resolution of non-linear problems is given in Appendix A.

Numerical results

In this part, we firstly investigate the accuracy and convergence of the proposed macroscopic or Fourier-based model for the instability phenomena in sandwich plates by comparing its results with those of other models, and then extend the model to simulate the antisymmetrical and symmetrical local wrinkling with larger wavenumbers. To verify the macroscopic model, reference results are obtained by constructing a 3D finite element model in ABAQUS, which are referred to as "Shell/Volume/Shell" or "3D FEM model". In this model, the 8-node shell element S8R is used for the skins and 20-node brick element C3D20R for the core layer, see [START_REF] Yu | A novel two-dimensional finite element to study the instability phenomena of sandwich plates[END_REF]. The skins are bonded to the core layer using "tie" constraints. For the 3D FEM solution, the nonlinear post-buckling analysis is performed by a loaddisplacement algorithm (RIKS) in ABAQUS. To trigger the interested buckling instability during simulation, the eigenmodes obtained from the linear perturbation analysis (BUCKLE) is introduced as the initial geometric imperfection, which this shortcoming, the multi-scale model will be a remedy that could deal well with balancing the local accuracy and the global computational efficiency, in which the enriched model such as the shell element is adopted near the boundary conditions to capture the boundary effects and the macroscopic model such as the Fourier model is used in the rest to reduce computational cost, the two models are bridged by the bridging domain method (Arlequin method).

• In this thesis, the dimensionless parameter is observed from the numerical aspects. Follow-up work could be further carried out by conducting experiments of thin film to determine the existence of this dimensionless parameter.

• The occurrence of wrinkles in the film/substrate is due to the compressive stress. With the increasing of the stress, the instability morphology may evolve from sinusoidal wrinkles into period-doubling [START_REF] Brau | Multiple-length-scale elastic instability mimics parametric resonance of nonlinear oscillators[END_REF]; Hutchinson ( 2013)), which would be of interest using the TSVFC to study a film/substrate system with a non-linear substrate that is responsible for the occurrence of the period-doubling phenomenon.