The Primary Bilayer Ruga-Phase Diagram II: Irreversibility in Ruga Evolution

Recent studies [1–9] on Ruga patterns, such as wrinkles, creases, ridges, and folds, directed us to construct, in our previous paper [10], a comprehensive PB-RPD of a neo-Hookean bilayer system under monotonic compressive loading on the $k¯$ – $ε$ plane. Here, $k¯=3μs/μf3$ denotes a normalized characteristic wavenumber, which represents a scaled relative shear modulus of the substrate (⁠$μs)$ to the film (⁠$μf)$⁠, and $ε$ is the compressive strain. We showed on the PB-RPD that the system develops various stable equilibrium configurations of Ruga-phases, depending on $k¯$ and $ε$⁠. Therein, we analyzed evolution paths of various Ruga-phases under gradually increasing compressive loading starting from the flat state to the ultimate global localization. The PB-RPD is expected to serve as a guiding map to engineering multifunctional Ruga materials which have a wide range of applications, such as stretchable electronics, soft robotics, and bio-inspired adhesion devices [11–20].

However, in practical applications, materials normally experience repetitive cycles of loading and unloading, for example, charging/discharging of dielectric elastomer actuators in soft robotics and mechanical or thermal cyclic loading of stretchable electronic devices. During these cyclic loading processes, surface morphologies and material response characteristics can change significantly. Therefore, it is crucial to evaluate irreversibility of the Ruga-phases on the PB-RPD. Recently, several studies have revealed that particular Ruga-phases are irreversible during unloading. For instance, Hohlfeld and Mahadevan [21] and Diab and Kim [22] analyzed the irreversibility of free-surface creasing of a neo-Hookean solid, and Zhao et al. [23] reported the irreversibility of certain multimode wrinkles. In this paper, we focus our attention on a broad spectrum of cyclically irreversible Ruga-phases by investigating configurational irreversibility of various Ruga-phase transitions on the PB-RPD. The configurational irreversibility of a Ruga-phase during a loading and unloading cycle implies loss of static equilibrium configuration typically at the critical points of either unstable crease localization or transition in the mode of wrinkle-Ruga configuration. We will simply call the former as “localization irreversibility” and the latter “modal irreversibility.”

An example of localization irreversibility is Ruga transition of a smooth free surface or interface to or from a state of a crease tip under subcritical or critical far-field strain. Here, we define the critical point as the instance at which the smooth surface snaps to form a crease tip, e.g., creasing of a flat free surface at the Biot strain, 0.46. Such a snapping transition of crease formation typically leads to a load–deformation hysteresis with crease-phase mode locking in reverse loading. Similar localization irreversibility is also observed in ridging. In contrast, two different types of modal irreversibility are observed in Ruga-phase transitions among the single-, double-, quadruple-, and folding-mode wrinkles. One is the rate-independent modal irreversibility caused by unstable snapping transition of a large-amplitude wrinkle mode to another large-amplitude wrinkle mode [24]. The other is the rate-dependent modal irreversibility triggered by inhomogeneous viscoelastic deformation of a wrinkle, which is typically promoted by local pop-out snapping of shallow valleys of the wrinkle. Such snapping switches the primary period of the wrinkle.

While the nonlinear finite deformation of Ruga-phases sets the energy landscape for rate-independent bifurcation of static equilibrium states, various modes of rate-dependent incremental deformation can trigger transitions from a Ruga-phase to various metastable Ruga-phases. Both types of modal irreversibility are also observed in ridging and crumpling, and the modal transitions are in general sensitive to lateral boundary conditions in a finite size specimen. Regarding rate-dependent irreversibility, while viscoelasticity effects in wrinkling have been studied extensively [25–28], viscoelasticity effects in general Ruga evolution have not been well understood yet. In this paper, we show that viscoelasticity promotes generation of diverse modes of metastable Ruga configurations, which instigates irreversibility of Ruga evolution in cyclic loading. In turn, it is found that cyclic instability or stability, i.e., divergence or shake-down, of Ruga configurations can be controlled by viscoelasticity. The rate-dependent irreversibility of Ruga configurations can be a critical behavior of bio-organ deformation [29–31].

