## Abstract

A wide variety of surface morphologies can be formed by compressing a bilayer comprising a thin film bonded to a compliant substrate. In particular, as the applied strain is increased, secondary instabilities are triggered and the initial sinusoidal wrinkles evolve into new complex patterns. Here, we propose a robust numerical analysis based on Floquet–Bloch boundary conditions to detect the primary and secondary instabilities triggered upon compression. Because the proposed method is based on unit cell simulations, it is computationally very efficient. Moreover, it accurately predicts not only the critical strains, but also the corresponding critical modes and their wavelengths, enabling us to follow the evolution of the surface morphology as the applied strain is progressively increased.

## 1. Introduction

Bilayers comprising a thin film bonded to a compliant substrate wrinkle under compression [1–5]. Interestingly, this observation has not only motivated the formulation of theories to predict both the critical strain and the wavelength of the wrinkles [6–9], but has also been instrumental for a wide range of applications. These include the design of tunable flexible electronics [10–13], novel micro/nanofabrication techniques [14,15], optical tunable devices [16,17] and adhesives [18,19], and the development of new approaches to measure the mechanical properties of thin films [20,21].

More recently, several experiments have revealed that for large values of applied compressive strain the wrinkles are no longer stable and new bifurcations occur, accompanied by changes of the wrinkled pattern. In particular, it has been shown that a second bifurcation can result either in period doubling [22,23] and tripling [23] or in the formation of chaotic patterns with mixed periodicities [24] and localized ridges [25–27]. Moreover, additional instabilities have been triggered by further increasing the applied deformation, leading to period quadrupling [22,23].

To understand these new bifurcations, both analytical and numerical approaches have been developed. On the analytical side, linearized stability analyses superimposed on the deformed wrinkled configuration have been used [28,29]. However, because the secondary bifurcation are triggered at high values of applied compressive strain when both geometric and material nonlinearities play an important role, these analytical approaches either result in limited accuracy or require to be complemented by numerical calculations [28].

On the numerical side, the nonlinear finite-element (FE) method employing both finite size [23,29,30] and unit cell [24,31] models has been successfully applied to capture the patterns induced by these secondary bifurcations. Although within the FE framework the stability of a structure can be easily examined using eigenvalue analyses [32], it is important to note that the wavelength of the pattern induced by the instability in a infinitely long bilayer system is a priori unknown. Therefore, the stability of models of different lengths capable of accommodating wrinkled patterns of any wavelength should be investigated. The mode associated with the minimum of the critical strains on all possible models should be then selected as the critical one. However, given the significant computational cost associated with this rigorous procedure, in recent studies the formation of complex patterns has been investigated either by comparing the elastic energies of modes with a limited set of potential periodicities [24,31] or by introducing an artificial damping coefficient and simulating directly the post-buckling response [23,29,30]. A rigorous and efficient numerical approach to fully track the sequential bifurcations observed during compression of an infinitely long bilayer is still lacking.

Here, we present a robust and efficient numerical approach to accurately determine all instabilities triggered in an infinitely long bilayer system upon compression. Our method is based on the Floquet–Bloch boundary conditions [33], which have been recently proved valuable to investigate the changes in periodicity induced by instabilities in infinite periodic structures [34–38]. In fact, it has been shown that in such structures buckling can be explored considering a given unit cell and investigating the propagation of small-amplitude perturbations of arbitrary wavelength superimposed on the current state of deformation. While a real natural frequency corresponds to a propagating wave, a complex natural frequency identifies a perturbation exponentially growing with time. Therefore, the transition between a stable and an unstable configuration is identified when the frequency vanishes and the periodicity of the associated mode is provided by the wavelength of the perturbation. While in previous studies this procedure has been successfully used to identify the first critical bifurcation, here we also extend it to capture the secondary instabilities, providing a robust tool to accurately investigate the formation of complex morphologies in bilayers.

The paper is organized as follows. After presenting the geometry (§2a), constitutive material model (§2b) and loading conditions (§2c) considered in this study, in §2d,e, we describe the numerical analysis that is used to investigate the primary and secondary bifurcations triggered upon compression. Finally, numerical results are presented and discussed in §3, highlighting the effect of the prestretch applied to the substrate on the formation of complex patterns.

