Royal Society Publishing

Ruga mechanics of creasing: from instantaneous to setback creases

Mazen Diab, Teng Zhang, Ruike Zhao, Huajian Gao, Kyung-Suk Kim


We present mechanics of surface creasing caused by lateral compression of a nonlinear neo-Hookean solid surface, with its elastic stiffness decaying exponentially with depth. Nonlinear bifurcation stability analysis reveals that neo-Hookean solid surfaces can develop instantaneous surface creasing under compressive strains greater than 0.272 but less than 0.456. It is found that instantaneous creasing is set off when the compressive strain is large enough, and the longest-admissible perturbation wavelength relative to the decay length of the elastic modulus is shorter than a critical value. A compressive strain smaller than 0.272 can only trigger bifurcation of a stable wrinkle that can prompt a setback crease upon further compression. The minimum compressive strain required to develop setback creasing is found to be 0.174. If the relative longest-admissible perturbation wavelength is long enough, then the wrinkle can fold before it creases, and the specimen can be compressed further beyond the Biot critical strain limit of 0.456. Various bifurcation branches on a plane of normalized longest-admissible wavelength versus compressive strain delineate different phases of corrugated surface configurations to form a ruga phase diagram. The phase diagram will be useful for understating surface crease, as well as for controlling ruga structures for various applications, such as designing stretchable electronics.

1. Introduction

Free-surface deformation instability of a soft solid, with its modulus varying in depth, often develops surface patterns, such as wrinkles, creases, ridges or folds, under in-plane compression. Here, a large-amplitude state of a wrinkle, crease, ridge or fold is collectively denoted by the Latin word ruga, and the study of motions associated with configurational changes among such states is called ruga mechanics. The motions include progressive amplitude growth of a ruga and the snap jump of the amplitude from one state to another. For example, it is commonly observed that a free surface of a homogeneous soft solid develops a crease with a snap-buckling instability, and a soft substrate with a stiff layer progressively grows stable wrinkles [1,2]. So far, much of ruga mechanics research has been limited to post-bifurcation analysis of relatively small-amplitude wrinkles as incipient of ruga structures. In this paper, we present post-bifurcation stability analysis of various ruga states up to very large amplitudes for a neo-Hookean solid [3,4], with its modulus decaying exponentially from the free surface [5,6].

Regarding small-amplitude stability analysis of a flat surface, Biot [7] showed, in his linear bifurcation analysis, that the flat surface of a homogeneous neo-Hookean solid becomes unstable at the Biot critical strain εB of 0.456 plane-strain compression. However, this critical strain can be hardly reached experimentally because the flat surface unstably creases at a strain lower than the Biot critical strain. Gent & Cho [1] observed in their experiments that a rubber-strip surface compressed by bending creased at a nominal strain in the range 0.28–0.42. Later, Hohlfeld [8] used finite-element analysis (FEA) to show that such a crease is typically formed through a snap-buckling process at the Biot critical strain, and that post-bifurcation unloading shrinks the crease depth until it disappears at a compressive strain of 0.35. Subsequently, Hong et al. [9] could trigger progressive crease growth by pre-pinching the surface in their FEA so as to follow the stable subcritical creasing branch of pitchfork bifurcation beyond the compressive strain of 0.35. Cao & Hutchinson [10] studied the imperfection sensitivity of creasing by analysing multiple-mode post-bifurcation stability, with a first-order displacement perturbation in the Koiter approach [11,12]. They found that the post-bifurcation wrinkle of an arbitrary wavelength is unstable, and the critical strain is sensitive to the imperfection amplitude. They also noted in their FEA that the local strain to set off creasing is always the Biot critical strain.

