## Abstract

We prove theoretically that when a soft solid is subjected to an extreme deformation, wrinkles can form on its surface at an angle that is oblique to a principal direction of stretch. These oblique wrinkles occur for a strain that is smaller than the one required to obtain wrinkles normal to the direction of greatest compression. We go on to explain why they will probably never be observed in real-world experiments.

This article is part of the themed issue ‘Patterning through instabilities in complex media: theory and applications.’

## 1. Introduction

When a soft solid is subjected to an extreme deformation, its surface (or a part of it) eventually buckles, provided it has not ruptured before. Guided by physical intuition and by experimental evidence (figure 1), we expect the formation of wrinkles or creases to be arranged orthogonally to the direction of greatest contraction, and thus aligned with a principal direction of deformation. However, as we show in this paper, the mathematical equations modelling the development of such surface instabilities do not necessarily predict that they should be such *principal* wrinkles (figure 2*a*). In fact, we find that for the simplest boundary value problem there is, i.e. that of a deformed semi-infinite solid, the theory predicts that *oblique* wrinkles (figure 2*b*) should appear on the free surface prior to the principal wrinkles. We present these results in §2, where we study the theoretical predictions of wrinkle orientation for the most general model of an isotropic, incompressible third-order elasticity soft solid or, equivalently, the Mooney–Rivlin model. These two equivalent models provide a good mathematical description of isotropic homogeneous soft solids such as rubber, silicone or gel, subjected to finite deformations. We establish that, for a certain range of material parameters, the preferred direction of wrinkle formation is not a principal direction of deformation, and that the *obliquity angle* can be as much as 33°, depending on the mode of pre-deformation.

In §3, we report experimental attempts to uncover oblique wrinkles on the surface of a homogeneous block of gelatin subjected to a shear-box deformation [1]. The search for such an occurrence was completely unsuccessful. First of all, we did not manage to generate actual sinusoidal wrinkles on the block, but we always obtained creases, which are indeed expected to occur prior to the wrinkles [2,3]. Moreover, the creases were not oblique, but aligned with the long diagonal of the shear-box, itself a direction of principal stretch.

We moved on to tests with a slightly different set-up, allowing us to finally observe the formation of sinusoidal wrinkles instead of creases. For that set-up, we let the free surface of the gelatin dry overnight to a certain extent and form a stiff film on top of a softer substrate. It is then well known [2] that, for such a layered structure, sinusoidal wrinkles occur first, and not creases (provided the stiffness contrast is large enough); however, we observed that the wrinkles always appeared in a principal direction of deformation.

In §4, we elucidate why oblique wrinkles do not appear experimentally, even in the case of a coated half-space. Using the Stroh formalism [4] to predict their onset, we find that oblique wrinkles should indeed appear, but only as long as the overlaying film is at most two times stiffer than the substrate. If it is stiffer than that, then the first wrinkles to appear are the principal ones according to our calculations. However, Hutchinson & Cao [2] used post-buckling finite-element simulations of neo-Hookean solids to show that creases dominate wrinkles as long as the film is less than about 10 times stiffer than the substrate. Hence, oblique wrinkles, although predicted by the theory, will probably never be captured in an experiment.

## 2. Oblique wrinkles on a half-space

To model isotropic, homogeneous, incompressible *soft solids*, we use a strain energy density *W* of the Mooney–Rivlin type. It is linear in *I*_{1}=tr(**C**) and *I*_{2}=tr(**C**^{−1}), where **C** is the right Cauchy–Green deformation tensor,
2.1where *μ*>0 is the shear modulus, and *β* is a material parameter such that −1≤*β*≤1. Note that these inequalities ensure the strong ellipticity of the incremental equations of equilibrium [5]. At *β*=−1, *W* recovers the so-called neo-Hookean material (linear in *I*_{1}), and at *β*=1, the so-called extreme-Mooney [6] material (linear in *I*_{2}).

As shown by Rivlin & Saunders [7], when we expand the Mooney–Rivlin strain energy density in powers of the Green–Lagrange strain **E**=(**C**−**I**)/2, and neglect terms of order higher than cubic, we find that (2.1) is equivalent to the most general model of isotropic, incompressible, third-order elasticity
2.2where *μ*_{0}>0 is the second-order Lamé coefficient and *A* is the third-order Landau constant. The connection between those constants is [8] *μ*_{0}=*μ*, *A*=−*μ*(*β*+3), which puts bounds on *A* in order to ensure strong ellipticity of the incremental equations: −4*μ*_{0}≤*A*≤−2*μ*_{0}.