## 2. Modelling

To study the response under uniaxial compression of a bilayer comprising a thin film bonded to a compliant substrate, we perform nonlinear simulations on unit cell models using the FE package ABAQUS/Standard. All models are constructed using eight-node quadratic quadrilateral plane strain elements with reduced integration (CPE8R element type). Each mesh is most dense towards the top (where the thin film is bonded) and gradually coarsens as we move towards the bottom, with approximately 9000 elements per unit cell.

Here, after briefly describing the material model and loading conditions, we present the Bloch wave analysis that we perform on a unit cell to investigate the stability of the structure. Interestingly, this analysis allows us to capture the strains and modes associated not only with the critical instability, but also with the subsequent secondary bifurcations.

### (a) Geometry

In this study, we consider a bilayer structure under plane strain conditions, as that shown in figure 1. The system comprises an infinitely long, thin and stiff film of thickness *h*_{f} perfectly bonded to a prestretched compliant substrate. In particular, we assume that a prestretch λ_{p} is applied to the substrate along the direction of the thin film prior to bonding. Finally, although in the analytical calculations developed to predict the onset of the first critical bifurcation the substrate is assumed to be infinitely thick [8], owing to computational limitations we set *h*_{s}=4600*h*_{f}, *h*_{s} denoting the substrate thickness. Note that this choice for *h*_{s} guarantees that the substrate thickness is much greater (approx. 30 times) than the expected wrinkling wavelength, so that the infinite substrate assumption is satisfied.

### (b) Material constitutive behaviour

In this study, both the substrate and the film are modelled as hyperelastic materials. In particular, their stress–strain response is captured using a nearly incompressible neo-Hookean model [39], whose strain energy is given by
*μ* and *K* are the initial shear and bulk moduli, respectively. Moreover, *I*_{1}=tr **F**^{T}**F** and **F**=∂**x**/∂**X** denotes the deformation gradient—a linear transformation which maps a material point from its reference position ** X** to its current location

**x**. For such material, the Cauchy stress,

*σ*, is given by

**I**is the identity tensor and

**B**=

**FF**

^{T}.

In all our analyses, we use *μ*_{s}=0.15 MPa and *μ*_{f}=6.67 GPa as shear moduli of the substrate and film, respectively. Moreover, we assume *K*_{s}=1000 *μ*_{s} and *K*_{f}=1000 *μ*_{f}, corresponding to Poisson's ratios *ν*_{s}=*ν*_{f}=0.4995. The resulting Young's moduli are *E*_{s}=2*μ*_{s}(1+*ν*_{s})=0.45 MPa and *E*_{f}=2*μ*_{f}(1+*ν*_{f})=20 GPa.

Finally, to describe the effect of prestretch in the substrate, we decompose the deformation gradient into a loading component, **F**^{l}, and a prestretch component, **F**^{p}, [40,41]

### (c) Loading conditions

A plane–strain deformation in the *X*–*Y* plane is considered, where the bilayer is compressed in the horizontal direction (i.e. in the direction of the thin film) by applying a strain *ϵ*. To subject the unit cells used in this study to such deformation, we define periodic boundary conditions on their two vertical edges (because the structure is periodic only along the horizontal direction). Focusing on pairs of nodes periodically located on the left and right edges, we relate their displacements (*u*_{x} and *u*_{y}) as [37,42]
*r* and *l* denote quantities associated with nodes on the right and left edges, respectively. Moreover, we assume traction free boundary conditions on the top and bottom edges. Note that equations (2.4) are implemented into ABAQUS/Standard via linear multi-point constraints and the virtual node approach [37,42].

### (d) Instability analysis: Bloch wave analysis

Although the wavelength of the wrinkled pattern triggered at instability is *a priori* unknown, buckling can still be rigorously studied on a unit cell of fixed length *L* by investigating the propagation of small-amplitude waves of arbitrary wavelength superimposed on the current state of deformation [37,43]. Because a wave motion can be described as *ω* correspond to propagating (stable) perturbations, whereas waves with imaginary natural frequencies *ω* correspond to exponentially growing (unstable) perturbations. Therefore, the onset of instability is identified when *ω* changes from real to imaginary (i.e. when a wave with vanishing frequency, *ω*=0, is detected). The periodicity and deformation of the associated mode is then provided by the wavelength and Bloch mode of the perturbation.