Now, questions arise. What distribution of the substrate modulus can stabilize the post-bifurcation wrinkle? If it can be stabilized, can we sustain the stability up to an arbitrary level of compression? If the stable post-bifurcation wrinkle becomes unstable at a large compression, does the second post-bifurcation ruga become stable or unstable? If unstable, can the unstable bifurcation develop period-multiplication modes [13]? To answer these questions systematically, we use a neo-Hookean solid model with an exponentially decaying modulus (EDM) for the study of ruga development. Then, analysis of the EDM model provides the deformation characteristics of a half-space from those of a homogeneous medium to those with a stiff film [1416] as a single-parameter family. The model was also used for different purposes by Biot [5] and Lee et al. [17]. Behaviour of the EDM model is analytically tractable, and the results are transformative for making interpretations of homogeneous material behaviour by taking the limit of the decay length to infinity. It is shown in this paper that the EDM model can capture not only the behaviour of instantaneous creasing, but also that of setback creasing delayed by progressive growth of wrinkles or folds. It is revealed in later sections that Koiter's approach with second-order perturbation analysis can determine whether the nonlinear bifurcation mode is a progressive wrinkle or an instantaneous crease. In turn, the Koiter analysis provides a criterion to delineate instantaneous creasing from setback creasing. However, once the wrinkle amplitude becomes large, nonlinear deformation processes in the post-bifurcation growth of wrinkles or folds are hardly tractable analytically. Therefore, FEA is used to analyse setback creasing.

Figure 1 shows the schematics of possible creasing processes on an EDM solid surface under isoperiodic compression. The creasing is found to be a snap-buckling process. The creasing process is, in general, sensitive to loading boundary conditions and/or imperfections. We thereby limit our analysis of bifurcation processes to those caused by uniaxial isoperiodic compression. In isoperiodic compression of EDM solids, the longest-admissible wavelength is typically set by the domain size in a numerical analysis or the specimen size in an experiment. In addition, the periodic boundary condition suppresses period multiplying bifurcations, because the wrinkle with the longest wavelength has the lowest energy at a given compressive strain. As a consequence of analysing isoperiodic compression of an EDM solid, we can construct a primary ruga phase diagram on a plane of the smallest admissible nominal wavenumber and the nominal compressive strain. In the phase diagram and as depicted in figure 1, the phase boundaries are constructed by analysing transitions from one phase to another, following isoperiodic transitions from a flat to a wrinkle (C12), a flat to a crease (C13), a wrinkle to a crease (C2c), a wrinkle to a fold (C24) and a fold to a crease (C4c).

Figure 1.

Schematic of ruga phases and their bifurcation processes in a neo-Hookean solid with its elastic modulus decaying exponentially with depth, under lateral compression: (a) flat phase (I) at a stress-free reference configuration in its material coordinate (X1, X2), where Q0 denotes the shear modulus at the surface; (b) progressive-wrinkle phase (II) formed by transition C12 under compressive strain εw, showing transition C2c to a wrinkle-setback crease phase (IIc) in the inset; (c) instantaneous-crease phase (III) made by transition C13 under compressive strain εc; (d) progressive-fold phase (IV) formed by transition C24 under compressive strain εf, showing transition C4c to a fold-setback crease phase (IVc) in the inset.

2. Perturbation analysis of creasing processes: instantaneous versus setback creasing

Using a linear perturbation analysis, Biot [6] examined bifurcation of a flat neo-Hookean EDM solid surface to a sinusoidal wrinkle under plane-strain uniaxial compression. He determined the bifurcation wavenumber as a function of the applied compressive strain. However, the linear perturbation analysis cannot determine the stability of post-bifurcation modes and their amplitudes. Here, we revisit the Biot problem in §2a and extend the study to nonlinear bifurcation analysis of the wrinkle by carrying out Koiter analysis with second-order perturbation solutions in §2b.

(a) Linear perturbation analysis of surface wrinkles on a neo-Hookean exponentially decaying modulus solid

