CELLULAR INSTABILITIES ANALYZED BY MULTI-SCALE FOURIER SERIES: A REVIEW

. The paper is concerned by multi-scale methods to describe instability pattern formation, especially the method of Fourier series with variable coeﬃcients. In this respect, various numerical tools are available. For instance in the case of membrane models, shell ﬁnite element codes can predict the de- tails of the wrinkles, but with diﬃculties due to the large number of unknowns and the existence of many solutions. Macroscopic models are also available, but they account only for the eﬀect of wrinkling on membrane behavior. A Fourier-related method has been introduced in order to modelize the main features of the wrinkles, but by using partial diﬀerential equations only at a macroscopic level. Within this method, the solution is sought in the form of few terms of Fourier series whose coeﬃcients vary more slowly than the oscillations. The recent progresses about this Fourier-related method are reviewed and discussed.


1.
Introduction. This paper deals with instabilities involving many spatial oscillations with a wavelength much smaller than the characteristic length of the considered domain. These problems are not new and nonlinear approaches have been presented since the late sixties for problems of fluid convection [33,38] or structural buckling [3] and alternative fields of physics and engineering [10,20,42]. Two typical experimental and numerical results are presented in Figure 1. These issues generated much interest recently, in particular for wrinkling in membranes [17,25] and in film-substrate systems [7,8,26,45]. All these phenomena can be modelized by boundary value problems whose set of solutions can be very intricate. As early established, there are many solutions that can be characterized by their wavenumber near the bifurcation point [9,35]. The response curves seem very intricate with a lot of bifurcations and limit points, what can be analyzed from various points of view: "dynamics" close to a heteroclinic orbit, numerical continuation techniques [23,41], envelope equations of Ginzburg-Landau type [9,11,24] ... According to [23], this apparently cumbersome behavior seems to follow an "ubiquitous scenario" occurring in the case of cellular instabilities with a 1D distribution of spatial oscillations. This point of view is assessed by its occurrence in many situations, as well in solid as in fluid mechanics. Moreover generic envelope equations have been established by asymptotic approaches for all these instability cases, what confirms the generic character of this scenario [11,24].
Here we are mainly interested in membranes that have a lot of applications, in spacecraft structures, civil engineering, biological tissues ... Despite of the aforementioned difficulties, classical shell finite element codes can permit to modelize in details membrane wrinkling and to describe the orientation and the amplitude of the wrinkles, but the computational cost is important and the management of the calculations can be difficult due to the large number of solutions. That is why many practical computations are based on simpler models that consider only the effect of wrinkling on membrane behavior by disregarding the description of the wrinkles [36].
In this paper, some multi-scale approaches are reviewed to describe concurrently the scale of the wrinkles (called microscopic) and the scale of the structure (called macroscopic). They are based on the concept of Fourier series with variable coefficients that has been introduced in [12] and has been developed and applied in many cases of instabilities in solid mechanics: beam on foundation, sandwich structures, 2D nonlinear elasticity, membrane wrinkling, film-substrate systems [4,13,14,15,16,21,22,27,28,29,31,44,45,46]. The goal is to define reducedorder models that can describe many oscillations with a relatively coarse finite element mesh. The deduced macroscopic models look like the envelope equation of Ginzburg-Landau, but contrarily to this asymptotic approach, the Partial Differential Equations are well-posed and can be associated with any domain shape and any boundary conditions.
2. About multi-scale approaches. Multiple scale methods are very classical in nonlinear science. The common idea is to distinguish two levels of variation, a local or microscopic level and a global or macroscopic level. For instance a nearly harmonic function can be described locally by a sine function and globally by two slowly variable functions: the mean value and the envelope of the oscillation, as pictured in Figure 2. The oldest approaches combine this two-level variation with asymptotic expansions, the main applications being in the field of nonlinear vibrations [6,32]. Typically an unknown function u(t, τ ) depends on a rapid (or local) variable t and a slow (or global) variable τ = εt. Moreover an asymptotic expansion with respect to the small parameter ε is performed. In this analysis, one starts from a model deemed to be a reference and the asymptotic process permits to deduce a sequence of differential equations alternatively at fine and coarse level. Homogenization of heterogeneous materials is another well known example of a multiple scale approach. In its famous asymptotic version [5,37], it mimics the double scale asymptotic used for nonlinear vibration by combining an asymptotic expansion of the unknown with the assumption of two independent levels of variation. For materials with a nonlinear constitutive law, the asymptotic expansion has to be given up and the nonlinear homogenization method introduced in [40] relies only on the framework of two independent levels. This leads to two nested continua, where the solutions at the local level provide an implicit constitutive law for the macroscopic level. These nested continuum models can be discretized in a classical way and this leads to the so called FE 2 technique [18,34,39]. Nowadays this two-level numerical procedure can be carried out with a personal computer in the 2D case, but this double discretization remains expensive. Such a heavy numerical tool may be not compulsory to account for cellular instabilities because the observed patterns seem to vary as a sine function. Thus it is proposed to represent the local behavior by few terms of a Fourier series. The method reviewed in this paper relies on this principle: it involves two independent levels of variations, but the local problem is represented by terms of a Fourier series so that one can get explicitly the macroscopic model to be solved with a much coarser mesh than the initial problem. Note that Fourier series have been successfully introduced [30] to solve the local problem for heterogeneous media with a small computational cost.
Despite of the common use of Fourier series, our approach is slightly different because only few terms will be sufficient in the case of multi-scale instabilities, which permits to define a-priori purely macroscopic models.
3. Fourier series with variable coefficients: Principle and main rules. Let us consider a 1D problem, where the unknown is a function U (x), x ∈ R, U ∈ R k . One assumes a nearly periodic behavior, as in Figure 2, which means that U (x) can be represented in the form of a truncated Fourier series: where the envelopes U n (x) vary more slowly than the harmonics e inqx . In other words, one keeps the assumption of a double level of variation as in the other multiscale approaches, the two levels being represented by the same symbol "x". The wavenumber q is given. Fully rigorous results about this decomposition (1) (uniqueness...) should be welcome. A discussion about the inverse transition U (x) → U n (x) can be found in [44]. The same idea can be extended in a multi-dimensional framework, where the wavenumber and the direction of the rapid oscillation is given [15]. The order of truncation is small (N = 1 or N = 2). As usual with Fourier series, the term of level 0 is real-valued while the other terms U n (x), n ≥ 1 are complex-valued. Nevertheless, to build up the most simple model, we can limit ourselves to real valued function U n (x), which leads to a model with only two envelope functions: the mean value and the envelope of the oscillation. In such a simplified framework, the amplitude of the oscillation can be modulated, but not its phase. Three main rules permit to deduce the differential equations of the macroscopic model from a given initial model. The first rule concerns the derivatives that are computed exactly in a first time. The n-th Fourier coefficient of a derivative can be computed as The second rule follows from the identification of Fourier series. According to the assumption of two independent levels, the envelopes are assumed to be constant in this identification. For instance, let us consider the nonlinear constitutive law of a beam in the Föppl-Von Karman approximation. This equation relates three unknowns, the beam tension N (x) and the two components of the displacement u(x), v(x), ES being a given constant: If one considers only the order n = 0 for the two first variables (mean value) and the order n = 0, ±1 for the transversal displacement v(x), the identification and the rule (2) allows to deduce a law for the envelopes, i.e. at the macroscopic level: One observes a coupling between the mean tension N 0 and the amplitude of the oscillations v 1 as it is expected in the cases of membrane wrinkling pictured in Figure 1. The last expression in (4) corresponds to the assumption v 1 ∈ R. This second rule is very easily applied in an energy framework or by starting from a weak form: indeed, because the envelopes are assumed constant on a period, the mean value of a harmonic term vanishes, except of course for n = 0. For instance with two harmonics u(x) = u 1 (x)e iqx +ū 1 (x)e −iqx , the Dirichlet integral can be rewritten in a reduced form as: The third rule concerns additional simplifications to avoid some spurious effects. These spurious effects can be understood by considering the eigenvalue problem given by the linearized Swift-Hohenberg equation [3], λ being a scalar parameter: Disregarding boundary conditions and considering only harmonic solutions, the smallest eigenvalue is λ 0 = 2 (minimum of the curve λ(q) = q 2 + 1/q 2 relating the parameter λ and the wavenumber q). The corresponding minimizer is q = 1 and the corresponding eigenmodes are e ix and e −ix . Next let us go back to general solutions of (6). According to the rule (2), the first envelope v 1 (x) should satisfy the following equation Because of the linearity of (6), the two equations (6) and (7) are equivalent via the change of variables v 1 (x) = v(x)e −ix so that the eigenmodes of the envelope equation (7) are v 1 (x) = 1, that is slowly varying, and v 1 (x) = e −2ix that is varying rapidly and therefore is not consistent with the previous assumptions. Thus the macroscopic model has some slowly varying solutions, as expected, but it has also other very oscillating solutions. There are several ways to disregard these unexpected oscillating solutions. A simple technique should be to discretize (7) with a coarse mesh. A finite element discretization of (7) by Hermite polynomials has been achieved in a case with 50 periods (the length of the domain is 100π) and with clamped boundary conditions (i.e. the unknown v(x) and its first derivative vanish at the ends of the interval), see [13]. A very accurate solution can be obtained with one or two elements as in Figure 3, but spurious oscillations are observed if the mesh is refined. A simple way to remove these oscillations is to drop the highest derivatives in the macroscopic equation (7). This leads to the second order equation that is nothing but the linearized Ginzburg-Landau equation classically obtained within the asymptotic multi-scale method. Within the latter, the highest derivatives are eliminated because their order of magnitude is smaller than the leading terms.
That is why it is consistent to drop the derivatives of order 3 and 4 in (7). The transition between microscopic and macroscopic descriptions is questionable. By starting from the macroscopic unknowns U n (x), it is easy to deduce the full description by the Fourier series formula (1). In the opposite sense, it seems natural  (6) obtained with a finite element solution of (8) and a single quadratic element. The same result can be obtained from the fourth order envelope equations (7) with a single element, but spurious oscillations appear by solving (7) with more elements.
to apply also a basic rule of Fourier series as: Nevertheless one is interested to nearly periodic functions while the formula (9) is exact only for periodic solutions. Numerical experiments have been done in [44] by considering numerical solutions of an ordinary differential equation. It was established that the formula (9) gives consistent envelopes, but with small unexpected oscillations, for the same reasons as previously. Fortunately, if one comes back to the initial function by the transition formula (1), the initial function U (x) is perfectly recovered. Thus the transition between the two levels seems rather easy from (1)(9), but to get an optimal transition micro → macro, a filter should be associated with (9) to erase the small and rapid oscillations. The question of boundary conditions is not simple within macroscopic envelope models, even with the oldest asymptotic multi-scale approaches. Indeed the transition from a fourth order EDO (6) to the second order EDO (8) leads to a loss of accuracy of the reduced model near the boundary, what generally implies the existence of boundary layers for the full model. Accounting for boundary layers may be considered, but this leads to tedious and generally non generic approaches as in [9]. For instance, let us look at a single function approximated by two harmonics n = ±1, v(x) = v 1 (x)e iqx +v 1 (x)e −iqx . Close to the end x = 0, the envelope can be considered as nearly constant: For a boundary conditions of the type "simply supported" v(0) = d 2 v dx 2 (0) = 0, one gets only a condition on the phase, but nothing for the amplitude |v 1 (0)|. On the contrary, for a "clamped" end v(0) = dv dx (0) = 0 and with an approximation [13,21], one gets v 1 (0) = 0, what is consistent with the numerical result in Figure 3. This condition of a vanishing amplitude seems reasonable due to the inconsistency between the approximate form and the boundary conditions, see for instance the experimental buckling results for long rectangular plates [43]. A more intricate, but efficient, numerical procedure has been proposed and applied to one-dimensional cases: one applies the chosen reduced model in the bulk, the full model near the boundary and the two numerical solutions are connected by a bridging technique [21,44]. 4. Few examples of reduced systems. The present multi-scale method does not keep two active levels of variation, as in FE 2 technique because the local variations are accounted by few terms of a Fourier series. Thanks to the properties of Fourier series, the local level can be disregarded and rather simple macroscopic models can be deduced explicitly from the full model by applying the three rules described in §3. Several reduced-order model have been obtained in this way and some of them are shortly described here. First we discuss a simple system of ODE's that represents a beam whose properties are denoted by the constants ES, EI. This beam is connected to a foundation of properties c, α, β. The equations and a corresponding envelope model are given on Table 1. The derivation of the reduced-order model is straightforward and follows the three rules described in §3, except for one point due to the presence of a quadratic term in the behavior of the foundation (α = 0): a second harmonic n = ±2 has been accounted in a simplified manner (see [14]) for a proper prediction of the first bifurcation, next this second harmonic has been dropped to get the model presented in Table 1, where three unknowns have been kept: the harmonic zero for the beam tension N 0 (x) and the axial displacement u 0 (x), the harmonic n = 1 for the transversal displacement v 1 (x) ∈ C. A more simplified model (v 1 ∈ R) or improved approximations are available in the literature. Note also that the last equation on the right column coincides with the famous Ginzburg-Landau if one fixes the beam tension at its critical level N 0 (x) = −2EIq 2 : this means that the envelope model is consistent with the asymptotic theory and permits to recover at least this initial bifurcation behavior. Nevertheless, the only assumption being the existence of two independent scales, it can be valid away from the bifurcation.