To model *finite strains*, we consider homogeneous deformations of a semi-infinite solid. We place ourselves in the Cartesian coordinate system (*x*_{1},*x*_{2},*x*_{3}) of axes aligned with the Eulerian principal directions (along the eigenvectors of **b**, the left Cauchy–Green deformation tensor). We denote the eigenvalues of **b** by , and , and call λ_{1}=λ<1 the principal stretch of contraction, λ_{2} the principal stretch along the normal to the free surface, and λ_{3} the third principal stretch. Because of the incompressibility constraint , we know that λ_{3}=(λ_{1}λ_{2})^{−1}. The homogeneous deformation is maintained by the application of a constant Cauchy stress ** σ**, and we take the stress component along the direction

*x*

_{2}to be zero, i.e.

*σ*

_{22}=0, as the boundary of the half-space is assumed to be free of traction. We will consider the following archetypes of deformation:

(i) uniaxial compression, λ

_{1}=λ, λ_{2}=λ^{−1/2}, λ_{3}=λ^{−1/2};(ii) plane strain, λ

_{1}=λ, λ_{2}=λ^{−1}, λ_{3}=1;(iii) simple shear, λ

_{1}=λ, λ_{2}=1, λ_{3}=λ^{−1};(iv) shear-box, λ

_{1}=λ, , .

We will use λ as our *parameter of bifurcation*, calling it λ_{cr} when the threshold for wrinkles is reached. Note that other quantities can be used, such as the material compressive strain [2] *ε*_{cr}=1−λ_{cr}, or the angle of tilting *ϕ* for shearing deformations (figure 3). For simple shear, *ϕ* is given by [9], and for the shear-box deformation, by [10].

We denote by *θ* the *obliquity angle*, i.e. the angle between the wavefront of the wrinkles and the direction of the largest stretch. At *θ*=0°, the wrinkles are normal to the direction of greatest compression, while at *θ*≠0°, they are *oblique* (figure 2).

For *principal wrinkles*, the wavefront is aligned with *x*_{3} and varies sinusoidally in the *x*_{1} direction, while the amplitude decays exponentially with *x*_{2}. For the Mooney–Rivlin model, the bifurcation criterion for principal wrinkles is independent of the material constants *μ* and *β* (or *μ*_{0} and *A*). It reads as
2.3where *σ*_{0} is the real root of the cubic
2.4This result can be traced back to Biot [11], Flavin [12] and Green & Zerna [13]. In table 1, we list the critical stretches and strains corresponding to the modes of deformation (i)–(iv) for principal wrinkles and for oblique wrinkles on an extreme-Mooney material.

Destrade *et al.* [15] found an explicit secular equation for surface waves in deformed Mooney–Rivlin materials, from which the bifurcation criterion here can be found by taking the wave speed to be zero. We used this explicit secular equation to predict the formation of wrinkles (it is too long to reproduce here). In parallel, we also used the bifurcation criterion based on **z**, the surface impedance matrix [16–19], which is a simple version of a method we will use for layered media in §4. It reads
2.5where **z** is the Hermitian, positive semi-definite solution of the algebraic Riccati equation,
2.6Here, **N**_{1}, **N**_{2}, **N**_{3} are the 3×3 blocks of the 6×6 Stroh matrix (see the appendix), and the superscript ^{†} denotes the Hermitian conjugate. To solve this equation numerically, we rely on Riccati solvers found in Mathematica and Matlab, which give the positive definite solution. The results from this method were the same (within machine precision) as those given by the explicit secular equation found in [15]. There are three short analytical expressions of the bifurcation criterion in some special cases.

For the neo-Hookean material (i.e. *β*=−1 in (2.1), or *A*=−2*μ*_{0} in (2.2)), Flavin [12] found the following bifurcation criterion:
2.7

For the extreme-Mooney material (i.e. *β*=1 in (2.1), or equivalently *A*=−4*μ*_{0} in (2.2)), we find the following explicit bifurcation criterion:
2.8with *σ*_{0} now given by equation (2.7). Using Mathematica, the reals roots of (2.8) can be obtained, and then substituted into equation (2.7) to get an explicit bifurcation criterion.

Finally, for a wrinkle-front aligned with the long diagonal of the shear-box (*θ*=0°), we have again *principal* wrinkles, and equation (2.3) applies.