For the perturbation analysis, the material coordinate (X1,X2) and the distribution of the shear modulus Embedded Image2.1ain the substrate are depicted in figure 1a, where Q0 is the shear modulus at the surface, and the coefficient a represents the reciprocal of the decay length of the modulus. In the deformation analysis, the plane-strain mapping from the stress-free reference configuration Xi to a deformed configuration xi is denoted by Embedded Image2.1bHere, λ(i) are the in-plane principal stretch ratios of the fundamental state aligned in (X1,X2) directions, and ui the plane-strain perturbation displacement components. For a neo-Hookean material, the strain energy of a plane-strain deformation is defined as Embedded Image2.1cfor variations of ui(X1,X2) and p(X1,X2). Here, Fij is the two-dimensional deformation gradient tensor defined as Fij=∂xi/∂Xj, J=det[Fij] the Jacobean, and p the hydrostatic pressure that plays the role of a Lagrangian multiplier function to enforce incompressibility in variational 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 displacement ui(X1,X2) and the hydrostatic pressure p(X1,X2) are expanded with respect to a perturbation parameter ζ up to the second order, Embedded Image2.2aand Embedded Image2.2bwhere Embedded Image and p(0)(X1,X2) are, respectively, the uniform displacement from the fundamental-state configuration and the hydrostatic pressure in the fundamental state. The parameter ζ is taken to be the amplitude of the first-order vertical displacement component at the top surface, Embedded Image. The superscripts in parentheses indicate the perturbation order. The first Piola–Kirchhoff stress, defined for an incompressible solid as Embedded Image2.3is used with the independent pressure expression in (2.2b) to obtain the explicit form of the stress field up to the second order as given in (A 12).

Then, the equilibrium equations up to the second order are reduced from (A 12) to Embedded Image2.4aand Embedded Image2.4band the incompressibility condition up to the second order to Embedded Image2.4cThe first-order displacement components and hydrostatic pressure in the perturbed configuration are determined by solving (2.4), as presented in appendix A, Embedded Image2.5a Embedded Image2.5b and Embedded Image2.5cwhere Embedded Image2.6aand Embedded Image2.6bwhere Embedded Image is the normalized wavenumber in the undeformed configuration. The exponents α3 and α4 have negative real parts. Thus, the constants C3 and C4 in (2.5) are set to zero to ensure that the perturbed displacement components vanish at Embedded Image.

Enforcing the traction-free boundary conditions at the top surface of the substrate up to the first order, we construct the following highly nonlinear eigenvalue problem: Embedded Image2.7Then, the characteristic equation provides the normalized bifurcation wavenumber as a function of the stretch ratio for λ<1 as Embedded Image2.8

Figure 2a shows the normalized bifurcation wavenumber Embedded Image in the undeformed configuration in terms of the applied nominal compressive strain, ε=1−λ. The solid curve in figure 2a indicates the formation of stable wrinkles for compressive strains smaller than 0.272, the transition strain at M. The dashed curve represents the formation of unstable wrinkles under compressive strains between εM (=0.272) and the Biot critical strain, εB (=0.456). The wavenumber Embedded Image diverges to infinity at the Biot critical strain, indicating that the homogeneous neo-Hookean half-space unstably bifurcates at the Biot critical strain for perturbations with an arbitrary wavelength. The chain curve shows a lower bound of the crease strain, which could be gradually reached when the bifurcated crease was unloaded in our FEA.

Figure 2.

(a) Traces of linear bifurcation points on the Embedded Image plane, exhibiting the critical M-point that delineates the instantaneous creasing (dashed curve) and progressive wrinkling (solid curve) traces. The chain curve represents the lower bound of crease strain for neo-Hookean EDM solids with its asymptote at εH. (b) Energy variation δφb normalized by ζ4 up to the second-order perturbation mode showing that δφb is positive for 0<ε<εM=0.272 for which a progressive wrinkle develops, and that δφb is negative for εM<ε<εB=0.456 for which an instantaneous crease forms.

(b) Second-order perturbation analysis for nonlinear stability of a wrinkle

In this section, nonlinear stability of the wrinkles along the solid and the dashed curves presented in figure 2a is examined by second-order perturbation analysis. A solution procedure of the second-order perturbation is described in appendix A, and only the final expressions of the second-order displacement components and the hydrostatic pressure are presented here for further discussions on energy variations, Embedded Image2.9a Embedded Image2.9b and Embedded Image2.9cwhere explicit expressions of Fm (m=1,2,3) are given in appendix A(c), and Bm,H1(X2) and H2(X2) are given in appendix B.