This paper is organized as follows: In Sec. 2, we will highlight generic irreversible Ruga-phases on the PB-RPD, employing the finite-element method (FEM) model developed in Part I [10]. Here, a generic irreversible Ruga-phase stands for the phase created by rate-independent bifurcation under quasi-static compression. Throughout this paper, we employ negligible incrementally linear viscoelastic damping [32,33] by setting the loss tangent (damping energy/elastic energy) $2×10−5$ to stabilize the incremental equilibrium-configuration search algorithm of FEM for various bifurcation processes, except for studying specific viscoelastic effects in Sec. 4. In Sec. 3, irreversible transitions of multimode wrinkle phases are presented in a PB with no mismatch strain. Then, we will discuss enhancement of Ruga irreversibility triggered by viscosity effect in multimode wrinkle transitions of PB with strain mismatch in Sec. 4. Finally, in Sec. 5, possible irreversible evolution pathways made by cyclic loadings and their corresponding morphological features are sorted and discussed.

Figure 1(a) shows the generic irreversible Ruga-phases represented as colored zones on the PB-RPD reported in Part I. While the wrinkling transition across the boundary T₁₂ in Fig. 1(a) is reversible during a loading–unloading cycle, irreversibility is observed in some other transitions; the irreversibility depends primarily on the compressive strain and the film/substrate moduli ratio. The zones of various Ruga-phases are denoted by Roman numerals on the PB-RPD as follows: Flat state (I), primary single-mode wrinkle [SM1; (II)], wrinkle-setback crease (IIc), global crease limit (III), double-mode wrinkle [DM1; (IV)], double-setback crease (IVc), quadruple-mode wrinkle [QM1; (V)], quadruple-setback crease (Vc), fold-setback crease (Vcf), fold (VI), and global fold localization limit (VII). The number 1 of multimode wrinkles, in SM1, DM1, and QM1, signifies the primary period of the mode being the lowest-energy wrinkle period of the PB. In the following, metastable multimode wrinkles with their primary periods being multiples of the lowest-energy wrinkle period are symbolized as SMi, DMi and QMi, i = 2, 3,…. In specifying the period, we distinguish two different types of Ruga characteristic period, primary period (P-period) and mode period (M-period). Their definitions are pictorially illustrated as P-period of single-mode [SM1; (II)] in Fig. 1(b1), M-period of double-mode [DM1; (IV)] in Fig. 1(b2), and M-period of quadruple-mode [QM1; (V)] in Fig. 1(b3). In M-period transitions under large compressive strain, its P-period remains the same as that of the primary wrinkle configuration.

As previous studies [21,22] revealed, an individual crease under cyclic loading makes a smooth free surface (or interface) snap-jump to and from subcritical crease states. As aforementioned, a flat free surface under monotonic compressive loading snaps to crease at the Biot critical strain. In a special case, a surface defect can cause singular perturbation on the fundamental state and trigger reversible creasing on the flat free surface of a neo-Hookean half space at a smaller compressive strain, approximately 0.35 [34]. However, in PB, all the crease tips along the curves of T₁₃, T_2c, T_4c, and T_5cf are under corresponding critical far-field strains and follow trajectories of subcritical crease states with mode locking upon unloading. In turn, the transitions across the curves are all irreversible under cyclic loading. The boundaries of crease unloading jump are located on the left side of the critical crease loading jump curves in PB-RPD (not shown in Fig.1(a)). These irreversible transitions make the Ruga-phases, IIc, IVc, Vc, and Vcf (in green) and those toward (III) (in pink) locally irreversible for their crease tips. Since the local irreversibility of creasing is relatively well understood by now, primarily the modal irreversibility of bilayer creasing and multimode wrinkling will be treated in the rest of this paper.