Full model Beam on foundation
Envelope model (2 harmonics) N, u: harmonic n = 0 v: harmonic n = 1 (complex) Table 1. Beam on foundation and a corresponding macroscopic model, from [14].
The second example is the famous Föppl-Von Karman plate system of equations that is the extension of the previous model in the 2D case. This plate model can describe buckling and wrinkling as in Figure 1, its main restriction being the moderate size of the rotation angles. It is recalled in Table 2, as well as the corresponding envelope model established in [15,16] with similar assumptions as the previous system of ODE's, the envelope of the transversal displacement w 1 (x, y) being real. The reduced system of Table 2 has been established for rapid oscillations in the direction of Ox. It looks like a nonlinear Poisson equation coupled with a nonlinear elastic plane model; nevertheless, the nonlinear part of the strain depends on the envelope w 1 (x, y) and not only on its gradient for the same reasons as in (4).  Table 2. Föppl-Von Karman plate and a corresponding macroscopic model, from [15].
The third example is nonlinear elasticity, with the so-called Saint-Venant Kirchhoff constitutive law that is a linear relation between the nonlinear strain tensor γ and the second Piola-Kirchhoff stress tensor s. The well known full system has the following form, where the balance equation is written in a weak form: The macroscopic models have been established in [13] and they have the same structure as the initial one (12). For instance the macroscopic stress is represented by an assembly of Fourier coefficients of the initial stress tensor s : The macroscopic strains and displacements are defined in the same way and the displacement gradient is obtained with the help of the first rule (2). The final macroscopic model is in the following form and we refer to [13] for more details: 5. Few numerical applications. One can easily solve the obtained reduced model, see (12) or the Tables 1 and 2. Even few analytical solutions have been found [16]. The presented computations have been obtained by classical finite element discretization.
The first numerical application is a sandwich beam made of one soft core and two stiff and thin symmetric skins. The structure is submitted to a pure compression, see Figure 4. In the presented case, the first instability is a global buckling of the structure, as for a homogeneous beam. This buckling tends to increase the compression in the top skin, what induces a local instability of this skin. Three calculations are presented in this figure. The first one is a full problem represented by the equations (11) for 2D nonlinear elasticity. This reference result is compared with two reduced models: the one is based on a reduction from 2D elasticity as in (12), the other one looks like the reduced model of Table 1, but it was built from a more advanced 1D system of ODE's designed to describe such instabilities in sandwich structures, see [27,46]. The two reduced models permit to obtain accurately the two successive bifurcations. After the second bifurcation, the response of the reduced model (12) remains very accurate, while the second reduced model diverges from the reference: clearly this is due to a weakness of the starting model and not of the Fourier reduction. This result shows the ability of this Fourier method to account for global instabilities (via the harmonics n = 0) coupled with wrinkling instabilities via the harmonic n = ±1. There are two bifurcation points, the second one corresponding to the appearance of wrinkles in the more compressed skin. Results from [27,28].
The second application is a classical traction problem of a rectangular membrane. Due to boundary conditions on the sides parallel to OX (especially u X = 0), a small compressive zone appears not too far from these sides, what induces wrinkles with the shape pictured in Figure 5-b. This problem has been analyzed by the reduced plate model of Table 2 and compared with well established shell finite elements. The reduced model predicts fairly accurately the deformed shape ( Figure 5-c), even in this case where the number of periods is small. Let us recall that this reduced model with a real envelope locks the period, while the full model predicts an increase of the wavelength when one goes from the center to the sides. However this very simplified model gives correct results, even if a model with a complex envelope (w 1 ∈ C) should be more accurate to account for an evolution of the wavelength. The interest of using a complex envelope is corroborated by a detailed numerical study [29] of the 1D system of Table 1.  The last application concerns a film/substrate system, i.e. a very thin and stiff film that relies on a thick and soft substrate, as in Figure 6. There is a wide and recent literature on the topics and various types of cellular instabilities can occur under the effect of strain mismatch or of loads applied on the boundary of the domain. It is possible to discretize the boundary value problem by the finite element method, but the computational cost increases dramatically if the studied domain includes many spatial periods. We applied a reduced system, the film being represented by the equations of Table 2 and the substrate by an envelope model deduced from linear 3D elasticity. Both were discretized by standard P1 or P2 finite elements. The intricate pattern of Figure 6 was obtained with only four elements along the length (respectively 20 along the width and 5 layers for the substrate). Roughly, with respect to a full model, the number of degrees of freedom is divided by 10 and the computation time by a ratio 50 ∼ 100. Similar performances have been obtained for some numerical simulations of membrane wrinkling [22].
6. Discussion. The method of Fourier series with slowly variable coefficients has been overviewed and discussed. It permits to define explicitly a family of macroscopic Partial Differential Equations, with a strong reduction of the computation cost. It yields rather simple macroscopic systems of Partial Differential Equations whose discretization leads to very interesting reduced-order numerical models. Up to now it has been applied to buckling and wrinkling of thin membranes and thin films on soft foundations, but applications in fluid mechanics and in other fields could be considered.
The weak point of the procedure is the need to choose a priori the period and the direction of the oscillations. Sometimes several computations with several wavenumbers have been performed with the hope that the best one corresponds to the lowest bifurcation, but this is not sufficient for the numerous cases where the local period and the direction of the wrinkles vary in the domain, for instance as for loss of flatness after a sheet rolling process, see [2] and Figure 1. The reduction technique should been improved to account for the evolution of the wave vector. Finding the local wavelength is difficult, but it is intrinsically related to the description of multi-scale instabilities and not to the Fourier reduction method: even if the material properties are locally periodic, possible instability wavelengths do not coincide always with the material wavelength, as explained in [19].
Finally the method is not based on solid mathematical results and some works in this direction will be welcome