Nonlinear stability of the wrinkle near the bifurcation point is then tested by the energy variation with respect to the displacement variations up to the second order. The energy variation, with λ and Embedded Image fixed at a bifurcation point on the bifurcation curves in figure 2a, can be obtained from (2.1), (2.2), (2.5) and (2.9) with respect to the displacement variation Embedded Image as Embedded Image2.10Here, φ 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. In (2.10) as well as in the following (2.11) and (2.12), the chain rule operation Embedded Image implies variation δφb of the functional φb, caused by variation of the displacement field, Embedded Image. Because the first- and the second-order displacement fields are orthogonal for the energy variation and the other terms up to the ζ3 term in (2.10) vanish for equilibrium at a bifurcation point, the expression of the nonlinear bifurcation energy variation is reduced to Embedded Image2.11The terms on the right-hand side of equation (2.11) are determined for a neo-Hookean EDM solid as Embedded Image2.12aand Embedded Image2.12bwhere i, j, k, l, m, n=1 or 2, and a linear operator Embedded Image means Embedded Image (U(m)W(n))=U(1)W(2)+U(2)W(1).

Equation (2.12a) exhibits the explicit coupling effect of Embedded Image and Embedded Image for the energy variation, whereas (2.12b) entails the effects of implicit coupling between the first- and second-order fields.

Figure 2b shows the nonlinear bifurcation energy variation along the solid and the dashed bifurcation curves in figure 2a. As shown in figure 2b, the energy variation is positive at every bifurcation strain below the nominal strain εM=0.272. This indicates that the bifurcated wrinkle is stable near the bifurcation point within 0<ε<εM, for which Embedded Image. However, the energy variation becomes negative for the bifurcation points within εM<ε<εB and Embedded Image. The result implies that the bifurcated wrinkle is unstable near its bifurcation point in this range, and its growth includes simultaneous unstable growth of the second-order mode of the wrinkle. Once the wrinkle becomes simultaneously unstable with the second-order mode, the wrinkle becomes unstable, regardless of all higher-order modes, because all the modes higher than the second order will contribute to the energy variation only up to O(ζ5). As we can observe in (2.5) and (2.9), the second-order mode escalates the local curvature and compressive strain at the coordinate origin, localizing the valley of the wrinkle. It is observed in the following sections that the nonlinear bifurcation instability up to the second-order mode implies creasing of the surface. Consequently, the M-point (εM=0.272; Embedded Image) breaks up Biot's bifurcation curve for 0<ε<εB in figure 2a into a progressive-wrinkle curve for 0<ε<εM and an instantaneous-crease curve for εM<ε<εB.

3. Finite-element analysis of post-bifurcation wrinkle growth and setback creasing

In §2b, second-order perturbation analysis predicted that isoperiodic compression of a flat EDM solid surface sets off progressive wrinkling or instantaneous creasing. In this section, we carry out finite-element method (FEM) simulations to trace the bifurcation modes and corresponding post-bifurcation processes. Isoperiodic compression is simulated for various wavenumbers in the range Embedded Image.

The FEM simulations are conducted as a two-dimensional plane-strain problem in ABAQUS/standard using the stabilized Newton–Raphson method. Hybrid elements CPE4RH are adopted to simulate the incompressible neo-Hookean material. For the isoperiodic compression process, periodic boundary conditions are enforced on the two sides of a planar unit cell in the (X1, X2)-plane, with width L and depth H=5L, as shown in figure 3a. To introduce surface imperfections, the simulations are first run by enforcing a sinusoidal displacement boundary condition at the top surface as Embedded Image3.1where L is chosen as the imperfection wavelength in the undeformed configuration, and Embedded Image as the relative imperfection amplitude. Then, the deformed geometry is imported, and the imperfection command in ABAQUS is used to create a perturbed mesh that conforms to the surface undulation of the half-space. A frictionless self-contact interaction is applied at the top surface of the half-space to capture the crease and fold behaviour without self-penetration.

Figure 3.

Finite-element analysis of ruga states: (a) finite-element meshes, full domain of L×H in the left inset, and a subdomain of (L/4)×(∼H/40) on the right; (bf) plots of nominal strain Embedded Image denoted by NE11 in the inset of (b), for various ruga states: instantaneous crease (III), progressive wrinkle (II), wrinkle setback crease (IIc), fold (IV) and fold setback crease (IVc).

The classical Lagrange multiplier method of constraint enforcement was adopted using the hard contact option in ABAQUS. The minimal mesh size is about 1/2000 of L. The shear modulus distribution (Q=Q0eaX2) is built up in the ABAQUS program by assigning a fictitious temperature distribution θ=θ0eaX2 with a linear relationship between Q and θ, and a vanishing thermal expansion coefficient.