Figure 2(a) shows the variations of the normalized crease-depth difference, ΔΩ, between near-neighbor creasing valleys in four P-period as a function of the nominal compressive strain, ε, for $k¯=0.59$⁠, across the initiation point of global creasing at $ε≈0.35$ during a loading–unloading cycle. Throughout this paper, ΔΩ is normalized by the critical onset-wrinkling (P-period) wavelength. Solid curve corresponds to loading and dashed curve to unloading, exhibiting apparent hysteresis. When the loading varies, the surface configuration maintains its translational symmetry through the SM1 wrinkle phase for $εw<ε<εd$⁠, where $εw$ and $εd$ represent the critical strains of wrinkling and period doubling, respectively. Once the period doubles, the DM1 wrinkle steadily grows upon further compression until $ε$ reaches the period quadrupling strain $εq$⁠, where every four other deep valleys crease, i.e., $εq=εc$⁠. Beyond this crease initiation, the periodic crease depths grow uniformly as shown in Fig. 2(b)- $ε1l$ until the compressive strain approaches an apparent bifurcation strain, ε₁. At this strain, periodic creases begin to deviate from the translationally symmetric configuration. Upon further compression, ΔΩ follows the upper bifurcation branch (solid line) in Fig. 2(a), until the compressive strain reaches creasing strain of the substrate, 0.35, beyond which ΔΩ grows rapidly, developing global crease localization. Upon unloading, traces of ΔΩ closely follow the solid branch until the strain reaches a value near 0.35 at which the system jumps down to the lower parity bifurcation branch (dashed line).

Figure 2(b) (⁠$ε2u)$ shows the parity bifurcation configurations observed during a loading–unloading cycle across the critical 0.35 strain. As shown in the insets of $ε1l$ and $ε2l$ in Fig. 2(b), symmetry of the cooperative crease mode is broken beyond $ε1$⁠, and ΔΩ grows rapidly across the 0.35 strain to reach the configuration of $ε3l$⁠. Upon unloading, variation of the configuration from $ε3u$ to $ε2u$ changes the sign of ΔΩ in $ε2u$⁠, compared to that of $ε2l$⁠. Then, further unloading recovers the periodic configuration of $ε1u$ same as $ε1l$⁠. We believe that this parity bifurcation is primarily caused by fictitious viscosity that we employed in our FEM model to stabilize convergence of the numerical solution. In the limit of no viscosity, the two traces may merge and transit vertically at the compressive strain of 0.35.

Regarding the initial crease formation, the critical strain, ε_c, strongly depends on the critical wavenumber $k¯$ in the range of $0.40<k¯<1.44$⁠; on the other hand, the transition of the cooperative film-crease mode toward the global crease localization of the substrate is consistently triggered at a fixed strain value of $ε≈0.35$ for $0.40<k¯<1.20$⁠. Starting from the homogeneous half space (⁠$k¯=1.44$⁠), for $1.20<k¯<1.44$⁠, the flat surface of the film develops instantaneous creases through transition T₁₃. For $0.40<k¯<1.20$⁠, an SM1 wrinkle develops periodic local setback creases on the free surface of the film in every four other wrinkle valleys. The setback creasing transitions are denoted as wrinkle-setback crease across T_2c for $1.00<k¯<1.2$⁠, double-setback crease across T_4c for (⁠$0.60<k¯<1.00$⁠), quadruple-setback crease across T_5c, and the fold-setback crease generated by further loading across T_5cf for $0.40<k¯<0.60$⁠. All these transitions lead to uniform stable growth of periodic film crease tips up to 0.35 compressive strain for $0.40<k¯<1.20$⁠, crossing the green zone in Fig. 1(a), followed by global localization of substrate creasing beyond the 0.35 strain (pink zone in Fig. 1(a)). Our FEM simulations reveal that the periodic film creasing occurs through unstable snap jump at a critical strain during loading cycle, while the crease phase locks in its mode across the critical strain during unloading cycle. Then, the periodic crease tip configuration maintains its mode down to a subcritical strain at which the periodic crease phase snaps back to the wrinkle phase [21,22,35]. In contrast, the global creasing exhibits a different irreversibility phenomenon—modal instability in growth or recession of the periodic crease tips. The modal irreversibility of creasing is often amplified by viscosity effect on the interactions among the film creases and the substrate crease.

