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 Q0 to a bulk value QB with a decay length of 1/a. Then, the normalized perturbation wavenumber, , and the compressive strain, ε, control formation of a wrinkle pattern and its evolution towards crease or fold patterns for every stiffness ratio η=QB/Q0. 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 . Competition among film-buckling, local film-crease and global substrate-crease modes of energy release produces diverse ruga-phase domains. Our analysis also reveals the subcritical crease states of the homogeneous half space. Our results are, then, compared with the behaviour of equivalent bilayer systems for thin-film applications.
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 . 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. , Brau et al. , Sun et al.  and Hutchinson . 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  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  studied the same problem using the formulation developed by Green et al.  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 . 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.  also obtained this critical strain using a pinching perturbation scheme in their finite-element analysis (FEA). Cao & Hutchinson  showed, using Koiter's approach  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.  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  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.  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.  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, , larger than 1.3 engenders the second-order instability of wrinkling with its critical strain ε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 prompts instantaneous creasing instead of wrinkling, whereas wrinkling sets off for and 0<ε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 . By contrast, for , the wrinkle requires further compression to commence a setback crease or fold. Here, the question still remains whether the second-order wrinkle instability can guarantee instantaneous creasing in other neo-Hookean material systems with different modulus distributions or under different boundary conditions. This question is treated in the following sections.
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 2.1 where η is the ratio of the shear modulus of the substrate bulk material QB to the shear modulus at the free surface Q0, 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 Xi to a deformed configuration xi is denoted by 2.2a Here, λ(i) are the in-plane principal stretch ratios of the fundamental state aligned in the (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 2.2b for 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 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 ui(X1, X2) and the net hydrostatic pressure p(X1, X2) are expanded about the fundamental state with respect to a perturbation parameter ζ up to the second order 2.3a and 2.3b where p(0)(X1,X2) 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, . The superscripts in parentheses indicate the perturbation order. Using equations (2.2) and (2.3), the first-Piola–Kirchhoff stress 2.4 is obtained up to the second order as expressed in appendix A.
Then, the equilibrium equations up to the second order are reduced as 2.5a and 2.5b and the incompressibility condition up to the second order to 2.5c Then, similar to the procedures in Diab et al. , we first solve the linear bifurcation problem to obtain the critical wrinkle bifurcation strain and its corresponding normalized wavenumber, () for various afv-EDM systems within 0<η<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 2.6 Here, φ 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 2.7 The terms at the right-hand side of equation (2.7) are determined for a neo-Hookean EDM solid as 2.8a and 2.8b where i, j, k, m, n=1/2, and a linear operator means .
Equation (2.8a) exhibits an explicit coupling effect of and for the energy variation, whereas equation (2.8b) entails the effects of implicit coupling between the first- and the second-order fields. The first- and second-order fields are obtained by solving the system of ordinary differential equations in equation (2.5) as shown in appendices A and B, and in the electronic supplementary material.
Figure 2a 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 2b, 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.  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 plane. The bifurcation curves have knee shapes for afv-EDM systems (0<η<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, . For afv-EDM systems, the second-order energy variation of the wrinkle state with respect to the flat state is used to determine the stability of every bifurcation state along the upper ( and lower () branches of the linear bifurcation curve. The critical states of stable ip-wrinkles are delineated from those of unstable ip-wrinkles by the M point as shown in figure 3 for η=0.1. There is a cut-off frequency above which the ip-wrinkle bifurcation becomes unstable and is displayed by dashed curves. The 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 than that of 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 . At η=0.006, an infinitesimally small critical zone appears at a point , in which the first-order wrinkle mode is unstable. The upper and lower limits of this critical zone are denoted by M′ and M′′ points. As η increases beyond 0.189, instantaneous creasing commences for all above 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 , perturbed deformation near the surface is likely to induce unstable creasing. However, as η 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 . The chained line in figure 3 represents the trace of the boundaries of the unstable wrinkle bifurcation states as a function of η as determined by the energy variation in equation (2.7). As η increases towards unity, i.e. a homogeneous half space, the critical of the limit point 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 is varied by changing the decay length of the temperature distribution. The quadratic plane-strain hybrid element with reduced integration (CPE8RH) is adopted to simulate the incompressible neo-Hookean material. As shown in figure 4a, 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, , at a stress free reference state. We used throughout our simulations. Details of generating the perturbed mesh are given by Diab et al. . 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 4b–f show nominal compressive strain (−E11) distributions in various ruga states created by isoperiodic compression at different 's for an afv-EDM with η=0.001. The surface stably wrinkles with a small ip- as shown in figure 4b . Upon further compression, the stable wrinkle can unstably crease to create a setback crease state as shown in figure 4e . Herein, the creasing strain is detected by evaluating the growth of the Fourier-series coefficients of the deformed surface that was shown by Diab et al.  to diverge at the creasing strain. The critical strain for the setback crease is computed and plotted in figure 5 as a function of for various afv-EDM systems under isoperiodic compression, as dashed curves below the 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- as shown in figure 4d . The trajectories of instantaneous crease states for various η 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 for 0.006≤η≤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 , and found that there is a critical wavenumber below which wrinkles may fold depending on η as shown in figure 4c for η=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 4f . For η>0.002, our FEM simulations show that wrinkles may never fold for any admissible wavenumber ip- and the crease curve makes a turn at low value as shown in figure 5.
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  or through unloading the critically or snap-bifurcated crease state . 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-. The arrows on the hysteresis loop of the loading and unloading path indicate the direction of the unstable snap jump of creasing at the critical strain ε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 . The plots show that the creasing is a typical subcritical pitchfork bifurcation, whereas wrinkling from a flat state is a typical supercritical pitchfork bifurcation. It is interesting to note that the lower-limit strain εL is still smaller than εM=0.272 at , indicating that the lower-limit strain of the subcritical state of the setback crease is smaller than the critical strain of wrinkling for slightly smaller than . It is well known that the wrinkle amplitude continuously grows following a typical supercritical bifurcation branch as the compressive strain increases beyond the critical bifurcation point. By contrast, creasing trails its stable branch of the bifurcation path after the amplitude of crease suddenly jumps at the critical bifurcation point. The jump is typically observed as a dynamic snap jump. In figure 6, all the stable branches of the bifurcation paths are represented by solid curves and unstable branches by dashed curves, whereas all the jumps are denoted by dotted lines. The results show that, as becomes large, the stable branch of the subcritical bifurcation merges to that obtained by Hohlfeld & Mahadevan  and Hong et al.  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 and ε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 . However, the height of the energy barrier from the flat to the crease states is positive 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 .
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=Qf/Qs, is conducted, as shown in appendix C. Here, Qf and Qs 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 , but with the factor β) 4.1 and 4.2 where β, and are, respectively, a geometrical factor of equivalent thickness, the effective film thickness and modulus of the graded boundary layer of the afv-EDM system. Using equation (4.2), the effective modulus ratio of the afv-EDM system, , is related to the modulus ratio of the bilayer system as 4.3 Our first-order stability analysis indicates that the effective boundary layer thickness of the afv-EDM system with β=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 , whereas much of the excessive energy is supplied by the substrate 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 , respectively, in figure 7a,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 . The afv-EDM system has the equivalent stiffness ratio (ηm=0.189), whereas the bilayer system takes the stiffness ratio Rm=1.71. In figure 3, it is shown that the critical creasing of an afv-EDM half space of has a self-selective wavelength and is an instantaneous creasing from a flat state. Regarding the self-selectivity of the wavelength, the analysis in §2 showed that self-selective wrinkle bifurcation becomes unstable to generate periodic instantaneous creasing for the afv-EDM system with η 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 7a shows that the critical strain of the afv-EDM system at the M point is almost invariant, εM≈0.272, for 0≤η≤0.1 () and decreases slightly down to εM=0.258 as η increases (i.e. decreases) within . However, for the bilayer system, εM increases monotonically from εM=0.187 (at ) to εM=0.33 (at Rm=1.71). Each pair of the dotted and chained curves in figures 7a,b represents the traces of the εM′′ and kM′′h for the afv-EDM and the bilayer systems as functions of and R, respectively. As noted earlier, for and 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 7a, 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 7b; 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 or of relatively low layer stiffness (η→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.
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 Q0 to the bulk value QB 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<η≡Q0/QB<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, , triggers cascade instability of infinite order manifested by localized creasing at the free surface, for 0≤η<1. The lower-limit wavenumber of creasing, , is greater than the self-selective wavenumber for 0<η<0.189, and jumps down to for 0.189≤η<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 for 0.006≤η<0.08. On the other hand, the wrinkling develops cascade instability of infinite order, i.e. creasing, within a finite wavenumber range, , for 0.08≤η<0.189. At η=0.189, merges with to disappear, and the lower-limit creasing wavenumber jumps down to . The parametric study also shows that self-selective periodic instantaneous creasing develops directly from a flat state for 0.148≤η<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 for a vai-EDM system. Taking the latter approach, the subcritical strains and the depths of creasing for are identified for a vai-EDM system by FEM simulations of loading–unloading cycles of the system with various . In this analysis, the lower limit of the subcritical creasing strain of the homogeneous half space, εHMHS=0.35, is obtained by taking . Our simulations also reveal an unstable subcritical bifurcation barrier hidden between the high-energy flat equilibrium state and the stable low-energy subcritical crease state for 0.35<ε<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.
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 A 1a A 1b A 1c A 1d where the traction-free boundary conditions along X2=0 is satisfied at the zeroth order by setting A e The first-order incompressibility condition given by (2.5c) allows us to relate the first-order displacement components through a function ψ(1)(X1,X2) as A 2a and A 2b Inserting (A 2a) and (A 2b) into the first-order expansion in (2.5a) and (2.5b), and eliminating the hydrostatic pressure by cross-differentiation, the system of partial differential equations becomes the following fourth-order homogeneous partial differential equation A 3 Here, we are interested only in solutions where the normal displacement at the surface is symmetric about X1=0, and thus we write A 4 Using (2.5a), the first-order hydrostatic pressure is determined as A 5a where A 5b Then, we insert (A 4) into (A 3) to obtain, A 6 The fourth-order differential equation in (A 6) is written in a matrix form as A 7a where A 7b and A 7c with Using (A 1c) and (A 1d), the traction-free boundary conditions at the top surface of the EDM system is written in a matrix form as A 8a where A 8b The displacement boundary conditions at the bottom surface of the EDM system in the perturbed configuration are written in a matrix form as A 9a where A 9b The eigenvalue problem in (A 7)–(A 9) can be solved analytically for η=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 B 1a and B 1b The second-order hydrostatic pressure is then determined by (2.5b) as B 2a where B 2b B 2c B 2d with .
Similar to the first-order solution in appendix A, the second-order perturbation problem is cast as B 3a where B 3b B 3c with B 3d
The components of the matrix A(2)(X2) are obtained by replacing by in equations (A 7).
The second-order traction-free boundary conditions at X2=0 and the second-order zero displacement boundary conditions in the perturbed configurations at are given, respectively, as B 4a and B 4b where C(2)=C(1) and B(2) is obtained from B(1) by replacing with and B 4c B 4d Finally, the first- and second-order eigenfunctions are obtained by solving recursively the system of differential equations in (A 7) and (B 3) and their corresponding boundary conditions in (A 8) and (B 4) using the standard two-point boundary value problems solver in Matlab (bvp5c).
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)(X2) and Ψ(2)(X2) are obtained as C 1a and C 1b where C 1c C 1d C 1e C 1f C 1g Using (A 2), (A 4), (A 5), (B 1) and (B 2), the first- and second-order displacements components and hydrostatic pressure are given as C 2a C 2b C 2c C 2d C 2e C 2f where C 2g C 2h C 2i C 2j C 2k and Γ(X2), E(X2) are deduced from (B 2c) and (B 2d).
The following first-Piola–Kirchhoff stress components are obtained up to the second order as follows C 3a C 3b where C 3c C 3d and C 3e C 3f The analytical expressions presented in this appendix are used to solve the nonlinear bifurcation problem of the bilayer system. Analytical solution of the linear bifurcation problem and closed form analytical solutions of the constants in the first- and second-order perturbations solutions are given in the electronic supplementary material.
- Received March 17, 2014.
- Accepted May 21, 2014.
- © 2014 The Author(s) Published by the Royal Society. All rights reserved.