As shown in the left inset of figure 3a, the isoperiodic compression process is simulated by compressing the domain with periodic boundary conditions and small incremental compression steps, Δε, less than 0.0005 or 0.0008. In the simulation, the normalized wavenumber Embedded Image is controlled by varying a for a fixed domain length L (=2π/k). Figure 3 shows five different ruga configurations, observed near various bifurcation points during isoperiodic compression of an EDM solid with three different Embedded Image Embedded Image, 0.471 and 0.100). The colour plot exhibits the level of the nominal (or Biot) strain Embedded Image. In each frame, a single period result is repeatedly displayed twice to show a configuration of two periods. Figure 3b shows the post-bifurcation configuration of an instantaneous crease near the bifurcation point for Embedded Image, whereas figure 3c presents the wrinkle near the Biot bifurcation point and figure 3d the setback crease at the wrinkle valley for Embedded Image. This is the normalized wavenumber with which a setback crease can be generated by the minimum applied strain, εT=0.174. Below this strain, a crease cannot be formed, regardless of whether it is an instantaneous or a setback crease. This critical strain contrasts εM=0.272, below which an instantaneous crease cannot be created. Figure 3e and f respectively display a folding configuration and setback creasing at the bottom of the fold channel for Embedded Image.

Because we use an imperfection on the surface profile as shown in (3.1), we inspect the evolution of the surface profile by taking the Fourier-series expansion of the deformed surface profile as a function of X1 at every step of the compression, as Embedded Image3.2for N up to 8, where y0 is the profile-base constant. Figure 4ac shows evolutions of the Fourier coefficients (A1, A2, A4, A8) as the domain is compressed with three different Embedded Image (Embedded Image, 1.30 and 2.51) held fixed, respectively. Figure 4d exhibits the growth rates of the lowest-mode coefficient, dA1/dεs, as a function of ε for the three different compression processes. Neglecting the coefficients less than 10−5, figure 4a shows that all the coefficients grow quickly near Biot's wrinkle bifurcation point (ε(i)w=0.157, Embedded Image); however, the growth rates are finite, and all are stable. The growth of the coefficients indicates sharpening of the wrinkle valley. The growth rate distinctly defines the bifurcation strain ε(i)w and thus the Biot wrinkle bifurcation point, as shown in figure 4d. The bifurcation strain ε(i)w is relatively insensitive to imperfections, and coincides closely with the prediction of the linear perturbation analysis. However, growth rates of all the coefficients diverge at the crease point, regardless of an instantaneous or setback crease, and the instantaneous crease point is imperfection sensitive. On the other hand, a setback crease requires substantial amplitude growth of the wrinkle before it locally sets off unstable creasing, and thus the setback crease is relatively insensitive to initial long-wavelength imperfections. However, the setback crease is sensitive to imperfections with wavelengths shorter than the final local zone size of the crease in the wrinkle valley.

Figure 4.

Evolution of Fourier-series coefficients (A1, A2, A4, A8) of surface profiles up to creasing, solid line for A1, dashed line for A2, dotted line for A4 and chained line for A8: (a) wrinkle-setback creasing under isoperiodic compression with Embedded Image; (b) instantaneous creasing at the critical M-point, compressed with Embedded Image; (c) instantaneous creasing with Embedded Image; (d) growth rate, dA1/dε, for creasing of (ac) showing distinct wrinkle bifurcation point at ε(i)w for progressive-wrinkle growth (i), and divergence of the rate at creasing points. The arrows in (ac) indicate the Fourier coefficients for the completed states of creasing.

4. Discussion

In figure 5, three bifurcation curves, C12, C13 and C2c define four basic domains of surface configurations on the Embedded Image-plane. The curves C12, C13 and C2c correspond to the progressive-wrinkle, instantaneous-crease and setback-crease bifurcation curves, respectively. The four domains on the Embedded Image-plane constitute the basic structure of the primary ruga phase diagram, on which Embedded Image implies the homogeneous half-space limit and Embedded Image the thin-film limit. The phase diagram displays possible equilibrium configuration phases of an EDM solid surface under isoperiodic compression. Phase (I) denotes flat, phase (II) wrinkle, phase (IIc) setback-crease and phase (III) instantaneous-crease configurations. Regarding the setback-crease curve, as the compressive strain isoperiodically increases beyond Biot's bifurcation point, nonlinear instability is pre-empted, and a localized setback crease develops at the valley of a wrinkle, preventing the wrinkle from growing further. Collection of the new setback-crease bifurcation points produces the setback-crease curve C2c on the phase diagram, the dashed line in figure 5, which delineates phase (IIc) from phase (II).