Besides the above-discussed creasing irreversibility, we have also investigated irreversibility of multimode wrinkling transitions: M-period doubling (II–IV) through transition T₂₄; M-period quadrupling (IV–V) through T₄₅; folding (V–VI) through T₅₆; and fold localization (VI–VII). Our FEM results show that reversibility of the various modes (DM1, QM1, and F) depends mainly on the modulus ratio of the bilayer system, $R (=μf/μs)$⁠. Therefore, we have run FEM simulations of the loading–unloading cycle for the whole range, $0<k¯<0.40$ (⁠$∞>R>45)$⁠, within which M-period multiplications occur.

Our simulations reveal two distinct characteristic wavenumbers, $k¯=0.29 (R=120)$ and $k¯=0.17 ( R=600)$ denoted as $QIR$ and $DIR$⁠, respectively, in Fig. 1(a). Within $0.29<k¯<0.40$⁠, i.e., above $QIR,$ the system undergoes reversible transitions from flat state (I) to fold (VI) through multimode wrinkling phases. We only observe pronounced rate-dependent irreversibility when folding localizes in (VII), which is triggered by creasing of the substrate. In other words, transitions from SM1, DM1, QM1, and F within this characteristic wavenumber region are all reversible with no hysteresis. This particular reversible corridor of $0.29<k¯<0.40$ and $0<ε<0.30$ allows controllability of the fold patterns. This feature is expected to be technologically very attractive for various applications, such as stretchable electronics and soft robotics [5,36]. For $0.17<k¯<0.29$⁠, M-period doubling is reversible, while M-period quadrupling and folding processes are irreversible. For very stiff film bilayers, $k¯<0.17$⁠, the transitions of M-period multiplication, $II→IV→V→VI→VII$⁠, are all irreversible, which engenders hysteresis in Ruga evolution under cyclic loading, with mode locking in unloading cycle. The irreversible multimode wrinkle phases are marked in blue, brown, and purple on the PB-RPD in Fig. 1(a). In Secs. 3, 4, and 5, detailed analyses of the irreversibility and their implications on the formation and evolution of various modes are presented.

Transitions between various pairs of Ruga modes occur through symmetry breaking of height variations in the wrinkle configuration. Therefore, we have traced the evolution of the height differences between a wrinkle valley and its two adjacent valleys, denoted by $ΔΩl$ and $ΔΩs$⁠, respectively, as shown in the insets of Fig. 3. Here, the superscripts $l$ and $s$ denote large and small, respectively. Figure 3(a) shows traces of $ΔΩl$ and $ΔΩs$⁠, during a loading–unloading cycle for a bilayer system with the stiffness ratio R = 90 (⁠$k¯=0.32$⁠). During loading, the SM1 wrinkle symmetry is conserved until the critical strain for period doubling is attained at $ε≈0.18$⁠, at which the SM1 symmetry is broken and a DM1 configuration emerges. Then, another higher-mode bifurcation occurs at a strain (⁠$ε≈0.26$⁠) at which the DM1 bifurcates into a QM1. As shown in Fig. 3(a), both traces of $ΔΩl$ and $ΔΩs$ coincide for loading (solid line) and unloading (dashed line), indicating reversibility of period doubling and quadrupling for the bilayer system with R = 90. However, for a bilayer system with a high stiffness ratio, R = 3000 (⁠$k¯=0.10$⁠), the traces of $ΔΩl$ and $ΔΩs$ in Fig. 3(b) do not coincide for loading (solid line) and unloading (dashed line), showing hysteresis of the Ruga configuration in the cyclic loading. This configurational irreversibility is clearly shown in Fig. 3(b) with the simulation frames taken at the same strains, $ε1,…, ε4$⁠, of loading and unloading. Although the small fictitious viscosity employed in our FEM model to stabilize convergence obscures the sharp jumps at the doubling and the quadrupling points of the mode, the distinct hysteresis loops are evident for a highly stiff-film PB without strain mismatch. The results indicate that the unloading transitions undergo mode locking as shown in the frames captured at particular strains $ε1$⁠, $ε2$⁠, $ε3$⁠, and $ε4$⁠. Now, questions arise: Does irreversibility always manifest itself through mode locking or it may lead to excitement of other modes that are not accessible during loading?

