## Abstract

We present an analysis of ruga-formation instabilities arising in a graded stiffness boundary layer of a neo-Hookean half space, caused by lateral plane-strain compression. In this study, we represent the boundary layer by a stiffness distribution exponentially decaying from a surface value *Q*_{0} to a bulk value *Q*_{B} with a decay length of 1/*a*. Then, the normalized perturbation wavenumber, *ε*, control formation of a wrinkle pattern and its evolution towards crease or fold patterns for every stiffness ratio *η*=*Q*_{B}/*Q*_{0}. Our first-order instability analysis reveals that the boundary layer exhibits self-selectivity of the critical wavenumber for nearly the entire range of 0<*η*<1, except for the slab (*η*=0) and homogeneous half-space (*η*=1) limits. Our second-order analysis supplemented by finite-element analysis further uncovers various instability-order-dependent bifurcations, from stable wrinkling of the first order to creasing of the infinite-order cascade instability, which construct diverse ruga phases in the three-dimensional parameter space of

## 1. Introduction

Material systems with tunable surface morphology are attracting attention recently owing to possibilities of constructing multifunctional materials with controllable effective surface properties of wetting [1–5], adhesion [4–7] as well as optics [8,9]. The common key ingredient of the material systems is the ability to control the surface patterns and the transitions among them by subjecting the material to a stimulus such as mechanical loading. A free surface or a layered structure of a soft material (e.g. polydimethylsiloxane (PDMS)), when compressed beyond a critical strain, becomes structurally unstable and may develop different modes of instability to create wrinkle, crease, ridge and fold patterns which are collectively denoted *Ruga* states [10]. Evolution of thin-film wrinkles on a soft substrate towards other ruga structures under increasing compressive strain has been recently studied by Song *et al*. [11], Brau *et al*. [12], Sun *et al*. [13] and Hutchinson [14]. A thick rubber slab with its modulus varying with distance from the free surface offers another example of using the elastic instabilities in constructing systems with different surface patterns [10,15].

Biot [16] showed, in his first-order perturbation analysis (named as incremental elasticity), that a homogeneous incompressible neo-Hookean material half space under plane-strain lateral compression becomes unstable at a critical strain of 0.456 for arbitrary wavelengths of sinusoidal surface-displacement perturbation, regardless of its elastic modulus. Later, Nowinski [17] studied the same problem using the formulation developed by Green *et al*. [18] for overall incompressible, isotropic hyperelastic materials. He found that, in general, the critical compressive strain for the instability is independent of the wavelength of the surface disturbance; however, it generally depends on the material elastic modulus in contrast to the behaviour of the neo-Hookean material. Biot's original work and Nowinski's later work were based on linear stability analysis. Biot's linear stability analysis predicted the critical compressive strain of 0.456 which turned out to be apparently inconsistent with the experimental observations of Gent & Cho [19]. They observed in their experiments that the free surface of a rubber close to a neo-Hookean material creased at critical strains in the range between 35%±7% which is below Biot's prediction. A clue to explain this discrepancy could not be found until recently when Hohlfeld & Mahadaven [20,21] conducted numerical and experimental stability studies on a thick slab under pure bending. They showed that the neo-Hookean material possesses a limit compressive strain (approx. 0.35) of subcritical crease instability below the Biot critical strain of 0.456. Hong *et al*. [22] also obtained this critical strain using a pinching perturbation scheme in their finite-element analysis (FEA). Cao & Hutchinson [23] showed, using Koiter's approach [24] assisted by finite-element analyses, that creasing instability of neo-Hookean half-space is imperfection sensitive and the critical strain of creasing may drop down to 0.27 for a sinusoidal imperfection amplitude.

Most of wrinkle or crease studies have been focused on a single- or multi-layer system with piecewise uniform material properties [23,25]. However, materials with graded stiffness such as human skins and geological structures are very common in nature. Rubber materials exposed to UV/ozone treatment offer another example of materials with graded mechanical properties that can be manipulated to design surfaces with predefined patterns. Recently, Guvendiren *et al*. [26] investigated wrinkle-to-crease transition in poly (2-hydroxyethyl) methacrylate gels with graded cross-linking density, and observed a dependency of the critical strains of wrinkling or creasing on the cross-link gradient. More recently, Torres *et al*. (J Torres, M Diab, M Monn, K-S Kim 2014, unpublished manuscript) have measured the modulus variation of a UV/ozone-treated PDMS along its depth with atomic force microscopy, and found a graded stiffness boundary layer with its elastic modulus decaying towards its bulk value along its depth from the free surface. Biot [15] performed linear bifurcation analysis on stability of a compressed neo-Hookean half-space with *exponentially decaying modulus vanishing at infinity* (vai-EDM) as well, and obtained a relation between the perturbation wavelength and the corresponding critical strain of instability. Later, Lee *et al*. [27] presented an analytical framework to determine the critical strain and the corresponding wavenumber for the wrinkling instability of a half space or thick layer loaded in compression, where the elastic modulus varies with depth exponentially or following a distribution of an error function.

Recently, Diab *et al*. [10] completed the second-order instability analysis and compared the results with those of FEA for a flat surface of vai-EDM neo-Hookean solid under lateral compression. In their perturbation analysis, they found that isoperiodic compression with a normalized wavenumber, *ε*_{c} in the range of 0.272<*ε*_{c}≤0.456, where the word ‘isoperiodic’ means that the wavelength is set by boundary conditions. On the other hand, the FEA showed that the compression with *ε*_{c}≤0.272. The results indicate that the second-order instability seems to trigger cascade instabilities of all the higher-order wrinkles to craft the instantaneous creasing for

Because a homogeneous half space of a neo-Hookean solid does not have a characteristic length, the creasing of the half-space has no self-selective periodicity. By contrast, a homogeneous thin film attached on a homogeneous half-space wrinkles with a self-selective critical wavelength, as it has a characteristic length—the thickness of the film. However, the self-selective critical wavelength of the wrinkling on the surface of a vai-EDM neo-Hookean half space diverges to infinity. Although the vai-EDM solid has a characteristic length—the decay length of the modulus—it is similar to a free-standing slab and buckles such as an Euler beam having no self-selective wavelength. Here, in this paper, we analyse ruga-formation instabilities of a half space of neo-Hookean material with an *exponentially decaying modulus approaching a finite* (*bulk*) *value* (afv-EDM). In this representation, the homogeneous and the vai-EDM neo-Hookean half spaces are two special limiting cases, *η*=1 and *η*=0, respectively, of an afv-EDM half space. The bifurcation behaviour of an afv-EDM solid is then compared with that of a homogeneous film on a homogeneous substrate.

The plan of this paper is as follows. In §2, we conduct the second-order stability analysis using the Koiter approach [24,28] for an afv-EDM half space of neo-Hookean material. In §3, FEA is presented for studying ruga formations in afv-EDM neo-Hookean half space under isoperiodic compression. Comparison between the afv-EDM system and a homogeneous thin film on a homogeneous substrate is provided in §4. Finally, a conclusion that summarizes the main findings and observations is presented in §5.

## 2. Bifurcation analysis of a neo-Hookean half space with a graded-stiffness boundary layer

Here, stability analysis of a neo-Hookean half space with a graded-stiffness boundary layer under plane-strain in-plane compression is conducted. Herein, we consider the case of a neo-Hookean system with afv-EDM in general (figure 1). The elastic modulus of the material is considered to be varying as
*η* is the ratio of the shear modulus of the substrate bulk material *Q*_{B} to the shear modulus at the free surface *Q*_{0}, and *a* represents the reciprocal of the decay length of the modulus. Here, the two limiting cases *η*=0 and *η*=1 represent vai-EDM and uniform-modulus half spaces, respectively.

The plane-strain deformation mapping from the stress-free reference configuration *X*_{i} to a deformed configuration *x*_{i} is denoted by
_{(i)} are the in-plane principal stretch ratios of the fundamental state aligned in the (*X*_{1},*X*_{2}) directions, and *u*_{i} the plane-strain perturbation displacement components. For a neo-Hookean material, the strain energy of a plane-strain deformation is defined as
*u*_{i}(*X*_{1},*X*_{2}) and *p*(*X*_{1},*X*_{2}). Here, *F*_{ij} is the two-dimensional deformation gradient tensor defined as *F*_{ij}=∂*x*_{i}/∂*X*_{j}, *J*=det[*F*_{ij}] the Jacobean and *p* as the hydrostatic pressure that plays the role of a Lagrangian multiplier function to enforce incompressibility in the analysis. Summation convention is applied to repeated indices throughout this paper. The principal stretch ratios of the fundamental state are set as λ_{(1)}=λ and λ_{(2)}=1/λ for incompressibility.

For the perturbation analysis, the perturbation displacement *u*_{i}(*X*_{1}, *X*_{2}) and the net hydrostatic pressure *p*(*X*_{1}, *X*_{2}) are expanded about the fundamental state with respect to a perturbation parameter *ζ* up to the second order
*p*^{(0)}(*X*_{1},*X*_{2}) is the hydrostatic pressure in the fundamental state, and *ζ* is taken to be the amplitude of the first-order vertical displacement component at the top surface,

Then, the equilibrium equations up to the second order are reduced as
*et al*. [10], we first solve the linear bifurcation problem to obtain the critical wrinkle bifurcation strain and its corresponding normalized wavenumber, (*η*<1. Throughout this paper, it is important to distinguish the *self-selective* (ss-) critical bifurcation occurring at the minimum critical strain from the bifurcation owing to the imposed *isoperiodic* (ip-) boundary conditions. Calculation details of the linear bifurcation problem are provided in appendix A. Second, stability of the wrinkle bifurcation mode is tested using the energy variation with respect to the displacement variations up to the second order
*φ* denotes the strain energy per period, which is the same as the potential energy per period for the variation, and the superscript *b* indicates that the quantity is evaluated at the bifurcation point. The first- and the second-order displacement fields are orthogonal for the energy variation and the other terms up to the *ζ*^{3} term in equation (2.6) vanish for equilibrium at a bifurcation point. Thus, the nonlinear-bifurcation energy variation is reduced
*i*, *j*, *k*, *m*, *n*=1/2, and a linear operator

Equation (2.8a) exhibits an explicit coupling effect of

Figure 2*a* and *b* shows, respectively, the wrinkle bifurcation strain of the self-selective wrinkle bifurcation state and its corresponding energy variation with respect to the flat state as a function of *η*. As shown in figure 2*b*, for *η*>0.148 energy variation changes from positive to negative, indicating that the self-selective bifurcation state is unstable for afv-EDM systems for *η*>0.148. Using finite-element methods (FEM), Diab *et al*. [10] have shown that once the second-order mode becomes unstable, higher-order modes become unstable simultaneously leading the material to crease instantaneously. Similar cascade instability is expected to occur for the general case when the wrinkle mode becomes unstable.

Next, we investigate instability of afv-EDM systems under isoperiodic compression. The results are presented in figure 3 as a single parameter (*η*) family of bifurcation curves on the *η*<1). The knee shape is a typical bifurcation-curve characteristic of a thin stiff film on a soft substrate. For such a thin-film layer, the film is easier to be buckled for smaller ip-k; however, the energy required to deform the substrate also increases. Competition between these two effects determines the knee shapes of the bifurcation curves in figures 3 and 5. The wavenumber at the front endpoint of the knee curve is the self-selective wavenumber, *M* point as shown in figure 3 for *η*=0.1. There is a cut-off frequency *M* point always appears on the upper bifurcation branch above the knee (ss-) point.

Our analysis shows that there are two more delineating points of wrinkle stability, denoted by *M*′ and *M*′′ as shown in figure 3 for *η*=0.1. In contrast to the *M* point, *M*′ appears near the knee point, sometimes below the knee point (0.006<*η*≤0.1), engendering the second-order instability for smaller *M*′ for 0.006<*η*≤0.189. Then, the *M* point and the *M*′ point merge together at *η*=∼0.189 and disappear together for approximately 0.189<*η*. For *η*<0.006, wrinkle bifurcation is stable for any isoperiodic compression with *η*=0.006, an infinitesimally small critical zone appears at a point *M*′ and *M*′′ points. As *η* increases beyond 0.189, instantaneous creasing commences for all *M*′′. The energy variation with respect to the fundamental state (equation (2.7)) is due to deformation in the boundary layer as well as deformation in the deep uniform substrate. For *η*<0.006, wrinkle instability is mainly determined by strain energy stored in the boundary layer, where for increasing *η* increases, the contribution of the deformation in the deep uniform substrate to the total energy variation of the whole system becomes significant, causing instabilities of the wrinkle state for an ip-wrinkle with *η* as determined by the energy variation in equation (2.7). As *η* increases towards unity, i.e. a homogeneous half space, the critical *M*′′ approaches approximately 0.37, and all critical strains of wrinkling merge to a single value, 0.456, the Biot critical strain.

## 3. Finite-element analysis

In §2, we have performed nonlinear bifurcation analysis to calculate the critical bifurcation strain and its corresponding wavelength for afv-EDM systems to check the stability of the critical wrinkle state. However, our second-order perturbation calculation is valid only near the critical bifurcation point of the fundamental state, and therefore cannot be used to analyse the wrinkle to crease transitions that usually involve bifurcation at a highly non-uniform large deformation state. Here, FEA is used to check the post-buckling deformation of afv-EDM systems with focus on the transition from a wrinkle state to a crease state, along an isoperiodic compression path.

The finite-element simulations are carried out as a two-dimensional plane-strain problem in Abaqus/standard with stabilized Newton–Raphson method. The exponentially decaying modulus of the half-space is modelled using the linearly temperature-dependent elastic modulus option in Abaqus, setting the temperature distribution exponentially decaying without thermal expansion. A unit cell model of width *L* and depth *H*=10 L is used in all of the simulations where the normalized wavenumber *a*, a mesh with uniform size along the span of the unit cell is adopted. Mesh size is increased with distance from the free surface where the six-node quadratic plane-strain triangle elements (CPE6H) are used to make the transition between the regions with different mesh density. Displacement components along the two vertical sides of the unit cell are tied with constraint equations to enforce periodic boundary conditions. The bottom surface is constrained to move only in the horizontal direction. A frictionless self-contact interaction is applied at the top surface of the half-space to prevent any self-penetration of the free surface in case the surface creases or folds to form a channel. Our numerical simulations involve two major steps. The first is to create a deformed mesh matched along a free boundary with a geometrical imperfection of sinusoidal mesh-distortion displacement, *et al*. [10]. Second, the unit cell of the stress-free deformed mesh is then subjected to lateral compression using a displacement-controlled boundary condition with small incremental compression steps, Δ*ε*, less than 0.0005.

Contour plots in figure 4*b*–*f* show nominal compressive strain (−*E*_{11}) distributions in various ruga states created by isoperiodic compression at different *η*=0.001. The surface stably wrinkles with a small ip-*b* *e* *et al*. [10] to diverge at the creasing strain. The critical strain for the setback crease is computed and plotted in figure 5 as a function of *M* points for various *η*. However, the wrinkle–setback–creases are observed only for *η*≤0.189. The flat state can also unstably crease instantaneously for large values of ip-*d* *η* are displayed as dashed curves above *M* points in figure 5. The FEM results are consistent with our theoretical analysis presented in §2. We checked wrinkle to crease transition for the finite range *η*≤0.189 and found that instantaneous crease occurs only for *η*>0.08. In the range of 0.006≤*η*≤0.08, the second-order instability could not guarantee creasing. We also investigated the evolution of a wrinkle towards folding for small *η* as shown in figure 4*c* *η*=0.001. The phase boundary that delineates the wrinkle region from the fold region is computed and plotted as a chained curve below the dashed crease curve in figure 5 for *η*=0.001. Under further compression, crease develops at the free surface inside the folded wrinkle causing extended self-contact fold-crease as shown in figure 4*f* *η*>0.002, our FEM simulations show that wrinkles may never fold for any admissible wavenumber ip-

Our analysis above is limited to the critical bifurcation which is generally revealed by linear perturbation. However, it was shown recently that for the limiting cases, *η*=1 [20,21], there is a subcritical bifurcation of instantaneous creasing. It was numerically demonstrated that the subcritical-bifurcation crease state could be reached by pre-pinching the material surface [22] or through unloading the critically or snap-bifurcated crease state [21]. Here, we compute the subcritical bifurcation points for vai-EDM systems (*η*=0) by conducting a numerical experiment of compressing the system until it critically creases, followed by releasing the compressive strain to have the surface flat again at the lower limit of the subcritical strain. The subcritical bifurcation curves are plotted in figure 6 for vai-EDM systems compressed by different ip-*ε*_{C} and unloading along a part of the subcritical bifurcation branch, followed by jump return to the flat state at the lower-limit strain *ε*_{L} of the subcritical bifurcation branch for *ε*_{L} is still smaller than *ε*_{M}=0.272 at *et al*. [22] for creasing of the homogeneous half space. However, our results additionally show that there exists an unstable subcritical bifurcation barrier hidden between the high-energy flat equilibrium state and the stable low-energy subcritical crease state for *ε*_{HMHS}<*ε*<*ε*_{B}, where *ε*_{B}=0.456. It is also known that the energy levels of the crease and the flat states are equal at *ε*=*ε*_{HMHS} and that the energy level of the flat state is higher than the energy level of the crease state at *ε*=*ε*_{B} for *ε*_{HMH}<*ε*<*ε*_{B} and vanishes at *ε*=*ε*_{HMHS} and *ε*=*ε*_{B}. Finally, we note that this type of crease hysteresis was observed experimentally on the surface of a polymer gel upon reducing the degree of swelling after creasing [29].

## 4. Comparison of an afv-exponentially decaying modulus system with a thin film on a homogeneous substrate

Using the nonlinear bifurcation approach developed in §2, stability analysis of a bilayer neo-Hookean system composed of a film on a deep substrate with modulus ratio, *R*=*Q*_{f}/*Q*_{s}, is conducted, as shown in appendix C. Here, *Q*_{f} and *Q*_{s} are, respectively, the constant moduli of the film and of the deep substrate. All the analytical expressions of the first- and second-order perturbation solution and the necessary coefficients for the analyses are provided in appendix C and the electronic supplementary material. To make a meaningful comparison between the bifurcation behaviour of the two different systems, the effective film modulus of the afv-EDM system and its corresponding thickness is determined by matching the area of the modulus distribution of the boundary layer and the first moment of the same area as (similar to [27], but with the factor *β*)
*β*, *β*=3 gives the closest self-selective wavenumber and the critical strain to those of the corresponding bilayer system. Moreover, our second-order stability analysis shows that, similar to the afv-EDM system, ip-wrinkle bifurcation of the bilayer system is also characterized by the three critical bifurcation points *M*, *M*′ and *M*′′. Appearance of such demarcation points delineating stable wrinkling and unstable creasing is presumably caused by competition between the strain-energy variations in the film and in the substrate as a result of the wrinkle bifurcation. It is believed that the excessive energy to upkeep snap bifurcation of creasing mainly comes from the stored strain energy in the film for

The nominal compressive strain, *ε*, and the normalized wavenumber, *kh*, of the *M* points of the bilayer (solid line) and the afv-EDM (dashed line) systems are shown as a function of the modulus ratio *R* and *a*,*b*. As discussed in §2, there is a critical stiffness ratio, *η*_{m}, above which the *M* point and the *M*′ point coalesce and ip-wrinkle bifurcation becomes unstable for *η*_{m}=0.189), whereas the bilayer system takes the stiffness ratio *R*_{m}=1.71. In figure 3, it is shown that the critical creasing of an afv-EDM half space of *η* greater than 0.148 for which the critical strain is 0.22. On the other hand, ss-wrinkle bifurcation of a bilayer system becomes unstable for *R* less than 1.73 for which the critical strain is 0.33 [14,30].

Figure 7*a* shows that the critical strain of the afv-EDM system at the *M* point is almost invariant, *ε*_{M}≈0.272, for 0≤*η*≤0.1 (*ε*_{M}=0.258 as *η* increases (i.e. *ε*_{M} increases monotonically from *ε*_{M}=0.187 (at *ε*_{M}=0.33 (at *R*_{m}=1.71). Each pair of the dotted and chained curves in figures 7*a*,*b* represents the traces of the *ε*_{M′′} and *k*_{M′′}*h* for the afv-EDM and the bilayer systems as functions of *R*, respectively. As noted earlier, for *R*<1.71, the cut-off frequencies of creasing instability for both afv-EDM and bilayer systems are determined solely by the *M*′′ points. As shown in figure 7*a*, the *M*′′-point critical strain traces of the two systems are very close, and both converge asymptotically to *ε*_{B}=0.456 *as R* approaches 1, i.e. the homogeneous half space. The *M*′′-point critical wavenumbers of the two systems are also close to each other as shown in figure 7*b*; however, they separate to different finite values as *R* decreases to 1, and then both values are supposed to jump down to zero at the singular point, *R*=1. In contrast to the *M*′′-point behaviour, the critical strain and the wavenumber of the *M*-point of an afv-EDM system are quite different from those of a bilayer system. These results indicate that wrinkle bifurcations of stiff-graded boundary layers with long wavelengths can be approximately represented by the behaviour of a bilayer system with its film thickness the same as the effective thickness defined by (4.1) with *β*=3. By contrast, crease bifurcations of a graded boundary layer with short ip-wavelengths *η*→1) cannot be well represented by those of a bilayer system; the stiffness gradient in the boundary layer significantly influences the crease bifurcation behaviour. As a final remark, it is noted that the case of *η*>1 (or *R*<1) is out of the scope of this paper, for which subsurface creasing and its competition with free surface creasing can be of interest.

## 5. Conclusion

In this paper, we present instability analysis of a *graded*-*stiffness boundary layer* near the free surface of an incompressible neo-Hookean material half space under in-plane compression. The boundary layer is represented by the stiffness distribution exponentially decaying from the surface value *Q*_{0} to the bulk value *Q*_{B} with the decay length of 1/*a*. Parametric studies in our first-order stability analysis reveal that such a system buckles with a self-selective wavenumber for the parameter range 0<*η*≡*Q*_{0}/*Q*_{B}<1, similar to the common bilayer system of a homogeneous film on a homogeneous half-space substrate. By contrast, the homogeneous half space (*η*=1) and the vai-EDM solid (*η*=0) do not exhibit such self-selectivity.

In the subsequent second-order perturbation stability analysis supplemented by FEM analysis, we uncover that any periodic sinusoidal perturbation with a wavenumber greater than a critical wavenumber, *η*<1. The lower-limit wavenumber of creasing, *η*<0.189, and jumps down to *η*<1. Such discontinuous behaviour of the lower-limit creasing wavenumber is caused by emergence of a peculiar isolated point on the bifurcation branch at which the wrinkling exhibits cascade instability only up to a finite order for an isolated wavenumber *η*<0.08. On the other hand, the wrinkling develops cascade instability of infinite order, i.e. creasing, within a finite wavenumber range, *η*<0.189. At *η*=0.189, *η*<1, whereas periodic setback creasing arises following progressive stable wrinkle growth for 0≤*η*<0.148.

In our parametric study, creasing of a homogeneous half space can be analysed by taking either *η*→1 for afv-EDM systems, or by taking *ε*_{HMHS}=0.35, is obtained by taking *ε*<0.456.

Comparison between the bifurcation behaviour of the afv-EDM system and its equivalent bilayer system discloses significantly different creasing behaviour between the two systems, whereas the two systems are closely matched for their stable–wrinkling behaviour. Thus, caution has to be exercised in calculating the creasing strain and the creasing wavenumber of neo-Hookean materials with a graded-stiffness boundary layer in terms of the equivalent bilayer system.

## Funding statement

This work was supported by the Korea Institute of Science and Technology, the Korea Institute of Machinery and Materials, and the US National Science Foundation (DMR-0520651) through MRSEC at Brown University.

## Appendix A. The linear bifurcation analysis

Using (2.3) and (2.4), the second-order expansion of the first-Piola–Kirchhoff stress is determined as
*X*_{2}=0 is satisfied at the zeroth order by setting
*ψ*^{(1)}(*X*_{1},*X*_{2}) as
*X*_{1}=0, and thus we write
*η*=0 and *η*=1 [10,15]. However, for 0<*η*<1, the authors could not obtain an analytical solution. Therefore, the well-established compound matrix method is used to determine the critical wavenumber as a function of the nominal strain [31,32,33]. Details on the numerical solution of (A 7) are presented in the electronic supplementary material.

## Appendix B. The second-order perturbation analysis

Given the second-order incompressibility condition in (2.5c), the second-order displacement components are written as

Similar to the first-order solution in appendix A, the second-order perturbation problem is cast as

The components of the matrix **A**^{(2)}(*X*_{2}) are obtained by replacing

The second-order traction-free boundary conditions at *X*_{2}=0 and the second-order zero displacement boundary conditions in the perturbed configurations at **C**^{(2)}=**C**^{(1)} and **B**^{(2)} is obtained from **B**^{(1)} by replacing

## Appendix C. First- and second-order solutions for the bilayer system

The first- and second-order perturbation solutions for the bilayer system composed of incompressible neo-Hookean materials with constant moduli are obtained using the formulations in appendices A and B by taking *η*=1. Analytical expressions for the first- and second-order functions, *Ψ*^{(1)}(*X*_{2}) and *Ψ*^{(2)}(*X*_{2}) are obtained as
*Γ*(*X*_{2}), *E*(*X*_{2}) are deduced from (B 2c) and (B 2d).

The following first-Piola–Kirchhoff stress components are obtained up to the second order as follows

- Received March 17, 2014.
- Accepted May 21, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.