Figure 5.

Primary ruga phase diagram for isoperiodic compression of the EDM neo-Hookean solid: solid curves C13 and C12 split by the critical M-point denote, respectively, linear bifurcation (LB) curves of instantaneous creasing and progressive wrinkling, and dashed curves C2c and C4c represent setback creasing (C). The dotted line C24 in the inset indicates folding (F). T is the strain-limit point below which a crease cannot form. R is the triple point of phases wrinkle (II), fold (IV) and crease (IIc and IVc). The dashed line superposed on the solid line C13 shows agreement between the FEM and analytical results for instantaneous creasing.

Furthermore, we found another phase of fold configurations, phase (IV). When we tracked the isoperiodic evolution of wrinkles with low wavenumbers, i.e. less than 0.471, we found that there is a critical wavenumber Embedded Image below which the wrinkles fold at C24 and form channels, as shown in figure 3e. Once the fold configuration is further compressed, a setback crease develops at the root of the fold channel (figure 3f). This creasing at C4c separates the phase of a fold-crease configuration (IVc), from that of a fold configuration (IV), as shown in the inset of figure 5. Beyond this point, a question still remains whether this crease would further propagate to collapse the channel, defining another phase on the phase diagram. It is worth noting that the phases (IV) and (IVc) sit past Biot's singular strain limit. An EDM solid can be isoperiodically compressed beyond the Biot singular strain limit with a low wavenumber Embedded Image because growth of a large wrinkle amplitude can keep the local compressive strain reduced below the Biot critical strain at the valley of the wrinkle.

The most notable point on the primary ruga phase diagram is the M-point in figure 5. As discussed above, the M-point splits the Biot bifurcation curve into two parts, the unstable instantaneous crease curve and the stable progressive-wrinkle curve. This M-point splitting has an important implication on creasing of neo-Hookean EDM solid surfaces, in general. It defines the range of applied compressive strain 0.272<ε<0.456 that can instantaneously crease the surface. The branch curve C13 in figure 5 implies that instantaneous creasing sets off when k/a and ε have values on the curve. Once the surface creases, the unloading process of the crease is irreversible. The crease configuration recovers the flat state only when the compressive strain is reduced down to the lower bound crease strain shown as a chain curve in figure 2a. By taking the limit a→0, or Embedded Image, we conclude that the lower bound of the crease strain converges to 0.35 for homogeneous neo-Hookean solids, which is consistent with the findings of Hong et al. [9] and Hohlfeld & Mahadevan [18].

Now, we consider what causes the nonlinear bifurcation to be stable on C12 and unstable on C13. Our energy calculation shows that the longer wavelength configuration has a lower energy level among admissible bifurcation modes with respect to the constraint of the boundary condition, typically the specimen size. Therefore, the longest-wavelength mode would be activated first, and the amplitude of the wrinkle would grow. Once the amplitude begins to grow, the wrinkle valley will encounter two competing effects. One is compressive strain relaxation owing to elongation of the cord length of the surface, and the other the local strain concentration caused by the new geometry of the surface. To follow these competing effects, we plot the local strain εl evolution against the applied compressive strain ε in figure 6 for isoperiodic compressions with different Embedded Image.

Figure 6.

Strain concentration diagram: horizontal and vertical axes represent, respectively, the applied nominal compressive strain and the local nominal compressive strain near the nucleation point of creasing. The solid and chained curves show traces of instantaneous creasing; the dashed, double-chained and long-dashed lines display loci of wrinkle-setback creasing; the coarse-dotted line exhibits a trace of fold-setback creasing. The mark (e) corresponds to the starting point of folding shown in figure 3e.