For cases (i)–(iv), we find that there is a range of values for *β* starting at *β*=−1 (neo-Hookean solid), where principal wrinkles appear first, and another range up to *β*=1 (extreme-Mooney solid), where we expect oblique wrinkles to appear. Figure 4 illustrates this occurrence for the shear-box deformation (case (iv)), in which we see that principal wrinkles are expected when , and oblique wrinkles when . For the extreme-Mooney material, the obliquity angle is *θ*≃33°, and the critical stretch is λ_{cr}=0.539, compared with λ_{cr}=0.471 for principal wrinkles. We focused on this deformation because of its ease of experimental implementation (see the next section), but the other deformations yielded a similar behaviour, with more or less marked differences between the values of λ_{cr} for principal wrinkles and for oblique wrinkles, as summarized in table 1.

In general, the numerical procedure to find the critical values of the compressive stretch and of the corresponding obliquity angle is straightforward and robust. After non-dimensionalization, the critical values turn out to be independent of *μ* (or *μ*_{0}), and so we can plot λ_{cr} against *θ* for a given solid characterized by the single material parameter *β* (or the ratio *A*/*μ*_{0}).

## 3. Experimental buckling with the shear-box deformation

To build a shear-box, we assembled four acrylic plates (10×10×1.0 cm) into a cube with four hinged edges. We then prepared 12 different gels, using commercial gelatin in several forms: solid leafs, powders or concentrated cubes to be dissolved in water. We varied the concentrations, from that prescribed by the manufacturers to three times the normal concentration. We poured about 300 ml of gelatin into the shear-box and let it set into a homogeneous solid phase, either at room temperature or at refrigerated temperature, from 4 to 12 h depending on the specimen.

The gelatin block was then subjected to a shear-box deformation (case (iv)), simply by applying manual pressure on two opposite hinged edges of the box (see figure 3*b* for a schematic representation of the deformed state). To ensure that the solid deformed homogeneously, we injected olive oil between the sides of the block and the walls to enhance the gliding. This precaution led to a central area where it was reasonable to consider that the deformation was homogeneous, i.e. where initially parallel lines remained parallel in the current configuration (figure 5*a*). We carried out the deformation manually until the free surface buckled as noticed by the naked eye; see examples in figure 5. The corresponding critical tilting angle *ϕ* was recorded.

Experimentally, we noted that we did not observe the formation of small-amplitude sinusoidal wrinkles but of creases instead (figure 5*a*). Moreover, the creases were always aligned with the long diagonal of the shear-box, and thus not oblique. We measured that the creases appeared at *θ*=38°±2°, irrespective of the concentration and type of gelatin used for the gel sample.

Hence, we conclude that the theoretical predictions of §2 were not verified experimentally. The main problem is that of crease formation trumping wrinkles; this was to be expected, as creases have been predicted to appear much earlier than wrinkles on homogeneous blocks. For instance, Hong *et al.* [20] showed through finite-element simulations that creases appear on a neo-Hookean half-space deformed in plane strain at a 35% compressive strain, so that the 46% required for sinusoidal wrinkles is never attained.

To bypass this problem, we focused instead on wrinkles appearing when a soft semi-infinite solid is coated with a thin stiff layer. Those are then remarkably stable and their onset is well predicted by incremental analysis; they develop into nonlinear patterns only with further compression; see Cao & Hutchinson [2] or Jin *et al.* [21] for details. Figure 6 displays experimental evidence of the formation of regular sinusoidal wrinkles for bending and torsion deformations.

Experimentally, we let the gelatin blocks dry out overnight at refrigerated temperatures, so that a skin formed on their surface due to surface dehydration during the longer cooling process. When deformed in the shear-box, the buckling surface now exhibited regular, sinusoidal wrinkles, at least in the centre of the block, where the homogeneous deformation was taking place (figure 5). We noted that the wavelength of the wrinkles was much smaller than the dimensions of the block, and thus that the half-space idealization was justified for the modelling (for an analysis of the influence of the finite depth of the block, see the recent work by Jin *et al.* [21]). However, oblique wrinkles were not observed, but principal wrinkles instead, again at a tilting angle of 38°±2° (corresponding to a strain of *ε*_{cr}=38%). In the next section, we explain why this absence of experimental oblique wrinkles is indeed in agreement with theoretical predictions.

## 4. Wrinkles on a coated half-space

### (a) Layered structure

We consider that the top surface of the gelatin block (in contact with air) has dried out and stiffened faster than the rest of the enclosed block during the solidification process. We model the resulting heterogeneous structure as a composite material: a thin film of thickness 2*h*, say, in the region −*h*≤*x*_{2}≤*h*, in bonded contact with a substrate in the region . The strain energy of each solid is expanded up to third order of Green strain as
4.1respectively. Now consider the mechanical response of these two solids to common modes of deformation such as uniaxial tension under Cauchy tensile stress *T* or simple shear under Cauchy shear stress *S*. Then, up to second order in the strain, we have for the film and for the substrate
4.2respectively, and similarly for other modes of deformation. (Here, the measure of strain *ε* is in turn the elongation 1−λ (from case (i)) and the amount of shear λ^{−1}−λ (from case (iii)).) A consistent way to ensure that the response of the film is always ‘stiffer’ than that of the substrate is to take
4.3where *Γ*>1 is defined as the stiffness ratio. For the constants of the Mooney–Rivlin material (2.1), this translates as
4.4Then the film is always stiffer than the substrate, in the sense that the magnitude of its response to any mechanical stress is *Γ* times that of the substrate.

### (b) Principal wrinkles

For our layered structure, the formulae we use below can readily be deduced by specializing the results of [22,23] to the present context. We write the bifurcation criterion for principal wrinkles as
4.5where **z**_{f} and **z**_{s} are non-dimensionalized surface impedance matrices for the film and the substrate, respectively. The substrate impedance is given explicitly by **z**_{s}=−i**BA**^{−1}, where
4.6and, as before, λ_{1} is the principal stretch of contraction and λ_{2} is the principal stretch along the normal to the surface. The film impedance is given by , where **M**_{1}, **M**_{3} are the respective top-left and bottom-left 2×2 submatrices of the matricant . Here,
4.7the overbar denotes the complex conjugate and *k* is the wavenumber of the wrinkle. It is easy to check that is equivalent to the bifurcation criterion (2.3), as expected. Also, note that here the bifurcation criterion is independent of the material parameters *β*_{s} and *β*_{f}; only the stiffness ratio *Γ*=*μ*_{f}/*μ*_{s} matters.

We may then follow the same strategy as Cao & Hutchinson [2] to find the critical amount of deformation signalling the onset of sinusoidal wrinkles. When *Γ*=1, the system is homogeneous and non-dispersive: the wrinkles appear when ; see equation (2.3). When *Γ*>1, the film/substrate structure is dispersive and to each value of *kh* corresponds a value of for which equation (4.5) is satisfied. By varying *kh*, we can find the maximal value for , which will be the critical value at which the structure wrinkles. Then by varying *Γ*, we can plot figure 7*a*, for all the considered cases of deformation. The routine described in the flowchart of figure 7*b* allows us to obtain a general curve of the maximum critical stretch ratio against *Γ*. Then, the critical value of is related to the critical value of stretch or strain according to each different mode of deformation, cases (i)–(iv).

Figure 7*a* recovers fig. 2 of Cao & Hutchinson [2] when the film and the substrate are both neo-Hookean and the deformation is plane strain. In the figure, we also report the analytical solution of Allen [24] for the buckling of a thin strut on a half-space in linear elasticity. As pointed out in [2], it becomes reliable when *Γ* is large, i.e. when the deformation allowed before wrinkling is small and close to the linear regime.

### (c) Oblique wrinkles

For oblique wrinkles, we must first integrate numerically the following differential Riccati equation for the film’s surface impedance matrix **z**_{f}=**z**_{f}(*x*_{2}),
4.8from its initial value [22] **z**_{f}(−*h*)=**0** to its value **z**_{f}(*h*). Then, we compute the substrate’s (constant) surface impedance matrix **z**_{s} as the positive semi-definite solution to the algebraic Riccati equation (2.6). Finally, we adjust the compressive stretch λ until the following bifurcation criterion is met:
4.9Note that here **z**_{f} and **z**_{s} are dimensional quantities, in contrast with the impedance matrices from §2 because **N**_{1}, **N**_{2}, **N**_{3} in (4.8) (in (2.6), respectively), depend on *μ*_{f}, *β*_{f} (*μ*_{s}, *β*_{s}, respectively).

Here, the analysis is further complicated by the increase in the number of parameters. For a given stiffness ratio *Γ* and a given deformation, we must vary not only *kh* to find the critical stretch, but also the obliquity angle *θ* and the material parameter *β*_{f}=*β*_{s}. Figure 8*b* displays the flowchart for the overall procedure used. (In particular, we used the Matlab code *fminsearch* to maximize the function .) For our purposes it suffices to record, for a given deformation, the smallest critical strain *ε*_{cr} obtained for an oblique wrinkle for a given *β*_{s} and to plot that curve. Then by varying *β*_{s} from −1 to 1, we obtain a family of curves. Here, we present the results for the shear-box deformation, to compare them to our earlier experiments (we found similar figures for the other deformations).

Figure 8*a* shows that oblique wrinkles (coloured curves underneath the thick black curve for principal wrinkles) precede principal wrinkles only when for the materials close to the extreme-Mooney model (*β*_{s}≃1.0) and for most materials (*β*_{s} away from 1.0). (Specifically, we found that the extreme-Mooney material predicts principal wrinkles when *Γ*≥2.20,1.50,3.62,2.28 for cases (i)–(iv), respectively.)

In conclusion, when the stiffness ratio *Γ* is small, creases precede wrinkles; when it is large enough to support wrinkles, principal wrinkles precede oblique wrinkles. So, even though oblique wrinkles exist as first bifurcation mode in theory, they will never be observed experimentally.

## 5. Discussion

To complete the paper, we must qualify our concluding statements above.

First of all our stiffness ratio *Γ* is one measure among others of the stiffness contrast between film and substrate. In practice, it is unlikely that the proportion between *μ*_{0f} and *μ*_{0l}, and *A*_{f} and *A*_{s}, should be exactly the same as in (4.3). So a whole host of possibilities has been missed in this investigation.

Another caveat is that the finite-element simulations of creases evoked in the paper all used the neo-Hookean strain energy density, which does not predict oblique wrinkles. Similarly, semi-analytical investigations into localized solutions, such as the one by Fu & Ciarletta [25], also rely on neo-Hookean models. It thus remains an open question whether oblique wrinkles can precede creases in simulations of compressed Mooney–Rivlin solids. However, our study indicates that they could theoretically only appear in a narrow range of low stiffness contrast, where experiments (our table-top ones, and the comprehensive set of controlled ones presented by Jin *et al.* [21]) have so far failed to exhibit any evidence of their existence.

Also, guided by our experiments, we adopted the half-space approximation for the substrate, because the observed wavelength was small compared with the dimensions of the block. If the thickness *h* of the stiff layer were comparable to that of the soft substrate, it might be possible that the finite depth of the substrate plays a role in the search for oblique wrinkles. The competition between creases and wrinkles is also affected in that regime, and the onset of one compared with the other is more complicated than here; see the analysis of Jin *et al.* [21] for neo-Hookean solids.

## Authors' contributions

A.L.G. carried out the experimental work and the calculations for surface instability, participated in the design of the study and drafted the manuscript. M.C. extended the numerical calculations to the case of oblique wrinkles on a coated half-space, helped in drafting the manuscript, and provided flowcharts, schemes and figures coming out of numerical calculations. M.D. conceived of the study, designed the study, coordinated the study and helped draft the manuscript. A.G. contributed with discussions about the study and helped in drafting the manuscript. All authors gave final approval for publication.

## Competing interests

The authors declare that there are no competing interests.

## Funding

The financial support of the Hardiman Foundation (NUI Galway) and of the Government of Ireland Postgraduate Scholarship Scheme (Irish Research Council) for A.L.G. and of Politecnico di Torino for M.C. is gratefully acknowledged.

## Acknowledgements

Part of this work is based on Chapter 1 of A.L.G.’s PhD thesis at NUI Galway, Ireland (2011–2015). We also thank Molly Aitken, Jorge Bruno, Juliet Destrade, Lucille Destrade, Sean Leen, Joanne McCarthy, Deirdre McGowan Smyth and Giuseppe Zurlo for their technical input.

## Appendix A

The submatrices **N**_{1}, **N**_{2} and **N**_{3} of the Stroh matrix can be found from the general expression in [15] as
where, for a Mooney–Rivlin material (2.1) or equivalently a third-order elasticity solid (2.2),

## Footnotes

One contribution of 13 to a theme issue ‘Patterning through instabilities in complex media: theory and applications.’

- Accepted August 9, 2016.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.