In Secs. 2 and 3, we showed that the synergistic effects of viscosity and large nonlinear deformation can substantially enhance irreversibility of Ruga-mode transitions in PB without strain mismatch. Here, in this section, we present the excitation of P-period switching jumps, triggered by minute increase in viscosity, that further enhances irreversibility in cyclic Ruga evolution of PB with strain mismatch. In our simulations, the strain mismatch is implemented by prestretching the substrate and subsequently attaching a stress-free layer on the surface of the substrate. Upon releasing the prestretch of the substrate, the film is compressed to make the flat state unstable and SM1 wrinkle mode emerges. The FEM-plot insets in Figs. 4(a1) and 4(b1) show typical multimode wrinkle configurations, DM1 and QM1, observed during loading cycle of a bilayer system with the stiffness ratio R = 10⁴ (⁠$k¯=0.07$⁠) and the substrate prestretch ratio λ_ps = 1.4. In the graphs of Fig. 4, the horizontal axis represents nominal compressive strain of the film, measured with respect to the prestretched configuration of the substrate. The vertical axis $ΔΩ$ denotes the depth difference between the neighboring shallow and deep valleys,normalized by the P-period wavelength. If the prestretch for the strain mismatch was greater than approximately 1.5 for the stiff film, ridging instability would trigger different types of irreversible Ruga evolution such as order–disorder transitions in ridging and crumpling [37,38].

During compressive loading on the film, the bilayer system develops Ruga-phase evolution (SM1 → DM1 → QM1 → F) similar to the case of high stiffness ratio PB with no strain mismatch, as shown with solid curves in Figs. 4(a1) and 4(b1). However, the critical strains for period doubling and quadrupling, 0.34 and 0.44, are, respectively, larger than their corresponding strains, 0.18 and 0.26, for the same stiffness-ratio bilayer system without strain mismatch. For this stiff film, the transitions from single-mode (SM1) to double-mode (DM1) and from double-mode (DM1) to quadruple-mode (QM1) both exhibit unstable bifurcation of snap buckling, similar to those of the same stiffness-ratio bilayer system with no strain mismatch. However, unlike the latter without strain mismatch which makes mode locking during unloading cycle, the former with strain mismatch exhibits P-period switching jumps promoted by viscosity upon unloading as described below in detail. In order to examine the sensitivity of unloading bifurcation on viscosity, in Fig. 4, we present unloading Ruga-phase evolution stimulated by a minute but not negligible viscoelastic loss tangent of $2×10−4$⁠.

The dashed lines in Fig. 4(a1) show traces of ΔΩ transitions, while the compressive strain of the film is unloaded from a state of DM1 at a strain very close to the mode-doubling critical strain. The Ruga configuration reversibly follows the (green) solid curve between the critical points of mode doubling and mode quadrupling, until it is unloaded to a strain within 1% from the mode-doubling critical strain, as marked with * on the solid curve in Fig. 4(a1). At this point, it starts to make a snap buckling transition to a completely new mode, SM2, which is a SM but it has its period twice the period of SM1. The P point in Fig. 4(a1) denotes the moment of local film buckling at which the shallow valley transits to a peak. At the transition point P, shallow valleys turn flat to have zero curvature locally. Subsequently, they bulge out and become peaks. The unloading process after P point is depicted in the dotted box, where ΔΩ represents the depth difference between the deep valley and the newly formed peak of the SM2 mode. Figure 4(a2) shows the variation of the total strain energy of the PB along the loading/unloading process near the mode-doubling point, which corresponds to the gray zone in Fig. 4(a1). Here, the total strain energy is denoted by W and the strain energy at the starting point (o on the curve) of P-period switching by W_(o). When the prestretch of the substrate, $λps=1.4,$ is fully relaxed, the film reaches its compressive strain at $ε≈0.286.$ Therefore, the substrate is compressed, and the total strain energy increases monotonically beyond $ε≈0.286$⁠. Surprisingly, it increases smoothly along the blue curve in the range of $0.32≤ε≤0.355$ in Fig. 4(a2) despite the transition from SM1 to DM1. When the PB is unloaded from DM1, the Ruga mode and the corresponding strain energy begin to deviate, at $ε≈0.353$ (marked by a black circle), from those of the solid blue curve to those of the dashed black curve in Fig. 4(a2), making P-period switching to SM2. As the shallow valley pops out at P, the dashed energy curve exhibits inflexion toward that of SM2. The energy level of SM2 is clearly higher than that of SM1, indicating that SM2 is a metastable state.