To investigate the stability of the considered bilayer, we first finitely deform the unit cell by applying a strain *ϵ* and then investigate the propagation of small amplitude elastic waves in the deformed configuration. To this end, we apply Bloch-periodic boundary conditions [33] on the two vertical edges of the unit cell,
*L* is the undeformed length of the unit cell and *k* is the magnitude of the wavevector **k**=*k***e**_{x}. In our analysis, we consider a large number of **k** vectors in the reciprocal lattice (i.e. *k*∈[0,2*π*/*L*]) and calculate the frequency *ω* of the waves associated with each **k** using a linear perturbation procedure (module *FREQUENCY in ABAQUS). An instability is detected at the smallest magnitude of compressive strain, *ϵ*_{cr}, where we find a wavevector **k**_{cr}=*k*_{cr}**e**_{x} associated with a wave with *ω*=0. The periodicity of the pattern induced by the instability is then obtained as

Finally, we note that to work with the complex-valued relations of the Bloch-periodic conditions (equation (2.5)) in a commercial software such as ABAQUS/Standard, all fields are split into real and imaginary parts and two identical FE meshes are used for the unit cell [37,44].

### (e) Post-buckling analysis: formation of complex patterns

Conclusively, to investigate the formation of complex wrinkled patterns in the bilayer system, we use the Bloch wave analysis in combination with post-buckling analysis. In particular, as shown in figure 1*b*,

(i) we perform the Bloch wave analysis on a unit cell of arbitrary length

*L*to detect the critical strain,*ϵ*_{cr}, and the associated mode with periodicity*L*_{cr};(ii) we create a new unit cell of length

*L*_{cr}and introduce an imperfection in the form of the critical mode;(iii) we perform a nonlinear static analysis and apply a strain

*ϵ*<*ϵ*_{cr}to capture the formation of the wrinkled pattern induced by the instability;(iv) we repeat the Bloch wave analysis on the wrinkled unit cell of length

*L*_{cr}to detect the critical strain,*ϵ*_{cr,2}, and mode with periodicity*L*_{cr,2}associated with the second bifurcation;(v) we create a new unit cell of length

*L*_{cr,2}and introduce an imperfection in the form of a combination of the critical modes associated with the first and second bifurcation;(vi) we perform a nonlinear static analyses and apply a strain

*ϵ*<*ϵ*_{cr,2}to capture the formation of the complex wrinkled pattern induced by the second bifurcation; and(vii) we repeat the same procedure to detect additional bifurcations encountered during loading.

## 3. Results

Here, we present numerical results obtained from the unit cell analysis discussed in §2d,e. To validate the proposed method, we first compare the numerically calculated critical strain and wrinkle wavelength associated with the first bifurcation to analytical predictions. Then, we show that our numerical analysis correctly captures the different types of secondary bifurcations that have been observed when the prestretch applied to the substrate is varied.

### (a) First bifurcation

It is well known that sinusoidal wrinkles form almost immediately when a compressive strain is applied to a bilayer and that the magnitude of the critical strain, |*ϵ*_{cr}|, decreases as the stiffness contrast between the film and substrate increases. In particular, for the case of a neo-Hookean substrate subjected to a prestretch λ_{p} and bonded to an elastic beam, it has been shown that *ϵ*_{cr} is given by [23]^{1}
*L*_{cr}, can be obtained by solving [23]
*Q*=*μ*_{s}(1−*ν*_{f})/(2*μ*_{f}) is the scaled stiffness ratio. Note that in the limit case of no prestretch applied to the substrate (i.e.

While in previous numerical analysis, these analytical results are taken as input and used to determine the size of the FE model [30], here we use them to validate the stability analysis discussed in §2d. We start by considering a stress-free substrate (i.e. λ_{p}=1) and choose the length of the unit cell to be one-fourth of the predicted wavelength (i.e. *L*=0.25*L*_{cr}). We then progressively increase the applied compression and after each strain increment Δ*ϵ*=−2.0×10^{−4} we check the propagation of waves characterized by a wavevector ranging from *k*=2*π*/(20*L*) to *k*=2*π*/*L* with an increment of Δ*k*=2*π*/(20*L*).

The results of the Bloch wave analyses are shown in figure 2*a*. In the plot, we report the evolution of the squared frequency, *ω*^{2}, as a function of the applied strain *ϵ* for seven different wavevectors. In the undeformed configuration (i.e. *ϵ*=0), all frequencies are positive. However, as *ϵ* increases in magnitude, the frequencies associated with each wavevector **k** gradually decrease and eventually become negative. The critical strain parameter associated with each **k** can be easily extracted from the plot, because it corresponds to the intersection point between each curve and the horizontal line (*ω*=0). Finally, the onset of instability for the bilayer is defined as the minimum critical strain on all the considered **k**. As expected, the frequency associated with *k*=2*π*/(4*L*) vanishes first at *ϵ*_{cr}=−4.14×10^{−4}, yielding a perfect match with the analytical predictions (equation (3.3)). Moreover, the critical mode can be reconstructed from the corresponding Bloch mode. As shown in figure 2*b*, the wrinkled pattern is characterized by a perfect sinusoidal profile, as predicted analytically.

Next, we repeat the same analysis on bilayers in which the substrate is subjected to different amounts of prestretch λ_{p} and compare the numerical results with the analytical solutions. The results reported in figure 2*c*,*d* show that the numerical analysis correctly capture the evolution of both the critical strain, *ϵ*_{cr}, and wavelength, *L*_{cr}, as a function of λ_{p}, further validating the proposed numerical analysis.

### (b) Second bifurcation

If the compressive strain is further increased after the formation of sinusoidal wrinkles, new instabilities can be triggered, resulting in significant changes of the surface morphology [22,23,25–27]. It has been recently shown that these post-wrinkling bifurcations are highly affected by the prestretch applied to the substrate [24,45]. In particular, (i) for 1.4≥λ_{p}≥0.7, the second bifurcation has been found to be associated with period doubling and to be followed by a third instability resulting in period quadrupling [22,23]; (ii) for λ_{p}≤0.7 (large precompression) the formation of chaotic patterns with mixed periodicities has been observed [24]; (iii) finally, for λ_{p}≥1.4 (large pretension) it has been shown that localized ridges emerges [25–27].

Here, we first focus on a stress-free substrate (i.e. λ_{p}=1.0) and show that in this case our numerical analysis correctly captures the second bifurcation that results in period doubling. Moreover, our stability analysis also detects a third bifurcation accompanied by period-quadrupling when the applied compression is further increased. Then, we investigate the response of a bilayer in which the substrate is largely precompressed (i.e. λ_{p}≈0.6). In this case, our numerical analysis indicates that several modes are triggered almost simultaneously. As a result, they interact upon loading and lead to the formation of chaotic patterns. Finally, we study the response of a bilayer in which the substrate is largely pretensioned (i.e. λ_{p}=2.0) and find that the secondary instability results in the formation of ridges.

#### (i) Stress-free substrate: period doubling and quadrupling

To investigate the response under large values of applied compression of a bilayer system with a stress-free substrate (i.e. λ_{p}=1), we use a unit cell of length *L*_{cr} (*L*_{cr}=154*h*_{f} being the wavelength predicted by the first bifurcation analysis). We first introduce a geometric imperfection into the mesh in the form of the critical sinusoidal mode with an amplitude of 0.05*h*_{f} and then compress the system by applying a uniaxial strain *ϵ*. As expected, for *ϵ*<*ϵ*_{cr}=−4.14×10^{−4} a sinusoidal wrinkled pattern emerges and becomes more accentuated as the magnitude of the applied strain is increased.

For *ϵ*<*ϵ*_{cr}, we monitor the stability of the wrinkled pattern by investigating the propagation of small-amplitude waves superimposed to the deformed configuration after each strain increment Δ*ϵ*=−0.01. Note that in this case we consider wavevectors *k*=2*π*/(m L_{cr}), where *m* is an integer number that varies from 2 to 20. In fact, while for the first bifurcation any wavelength is possible (so that, in principle, any wavevector within the reciprocal lattice should be considered), for this second instability the periodicity of the new pattern should be a multiple of the wavelength of the sinusoidal mode, *L*_{cr}.

The results of this stability analysis are reported in figure 3*a*, where only the frequencies associated with five wavevectors are shown for the sake of clarity. The numerical data clearly show that at *ϵ*=*ϵ*_{cr,2}=−0.187 the mode resulting in period doubling (i.e. *k*_{cr,2}=2*π*/(2*L*_{cr})) becomes unstable, in good agreement with previous studies [24,29]. The mode triggered by this second bifurcation is then reconstructed from the corresponding Bloch mode and is reported in figure 3*b*, showing that neighbouring valleys alternatively sink and rise.

To capture the formation of this complex pattern, we then use a unit cell with length of *L*_{cr,2}=2*L*_{cr} and introduce a geometric imperfection in the form of the mode shown in figure 3*b* with an amplitude of 0.05*h*_{f}. The model is then uniaxially compressed and snapshots of its deformed configurations at different levels of applied strain *ϵ* are shown in figure 3*c*. Clearly, for *ϵ*>−0.187, the bilayer is characterized by a sinusoidal wrinkled pattern. However, as *ϵ* is decreased below *ϵ*_{cr,2} one of the valleys rises, whereas the other becomes more pronounced, resulting in a pattern similar to that shown in figure 3*b*.

Again, for *ϵ*<*ϵ*_{cr,2}, we monitor the stability of this new, more complex surface pattern after each strain increment Δ*ϵ*=−0.01 by investigating the propagation of small amplitude waves with *k*=2*π*/(mL_{cr,2}) (*m*∈[2,20]) superimposed on the deformed configuration. The results reported in figure 4*a* indicate that at *ϵ*=*ϵ*_{cr,3}=−0.256 a wave with vanishing frequency exists for *k*=2*π*/(2*L*_{cr,2}). Therefore, at *ϵ*_{cr,3}, a third instability is triggered resulting again in period doubling (period quadrupling with respect to the initial sinusoidal wrinkled pattern). The reconstructed mode shown in figure 4*b* shows that neighbouring valleys alternatively sink and rise, resulting in the formation of pronounced folds. Finally, we use a unit cell with *L*_{cr,3}=4*L*_{cr} to investigate the post-buckling deformation. Snapshots reported in figure 4*c* show how the deformation evolves up to a compressive strain *ϵ*=−0.3. As predicted by the stability analysis, for *ϵ*<*ϵ*_{cr,3} the wavelength of the pattern becomes 4*L*_{cr} and pronounced folds emerge. Note that similar deformed configurations have been observed in experiments when large compressive strains are applied [22,23].

The above-reported post-buckling results are obtained performing unit cell simulations with periodic boundary conditions. To further validate the results, we also conduct dynamic implicit simulations using ABAQUS/Standard. For this set of simulations, we use a finite size model of length 12*L*_{cr} and, to facilitate convergence, we introduce some artificial, numerical damping (by setting the parameter *α* in the Hilber–Hughes–Taylor time integration algorithm equal to −0.1). Moreover, quasi-static conditions are ensured by monitoring the kinetic energy, and no geometric imperfection is introduced (because the numerical imperfection introduced by the integration algorithm is enough to trigger the formation of wrinkles). Considering the tremendous computational cost, we use four-node linear quadrilateral plane strain elements with reduced integration (CPE4R element type) and reduce the mesh density, resulting in a total of 13 000 elements. Finally, in these simulations, we uniformly displace horizontally all the nodes on the right edge, while constraining the horizontal displacement of the left edge (note that the vertical displacement of all the nodes on both the left and right edges is unset).

In figure 5, we report three snapshots from the dynamic implicit simulation at *ϵ*=−0.08, −0.21 and −0.271. The results perfectly agree with the predictions of our stability analysis and confirm that for bilayers comprising moderately prestretched substrates the surface morphology evolves from a perfectly sinusoidal pattern (see snapshot at *ϵ*=−0.08) to a pattern with double (see snapshot at *ϵ*=−0.21) and then quadruple (see snapshot at *ϵ*=−0.271) wavelength.

#### (ii) Largely precompressed substrate: chaotic patterns

Next, we investigate the response of a bilayer in which the substrate is largely precompressed and assume λ_{p}=0.63. As for the case of the stress-free substrate, we choose a unit cell of length *L*_{cr} (in this case *L*_{cr}=174*h*_{f}) and monitor the stability of the sinusoidal wrinkled pattern that emerges upon compression by conducting the Bloch wave analysis. The results reported in figure 6*a* show that at *ϵ*=*ϵ*_{cr,2}=−0.0125 a second bifurcation is triggered, this time resulting in period tripling (because *k*_{cr,2}=2*π*/(3*L*_{cr})). However, differently for the case of the stress free substrate where the critical strain associated with higher modes are well separated, the modes resulting in period tripling, doubling and quadrupling have very close critical strains (−0.0125, −0.013 and −0.013, respectively). Therefore, we expect that these modes interact upon loading, resulting in the formation of a chaotic pattern.

As a result, although the snapshots reported in figure 6*c* show that for a unit cell of length 3*L*_{cr} with periodic boundary conditions period tripling occurs when *ϵ*<*ϵ*_{cr,2}, we do not expect to observe such periodic pattern in a long bilayer. To substantiate this point, we perform dynamic implicit simulations on a finite size model of length 12*L*_{cr} (note that these simulations are performed exactly as described in §3b(i)). The snapshots reported in figure 7 confirm that a chaotic pattern emerges when the bilayer is largely compressed. In fact, while at *ϵ*=−0.011, a sinusoidal pattern is observed, at *ϵ*=−0.013, the periodicity is lost and we do not see any sign of period tripling.

Finally, we perform the Bloch wave analysis also for bilayers characterized by λ_{p}=0.61 and λ_{p}=0.59. The results shown in figure 8 are similar to those reported in figure 6, except that the critical mode in this cases results in period quadrupling and quintupling, respectively. However, in both cases, the critical strains associated with the first few modes are still very close, so that we expect again the formation of chaotic patterns when the system is largely compressed.

#### (iii) Largely pretensioned substrate: ridges

Finally, we study the response under uniaxial compression of a bilayer comprising a largely pretensioned substrate (λ_{p}=2.0). In this case, sinusoidal wrinkles with wavelength *L*_{cr}=114*h*_{f} start to form at *ϵ*_{cr}=−7.6×10^{−4} and we investigate their stability by conducting a Bloch wave analysis. The results reported in figure 9*a* indicate that a second bifurcation is triggered at *ϵ*_{cr,2}=−0.0264. The corresponding mode shown in figure 9*b* has wavelength *L*_{cr,2}=7*L*_{cr} and is characterized by the flattening of several wrinkles and the rising of a few ridges.

We then investigate how the surface morphology evolves for *ϵ*<*ϵ*_{cr,2} by uniaxially compressing a unit cell of length 7*L*_{cr} with an imperfection in the form of the critical mode shown in figure 9*b*. Because it has been recently shown that the ridges are the result of a snap-through instability [28,45], for this nonlinear static analysis, we use the pseudo-dynamic method and introduce a viscous force in the equilibrium equations (such force is defined as a damping factor, *δ*=0.0002, times the nodal velocity times an artificial mass matrix with unit density) [45]. Note that such velocity-dependent damping provides numerical stabilization when loss of stability occurs and allows the unstable transition from wrinkles to ridges to take place (differently from the case studied in previous sections, without this damping term, the simulations are found to stop at the onset of the second instability).

In figure 9*c*, we report numerical snapshots of the deformed bilayer at different levels of applied strain. The amplitude of the sinusoidal wrinkles is found to grow up to a strain of *ϵ*=−0.027. At that point, the wrinkles become unstable and a localized ridge abruptly forms, while the neighbouring wrinkles almost vanish. Note that the pattern emerging from the nonlinear static analysis comprises a single ridge and is different from that predicted by the stability analysis (figure 9*b*). This is because the transition from wrinkles to ridges is unstable.

Finally, we conduct an implicit dynamic simulation and study the response under compression of a finite size model of length *L*=12*L*_{cr}. The results reported in figure 10 indicate that also in this case an abrupt transition from wrinkles to a single ridge occurs. Note that the strain at which it occurs (*ϵ*=−0.0251) is higher than that predicted by the nonlinear static analysis (*ϵ*=−0.027). We attribute this discrepancy in the critical strain to the damping introduced in the simulations [45].

## 4. Conclusion

In summary, we have proposed a robust numerical approach based on the Bloch wave analysis to capture the successive bifurcations triggered in an infinitely long bilayer structure upon compression. We have shown that using computationally efficient unit cell simulations the critical loads, modes and wavelength associated with the primary and secondary bifurcations can be detected and a deeper understanding of the highly nonlinear response of the system can be gained. First, we have validated the proposed approach by comparing the results numerically obtained for the first bifurcation with available analytical solutions. Then, we have used the analysis to investigate the effect of the prestretch applied to the substrate on the secondary bifurcations. For the case of a stress-free substrate, we have found that upon compression two secondary instabilities are triggered. In both cases, the eigenvalues associated with the different modes are well separated, so that the one with the lowest eigenvalue invariably grows, resulting in period doubling followed by period quadrupling. Differently, for largely precompressed substrates, our analysis indicates that the secondary bifurcation is characterized by the appearance of closely spaced modes. This observation explains the formation of chaotic patterns that have been experimentally observed [24]. In fact, we expect the closely spaced modes to interact during loading, resulting in a pattern without a well-defined periodicity. Finally, for largely pretensioned substrates we find that the second instability results in the formation of ridges. Different from the other cases, the transition from wrinkles to ridges is unstable, resulting in an abrupt change of the surface morphology.

Compared with standard nonlinear FE simulations performed on finite-size models, our methods reveal new aspects of the nonlinear response of bilayers upon compression, facilitating the design of system with target behaviour. Interestingly, although the proposed method assumes that the pattern triggered at the onset of the instability is periodic, the emergence of irregular and chaotic surface morphologies can be still detected when closely spaced modes appear (but classical FE algorithms such as the dynamic implicit method are then required to capture the details of the emerging non-periodic pattern). Moreover, it is also computationally more efficient than the approach typically used to numerically investigate the stability of infinitely long bilayers. While traditionally a large number of models of various size are constructed and the critical strain of the infinitely long structure is then defined as the minimum of the critical strains on all possible models (obtained using a linear perturbation procedure), here we focus on a single, minimal unit cell, significantly reducing the number of degrees of freedom (note that the computational time is approximately cubic on the number of degrees of freedom). We also note that in this study to detect an instability we progressively increase the applied compressive strain by a fixed amount and then perform the Bloch wave analysis on the deformed configuration. While this approach is straightforward to implement, its convergence is slow. However, if only the critical bifurcation point is of interest, an efficient root-finding algorithm (such as the bisection or Newton's methods) can be used to more quickly detect the critical point. Finally, although in this paper we consider the simple case of a two-dimensional bilayer structure with both layers made of a neo-Hookean material, the proposed approach can be extended to study more challenging three-dimensional geometries as well as different material models.

## Ethics

The authors declare no misconducts of ethics.

## Data accessibility

All ABAQUS python files and user defined subroutines used to generate the data included in the manuscript can be downloaded from http://bertoldi.seas.harvard.edu/downloads. The access to the commercial FE package ABAQUS can be obtained through Dassault Systèmes.

## Authors' contributions

J.L. performed all the simulations and drafted the manuscript. K.B. coordinated and supervised the research and revised the manuscript. All authors gave final approval for publication.

## Competing interests

The authors declare no competing interests.

## Funding

This work was supported by the Materials Research Science and Engineering Center under NSF award no. DMR-1420570. K.B. also acknowledges support from the National Science Foundation (CMMI-1149456-CAREER).

## Acknowledgements

The authors thank Prof. John Hutchinson for helpful discussions.

## Footnotes

- Received July 23, 2015.
- Accepted September 21, 2015.

- © 2015 The Author(s)