In figure 6, the 45° straight line indicates stable isoperiodic compression keeping the surface flat. In all cases of different Embedded Image, the lines initially turn upwards as they are bifurcated, meaning that every local strain at the valley of the wrinkle is amplified as soon as the surface buckles. The bifurcated curves have finite slopes soon after the bifurcations for Embedded Image; the steeper the slopes, the larger the Embedded Image is. However, the post-bifurcation slopes immediately diverge to infinity for Embedded Image implying instability. For stable bifurcations of Embedded Image, the strain concentration rates are reduced as they are further compressed. At this stage, the local strain relaxation rate caused by elongation of the surface cord length takes over the strain concentration rate owing to the geometry of the surface. Then, growth of the local strain at the valley of the stable wrinkle accelerates until the local strain reaches close to εB. As Hohlfeld & Mahadevan [18] and Cao & Hutchinson [10] noted in their FEM observations, the plane-strain crease-tip strain is believed to be the Biot critical strain εB. Regarding the curve for Embedded Image in figure 6, the local strain relaxation mechanism prevails, and the specimen can be compressed with its nominal strain much beyond the Biot critical strain limit. In addition, as described earlier, isoperiodic compression makes the progressive wrinkle eventually fold for Embedded Image. Once it folds, as shown in figure 3e and marked as (e) in figure 6, the strain concentration rate changes substantially, and creasing is delayed.

Returning to figure 5, the setback-crease curve C2c has two peculiar points. One is the T-point (ε=0.174, Embedded Image) at which the applied compressive strain εT takes the absolute minimum value 0.174 for setback creasing. While spontaneous setback creasing is possible in neo-Hookean EDM solids for εTε<εM, setback creasing is only possible with initial imperfections in homogeneous neo-Hookean solids. The other is the R-point (ε=0.593, Embedded Image); compression beyond this point folds the wrinkle. The folding is a characteristic of a thin stiff film attached on a soft substrate. This result shows that a compressive strain of more than 0.60 is required to fold a neo-Hookean EDM solid. Once the fold is further compressed, a setback crease sets off at the root of the fold. When we rescale the solid surface at the root of the fold, the renormalized effective compressive strain and the effective admissible minimum wavenumber satisfy the condition to sit on the instantaneous-crease curve C13 to set off the final setback creasing.

Moreover, it is found that, unlike a soft substrate with a stiff film, a neo-Hookean EDM solid does not activate period multiplication during its wrinkle growth because the critical wavenumber monotonically increases with the compressive strain. Creasing in a neo-Hookean EDM solid hardly develops period multiplication as well, although there is a setback-creasing branch along which the critical wavenumber decreases for increasing compressive strain. It is because setback creasing requires growth of a wavenumber of a preset wrinkle to trigger the unstable snap jump of creasing.

5. Conclusions

In this paper, we have explored creasing processes of a neo-Hookean EDM solid surface under isoperiodic plane-strain compression, with results leading to the construction of a primary ruga phase diagram. In constructing the phase diagram, our nonlinear bifurcation analysis revealed a critical point, the M-point (ε=0.272, Embedded Image), on the phase diagram. The M-point splits Biot's bifurcation curve [6] such that a flat surface wrinkles stably under a compressive strain in the range 0<ε<εM=0.272, whereas it creases unstably for εM<ε<εB=0.456. The strain at which the flat surface starts to wrinkle or crease depends on the longest-admissible perturbation wavelength normalized by the decay length of the exponentially graded stiffness.