Similarly, in unloading the system from QM1, shown in Fig. 4(b1), the system experiences a snap buckling transition from QM1 to DM2 as the strain approaches the critical quadrupling strain within 1%, $G+$ and $R+$⁠. The change in Ruga configuration at the P point in Fig. 4(b1) is similar to that in Fig. 4(a1); shallow valleys transit to peaks. Upon further unloading, another transition from DM2 to SM2 is observed. Figure 4(b2) shows variation of the strain energy along the loading/unloading process near the mode-quadrupling point within the strain range of the gray zone and the dashed-line box in Fig. 4(b1). The dashed-line energy variation curve shows that DM2 is indeed a metastable state. The FEM simulation configurations of the snap buckling transitions (SM1 $→$ DM1 $→$ QM1 $→$ DM2) are exhibited in the frames (c1–c4) of Fig. 4(c). This unloading bifurcation process of QM1 $→$ DM2 is made in competition between the local bifurcation process of shallow valley pop-ups, $QM1→DM2$⁠, and the modal bifurcation process of the long-range interaction mode, $QM1→DM1$⁠. In this competition, a small fluctuation of the configuration evolution caused by material inhomogeneity or viscosity can bias the bifurcation. In particular, it is found that material viscosity plays a significant role to bias the bifurcation toward $QM1→DM2$⁠. Once QM1 bifurcates to $DM2$⁠, the amount of unloading strain to reach $SM2$ is much smaller than the difference between the critical strains of $SM1→DM1$ and $DM1→QM1$ transitions. Then, upon further unloading, SM2 maintains its mode in a metastable state. Although SM2 is a higher energy mode than SM1, the energy barrier from SM2 to SM1 by doubling the frequency is not negligible and the SM2 sustains the mode as a metastable state. Direct doubling of the frequency from SM2 to SM1 for every crust has very high activation energy barrier, while frequency doubling through snaky adjustment by propagation of wrinklons [39–41] is considered to have lower energy barrier. However, the latter energy barrier is considered to depend on the lateral boundary conditions and the specimen size.

We discussed in Sec. 3 that mode locking in reverse loading near a critical point can instigate Ruga irreversibility, and in Sec. 4 that P-period switching jumps can dramatically boost cyclic irreversibility of Ruga evolution. In this section, we examine how viscosity triggers P-period switching jumps to enhance the cyclic irreversibility, and what are the consequences of the mode-locking and the P-period switching mechanisms on attaining cyclic stability of Ruga configuration hysteresis.

As discussed in Secs. 3 and 4, different bifurcation mechanisms such as mode-locking and P-period switching can compete near a certain critical strain during unloading cycle. Bifurcation of an autonomous local deformation mode, e.g., local curvature flipping of shallow wrinkle valleys, can occur while a transition is made in a long-range interaction cooperative mode, e.g., DM or QM, etc. Then, such an autonomous local mode bifurcation can substantially enhance irreversibility of cooperative Ruga modes in cyclic loading of PB with strain mismatch. In particular, material viscosity can delay relative evolution speed of a certain mode to trigger an unexpected autonomous mode which would not be activated if the materials were purely elastic. Here, we show that such an autonomous mode bifurcation can trip cyclic instability of Ruga evolution, leading to cyclic divergence or cyclic shake-down of Ruga modes.

Figures 5(a) and 5(b), respectively, show nearly reversible unloading and loading states of surface velocity V and configuration Ω of a DM1 Ruga surface of a low viscosity system with $2×10−5 loss tangent$⁠, while Figs. 5(c)–5(e) display those of irreversible unloading process in a relatively high viscosity system with $2×10−4 loss tangent$⁠. Here, the surface velocity V is defined as change in vertical surface displacement normalized by the P-period, with respect to global strain increment. The configurational deformation amplitude Ω is also normalized by the P-period. As illustrated in Fig. 5(a) for unloading of the low viscosity system, the shallow valley is moving upward while the neighboring deep valleys are moving downward, exactly opposite to the loading case (Fig. 5(b)). For the low viscosity system, materials are considered to experience nearly reversible deformation. In contrast, the unloading surface velocities shown in Figs. 5(c)–5(e) for the high viscosity system are out-of-phase with respect to the velocity shown in Fig. 5(a). The criticality of the velocity bifurcation is noticed in Fig. 5(d); the velocity of the shallow-valley surface is just reversed to that in Fig. 5(a). Figures 5(e) and 5(f) show that upon further unloading the shallow valley snaps out while its adjacent valleys reverse their motion albeit at smaller speed. This local snapping is believed to be caused by viscosity effect for different deformation rates along the corrugated surface.

Figure 6 collectively shows the schematic configurations of various Ruga-phases that can appear on the surface of a PB system under repeated cycles of loading and unloading. The three major PB Ruga localization types are illustrated in Fig. 6(a)—the crease, the fold, and the ridge localizations. The PB crease localization is a global localization. It commences at the Biot strain (0.46) on a smooth surface of a homogeneous incompressible neo-Hookean half space, and at 0.35 strain for the half space with a surface imperfection. Such an imperfection includes single or multilayer thin films, graded modulus distribution, surface singularity, and defects. The PB fold localization is induced by substrate crease localization, but the onset strain widely varies depending on self-contact conditions of local folds elicited by the folding process. For a PB system with $0<k¯<0.60$⁠, the surface layer generates periodic local folds at every fourth wrinkle valley, but the critical strain of global fold localization varies depending on smoothness of the local fold fronts. (i) When $0.40<k¯<0.60$⁠, each local fold generates a crease tip on its advancing front. The crease tip tends to stabilize cooperative parallel growth of the local crease tips until the substrate globally creases at the compressive strain of 0.35. (ii) When $0<k¯<0.40$⁠, the local folds grow with a smooth fold front, which easily induce instability with high imperfection sensitivity in cooperative growth of the local folds. Such instability can occur at any compressive strain between the local-folding and the Biot strains. These two kinds of folding localizations—crease tip guided localization and smooth fold localization, have different characteristics in configurational energetics. Crease tips on a fold front make lateral configurational sliding of the tips difficult, in contrast to smooth fold fronts. Easy lateral configurational sliding of smooth fold fronts causes high imperfection sensitivity in triggering global fold localization and thus large uncertainty in the critical strain of global fold localization. The PB ridge localization is a growth-limiting localization which sprouts out instead of sinking into the substrate. It is usually caused by strain mismatch between the film and the substrate, and growth of ridge height is intrinsically limited.

For a PB system of $k¯<0.40$ without strain mismatch, no crease manifests on the film surface, and the Ruga-phase evolves along the pathway (SM1 → DM1 → QM1) as shown in the bottom row of Fig. 6(b). For PBs with negligible viscosity, unloading follows the bottom row configurations (QM1 → DM1 → SM1) either with or without hysteresis, depending on the stiffness ratio. However, for an appreciably viscoelastic PB system of $k¯<0.40$ with strain mismatch, P-period switching makes the ground modes, DM1 and QM1 set off to SM2 and DM2, respectively, during unloading. Various possible P-period switching transitions are illustrated in Fig. 6(b). We denote these modes by SMi, DMi, and QMi, $i=1,2,3,…,$ for which modal translational symmetry holds in every single, double, and quadruple primary period(s). The index, i, stands for the ratio of the primary period of the mode with respect to the primary period of the ground mode. For example, the set-off modes SM2, DM2, and QM2 have the same geometric mode configurations as the ground modes SM1, DM1, and QM1, correspondingly; however, the primary period is doubled. When SM2 generated by P-period switching of DM1in the unloading cycle is reloaded by compression, the mode will bifurcate to DM2 instead of returning back to DM1. The energy barrier is too high to directly flip the convex curvature of the SM1 wrinkle peak to a concave curvature of a wrinkle valley. Therefore, it is likely that if the viscosity is properly matched for the loading rate, the wrinkle modes would diverge under cyclic loading, following SM1 → DM1 → SM2 → DM2 → SM4 → DM4. This mechanism can excite the ground mode (SM1, DM1, and QM1) to a higher energy modes (SMi, DMi, and QMi, i = 2 or 4,…) by tuning the loading rate for a given viscosity.