Moreover, our finite-element post-bifurcation analysis uncovered that the local compressive strain at the valley of the wrinkle is always amplified initially once the flat surface buckles. However, as the compression progresses further, growth of the large wrinkle amplitude tends to relax the local strain due to elongation of the surface cord length. For instantaneous creasing with εM<ε<εB and Embedded Image, the initial local strain amplification mechanism dominates to escalate the strain concentration to the level of strain localization of creasing. On the other hand, for stable wrinkle bifurcation with 0<ε<εM and Embedded Image, the local strain relaxation mechanism leads to progressive growth of the wrinkle amplitude until further compression sets off setback creasing. In tracing the setback creasing on the ruga phase diagram, we also found another critical point, the R-point (ε=0.593, Embedded Image, beyond which the wrinkle folds before it creases. In addition, another limit point, the T-point (ε=0.174, Embedded Image defines the minimum compressive strain required to trigger setback creasing. The FEM analysis also revealed that regardless of its type, instantaneous or setback, creasing is always accomplished when the local crease-tip compressive strain reaches the Biot strain, ε=0.456.


This work was supported by the Korea Institute of Machinery and Materials, the Korea Institute of Science and Technology and the US National Science Foundation (DMR-0520651) through MRSEC at Brown University. The authors acknowledge that Eunice Kim in Classics introduced the Latin word ruga for this research.

Appendix A. Perturbation analysis

(a) The zeroth-order solution

The zeroth-order solution is the affine deformation of F11=λ, F22=1/λ and F12=F21=0 with zero traction along the boundary X2=0, for which the only non-vanishing stress components are Embedded ImageA1aand Embedded ImageA1b

(b) The first-order solution

The first-order incompressibility condition given by equation (2.4c) allows us to relate the first-order displacement components through a function ψ(1)(X1,X2) as Embedded ImageA2aand Embedded ImageA2bInserting equations (A 2) into the first-order expansion in (2.4a) and (2.4b), and eliminating the hydrostatic pressure by cross-differentiation, the system of partial differential equations becomes the following fourth-order homogeneous partial differential equation: Embedded ImageA3The general solution of (A 3) is obtained analytically as Embedded ImageA4where the coefficients Ci are determined only from the boundary conditions. Using equations (A 4) and (A 2), we can now determine the first-order displacement components as given in (2.5a) and (2.5b). The hydrostatic pressure is determined by the tractions in (2.3) with (2.5a) and (2.5b) and the traction-free boundary condition as Embedded ImageA5

(c) The second-order solution

Similar to the technique used to solve for the first-order displacement components, the second-order incompressibility condition is used to relate the second-order displacement components as Embedded ImageA6aand Embedded ImageA6bwhere Embedded ImageA6cand Embedded ImageA6dFollowing the same procedure as in the first-order solution, the system of partial differential equations becomes the following fourth-order inhomogeneous partial differential equation: Embedded ImageA7Then, (A7) admits the following solution: Embedded ImageA8where the exponents βm (m=1,…,4) are obtained from the exponents α1 and α2 by replacing Embedded Image by Embedded Image in equation (2.6a). In (A 8), the coefficients B3 and B4 are set to zero to ensure the solution decays to zero at infinity.

By inspection of the right-hand side of the differential equation in equation (A 7), we can deduce that the particular solution ψ(2)p(X1,X2) has the following form: Embedded ImageA9aThe coefficients Fi (i=1,2,3) are determined by inserting (A 9a) into (A 7) followed by matching the coefficients of the different eigenfunctions to obtain Embedded ImageA9band Embedded ImageA9cwhere Embedded ImageA9dwith υ1=λ4+1, Embedded Image and Embedded Image.

We then insert equation (A 8) into equations (A6a) and (A6b) to determine the second-order displacement components as (2.9a) and (2.9b) and the second-order hydrostatic pressure as (2.9c) given in the main text.

Enforcing the traction-free boundary conditions at the top surface of the half-space with the second-order Piola–Kirchhoff stress, we obtain the following two equations: Embedded ImageA10aand Embedded ImageA10bwhich constitute a system of linear equations for B1 and B2 given in appendix B. Normalizing the amplitude of the first-order displacement component Embedded Image to unity in (2.5b), and using (2.7), the constants C1 and C2 are determined as Embedded ImageA11aand Embedded ImageA11bThe second-order expansion of the first Piola–Kirchhoff stress is then given as Embedded ImageA12a Embedded ImageA12b Embedded ImageA12c and Embedded ImageA12d

Appendix B. H1(X2), H2(X2), B1 and B2

Embedded ImageB1a Embedded ImageB1b Embedded ImageB2a Embedded ImageB2b and Embedded ImageB2cwhere Embedded ImageB3aEmbedded Image3bwith Embedded ImageB4a Embedded ImageB4b Embedded ImageB4c Embedded ImageB4d Embedded ImageB4e Embedded ImageB4f Embedded ImageB4g Embedded ImageB4h and Embedded ImageB4i

  • Received December 28, 2012.
  • Accepted June 4, 2013.


View Abstract