As discussed above, if the viscosity and the loading rate are properly matched, cyclic loading can make the Ruga-phase diverge to those of high-energy metastable configurations or shake-down to a certain energy-level metastable configuration. Furthermore, this mechanism can generate variety of coexistent Ruga-phases of PB with spatial imperfections in its property, geometry, or loading conditions.

We have revisited the PB-RPD [10] to investigate configurational reversibility of various Ruga-phases under cyclic loading. Our extensive finite-element simulations of single loading–unloading cycle reveal that many metastable Ruga states exist, and many phases (SM1, DM1, QM1, and F) on PB-RPD are irreversible in their configurational variations under the cyclic loading, despite the system is locally elastic everywhere. It is also found that the cyclic irreversibility is enhanced by initial film/substrate strain mismatch and material viscosity. The cyclic irreversibility is classifiable, in large, by localization irreversibility and modal irreversibility.

The localization irreversibility is found to be associated with unstable creasing of either the film surface or the substrate. The film surface creasing typically emerges as periodic setback creasing of wrinkle or fold valleys within a range of the normalized critical wavenumber, $0.40≤k¯≤1.2$⁠, while the substrate creasing induces global localization. When the setback creases globally localize later, sharp tips of the periodic film creases stabilize uniform cooperative parallel growth of the tips up to the subcritical creasing limit strain, 0.35, of the neo-Hookean half space. On the other hand, easy lateral configurational sliding of smooth fold valleys erratically prompts instability in the cooperative growth of the fold fronts, depending on the self-contact friction of the fold surfaces. Such erratic instability elicits large uncertainty in determining the starting strain of fold localization.

On PB-RPD of $k¯<0.40$⁠, the film surface does not crease, and period multiplications and folding of a wrinkle are found reversible in the corridor of $( 0.29<k¯<0.40; 0<ε<0.30)$⁠, which is considered useful for various technological applications. Modal irreversibility is observed for $0<k¯<0.29$⁠; the critical $k¯$ of irreversible DM1 and QM1 period multiplications is determined as $DIR (k¯=0.17; R=600)$ and $QIR (k¯=0.29; R=120)$ below which DM1 mode and QM1 mode become irreversible, respectively. Mode locking during unloading is found to be the main irreversibility characteristics of the PB system without film/substrate strain mismatch. In contrast, a PB of high film/substrate ratio (e.g., $R=105)$ with sufficient strain mismatch (e.g., $εps=0.4)$ and viscosity (e.g., loss tangent $=2×10−4$⁠) experiences P-period switching from DM1 and QM1 to SM2 and DM2, respectively, during unloading. It is found that material viscosity inhibits relative deformation rate of deep valleys of multimode wrinkles to that of shallow valleys, and in turn, it promotes flipping the curvature of shallow valleys during the unloading cycle. Such flipping triggers P-period switching which provides a jumping mechanism of a multimode wrinkle toward a higher energy metastable wrinkle configuration. Here, we expect that excitation of various P-period switching can control stability of cyclic multimode configuration hysteresis—cyclic divergence or shake-down. The results may inspire designing new systems with unprecedented surface modes and lead to a new way for period control of various Ruga-phases.

The authors gratefully acknowledge the supports provided by the U.S. National Science Foundation (Award No. 1462785), Institute of Molecular and Nanoscale Innovations (IMNI) at the Brown University, and Korea Institute of Science and Technology (KIST). The simulations were performed at the Center for Computation and Visualization (CCV) at the Brown University.