Multimodal Surface Instabilities in Curved Film–Substrate Structures

Structures of thin films bonded on thick substrates are abundant in biological systems and engineering applications; examples include various biological tissues [1,2], flexible electronics and wearable devices [3–6], strain-controlled superhydrophobicity [7,8], and self-assembled patterns [9]. Mismatch strains due to expansion of the films or shrinkage of the substrates can induce various modes of surface instabilities including wrinkling [10–15], creasing [16–20], period doubling [21–23], folding, ridging [24,25], delamination [26], 2D-buckling [27–30], 3D-buckling [31–33], and hierarchical wrinkles [34–36]. Some examples of surface instability patterns have been given in Fig. 1(a). Over last few decades, various modes of instabilities in the film–substrate structures have been studied individually [23,25,37,38] and systematically by categorizing them on phase diagrams [39–43]. It is well known that these instabilities can be controlled by mechanical properties such as rigidities of films and substrates, the adhesion energies between them, and their mismatch strains or compressive strains applied on the structures.

Many film–substrate structures found in nature, biological systems, and engineering applications are not flat but curved, although previous studies have been focused on flat cases. For example, Fig. 1(b) illustrates different surface instability patterns on concave film–substrate structures. From left to right (Fig. 1(b)), mucosa-growth induced creases are seen in colon; the anatomic ductus deferens in male reproductive system shows ridges on surface; muscular arteries manifest doubling patterns along the lamina; and wrinkling patterns are observed on the surface of bronchus lumen. Film–substrate structures with convex surface are commonly seen in plants, as shown in Fig. 1(c). From left to right, a bitter melon shows creases on its surface; cactus develops surface ridges; a tree stump reveals irregular doubling patterns on its bark; and a cereus presents uniform wrinkling on its skin. Other surface patterns on curved biological structures that have drawn much attention in the biochemistry and biomedical fields include morphological development of brain [44,45], growth-induced creases and folds in tubular tissues [46,47], crumpled gut surfaces [48–50], deformation of cell-culture substrates [51–54], and morphological evolution of tumor growth [55–58]. In engineering, curved film–substrate structures have been widely used in design of flexible and stretchable electronics [3,59], active detachment of biofouling [60–62], and control of aerodynamic drag [63,64].

A few theoretical and computational models have been developed to investigate the morphological instabilities and analyze the bifurcation conditions on curved film–substrate such as cylindrical [65–68] and spheroidal structures [9,69–72]. While previous studies have been focused on individual mode of instability such as wrinkle, crease, and fold, the effects of the structures’ curvature on multiple modes of surface instabilities have not been well studied or understood. In this paper, we provide a systematic study on the formation of multimodal surface instabilities on film–substrate structures with different curvatures through both theoretical analysis and numerical simulation. We will focus our study on tubular film–substrate structures with both convex and concave curvatures (Fig. 1), as these structures appear abundantly in biological systems and engineering applications. We will apply the FFT method to quantitatively categorize various instability patterns on curved structures by analyzing their wave frequencies. Thereafter, we will develop a theoretical model to predict the critical mismatch strain and wavenumber of wrinkling instability on curved film–substrate structures, followed by finite element simulation of the formation of advanced modes of instabilities. We will show that the curved film–substrate structures can delay the critical mismatch strain for wrinkling for systems with relatively large film–substrate modulus ratio, compared with flat ones with otherwise the same properties. Furthermore, concave structures promote creasing and folding and suppress ridging, but convex structures promote ridging and suppress creasing and folding. A set of phase diagrams will be calculated to guide future design and analysis of multimodal surface instabilities on curved structures.

The outline of this paper is as follows. Section 2 explains the definition of configurations and quantities in the problem. In Sec. 3, the FFT method is applied to quantitatively categorize the surface instability patterns. In Sec. 4, basic equations of continuum mechanics and incremental method in analyzing the deformation and instability of film–substrate structures are summarized. Section 5 presents the instability analysis to search for the critical mismatch strain and wavenumber for wrinkling on curved film–substrate structures. Section 6 accomplishes the finite element simulations of various modes of instabilities on film–substrate structures with various curvatures, mismatch strains, and modulus ratios. A set of phase diagrams for predicting various modes of instabilities in curved structures will be provided with the curvature’s effects discussed. Final conclusions are made in Sec. 7.

Consider the plane-strain, incompressible deformation of a tubular film–substrate structure as illustrated in Fig. 2. The film can be bonded on the inner (Fig. 2(a)) or outer (Fig. 2(d)) surface of the tubular substrate, which gives the negative or positive curvature for the film–substrate structure, correspondingly. At the reference state (Figs. 2(a) and 2(d)), there is no mismatch strain between the film and the substrate, and a material point in the structure is denoted by its cylindrical coordinates (R, Θ, Z). The structure has the film radius A and interfacial radius B at the reference state (note that for negative-curvature structure, far-end of the substrate has radius C and C ≫ B). The film thickness H at the reference state is defined as $H=|A−B|$⁠. A mismatch strain may be applied to the system by either shrinking the substrate or expanding the film. In the current study, we will shrink the substrate to introduce the mismatch strain.

where $A≤R≤B$ for concave structures and $B≤R≤A$ for convex structures. At the film–substrate interface, Eq. (4) can be reduced to $εf(R=B)=εM$⁠. As R moves from the interface $(R=B)$ to film surface $(R=A)$⁠, the compressive strain $εf$ increases in the concave structures and decreases in the convex structures.

While different modes of surface instabilities have different morphological patterns, categorizations of these surface instabilities patterns have mostly followed qualitative approaches. For example, the distinctions between wrinkles and advanced modes of instabilities including ridges, period doubles, and folds have only been qualitatively prescribed by previous researchers [15,16,24,37,39–43]. In addition, the introduction of curvatures in film–substrate systems will make the categorization of different modes of instabilities more complicated, as the curvature will affect the morphology of the same mode of instability (e.g., Figs. 3(a)–3(c)). As different modes of surface instabilities have different types of spatial periodicities, we propose to use the fast Fourier transformation method to quantitatively categorize various modes of surface instabilities by their frequency distribution. While this method has been used to analyze instability transition from wrinkle to crease of a flat graded material [37], here we further apply this method to distinguish various modes of instability patterns on curved surfaces. The frequency is in the space domain, and it is described by the wavelength of the instability pattern.

where i represents the imaginary unit, and p and j iterate from 0 to m − 1. The absolute value of y denotes the amplitude of each frequency in the current deformed surface. Here, we normalize the amplitude of the decomposed wave frequencies by the amplitude of the dominant frequency for the surface deformation at current state. As shown in Fig. 3, the deformation configuration was categorized by FFT into different states, which are patternless, wrinkle, crease, fold, double, and ridge. (Note that the film and substrate are assumed to be perfectly bonded without delamination in the current study.) By performing FFT, we find the characteristics of the deformation states are quantitatively distinguishable in the frequency domain. The detailed FFT characterization of various patterns on curved structures in Fig. 3(c) is discussed as follows.

• (i)

  Patternless state: Before the film–substrate system manifests instabilities on surface, the surface remains patternless under homogeneous deformation, and no wave frequency is detected by FFT.

• (ii)

  Wrinkling state: When the critical strain for wrinkling is reached, the initially patternless film surface develops sinusoidal wave, and a single frequency is detected by FFT.

• (iii)

  Creasing state. Unlike wrinkling, the creasing patterns are localized self-contacts developed on the initially patternless surface. The FFT spectrum of the creasing pattern consists of a wide range of wave frequencies with different amplitudes, as shown in Fig. 3(c) (crease).

• (iv)

  Folding state: After the film wrinkles during surface deformation, the valley regions of the wrinkles may further form self-contacts, giving the folding state [41]. It has a dominant wrinkling frequency coupled with many low amplitude wave frequencies, which represent the wave components close to the tip of the fold.

• (v)

  Doubling state: Doubling is another post-buckling transition from the wrinkling state [22]. During doubling evolution, initial wrinkles bifurcate into two branches, and the branching process results in two dominating frequencies for the doubled state. The large-amplitude (dominant) frequency represents the wrinkling instability, and the other frequency with smaller amplitude represents the doubling instability, as shown in Fig. 3(c) (double).

• (vi)

  Ridging state: Ridge is a height-limited regional localization that has a higher amplitude compared with wrinkle. For a curved film–substrate structure, it is usually difficult to distinguish ridge from wrinkle qualitatively. However, the wave frequency distribution of the FFT of ridge quantitatively shows that the ridging profile is not perfectly sinusoidal. There is a dominant wave frequency which corresponds to the wavelength of ridging, and this main frequency is accompanied by many small-amplitude wave components, as illustrated in Fig. 3(c) (ridge).

The surface stability of various film–substrate structures have been theoretically investigated with bifurcation and perturbation analysis, such as creasing on flat homogeneous half space [10], creasing and folding on half space with exponentially decayed modulus [37,73], wrinkling on flat film–substrate structure [13,16], instability bifurcations of thick-walled circular cylindrical [65,66,74–76], and buckling instability on cylindrical multilayer structure along circumferential and longitudinal directions [68,77]. In this paper, we will adopt the incremental theory to analyze the critical mismatch strain and critical wavenumber of wrinkling instabilities on curved film–substrate structures.

where p is the hydrostatic pressure to accomplish the incompressibility condition in the calculation.

where n denotes the normal unit vector of the surface or interface at the current state, and $‖‖$ represents the jump of a quantity.

with $ur,uθ$ the incremental displacement in the radial and circumferential directions, respectively.

where $Wi=∂W/∂λi$ and $Wij=∂2W/∂λiλj$ (no summation of index).

where $μ=μf$ when A < R < B for concave structures (B < R < A for convex structures) and $μ=μs$ when R > B for concave structures (0 < R < B for convex structures). Here, $μf$ and $μs$ are the initial shear moduli of the film and substrate, respectively.

where n is the wrinkling wavenumber on the film surface.

where $‖‖$ denotes the jump of certain quantity at the interface as defined previously.

where $λ$ and $1/λ$ are the principal stretch along the circumferential $θ$ and radial r directions, respectively, as defined in Eq. (11). We will specify $λ=λf$ for the film and $λ=λs$ for the substrate, when it is needed in later parts of the analysis.

where B is the interfacial radius at the reference state, as illustrated in Fig. 2.

The fourth-order ODE in Eq. (47) is applicable to both the film and the substrate. There are in total eight boundary conditions to be considered when solving the problem: two traction-free conditions at the film surface (Eq. (43)), two displacement conditions at the far-end of the substrate (r = c for concave structure and r = 0 for convex structure) (Eq. (44)), and four continuity conditions for traction and displacement at the interface (Eq. (45)). To solve the system, we first integrate the ODE from both the film and the substrate sides to the interface by guessing two sets of solutions which satisfy the four boundary conditions in Eqs. (43) and (44) for the film and the substrate, respectively. We further match the four continuity boundary conditions (Eq. (45)) at the film–substrate interface. The explicit matrix form of Eq. (45) is shown in the Appendix. To have nontrivial solutions, the determinant of the interface condition matrix must vanish. We can then numerically solve the eigenvalue problem to search for the critical mismatch strain $εMWrinkle$ for wrinkling of the curved film–substrate structure.

To explore the curvature effect on wrinkling, we examine critical mismatch strain $εMWrinkle$ and critical wavenumber $nc$ for the film–substrate system with different normalized curvatures $ρ=−0.1,−0.04,−0.02, 0, 0.02,0.04, and 0.1$⁠, as shown in Fig. 4(a). We plot critical mismatch strain for film–substrate system with modulus ratio $μf/μs$ ranging from 10 to 10⁴ in Fig. 4(a). The solid curves from left to right denote the critical mismatch strains for wrinkling of the film-substrate system with normalized curvature ρ = 0, 0.02, 0.04, and 0.1, respectively [15,37,41]. The dotted curves from left to right are for structures with the normalized curvature ρ = −0.02, −0.04, and −0.1, respectively. From Fig. 4(a), we can see that for different normalized curvature amplitudes, the critical mismatch strains for wrinkling deviate significantly when the modulus ratio is beyond 100. Under this condition, the curvature of film–substrate structure delays wrinkling compared with the flat film–substrate structure. In addition, Fig. 4(b) shows the critical wavenumber for the three examined normalized curvature amplitudes. Our calculation shows that the critical wavenumber for wrinkling is the same for concave and convex structures with the same curvature magnitudes.

Wrinkling is a mode of surface instability that occurs under relatively low mismatch strains and can be predicted analytically as discussed in Sec. 5. As the mismatch strain further increases, advanced modes of surface instabilities such as creasing, folding, doubling, and ridging may occur. To quantitatively understand and predict various modes of surface instabilities in curved film–substrate structures, we adopt the finite element method to systematically investigate the occurrence and evolution of multimodal surface instabilities and, in particular, the curvature’s effects on these instabilities.

We implement the two-dimensional plane-strain tubular film–substrate structures with negative and positive curvatures illustrated in Figs. 2(a) and 2(d) into finite element models using software package abaqus/standard. The 2D plane-strain four-node bilinear element with reduced integration and hourglass control (CPE4RH) is used for the negative-curvature case. For the positive-curvature case, in order to satisfy the incompressibility requirement, the 3D eight-node linear brick with reduced integration, hourglass control, and hybrid with constant pressure (C3D8RH) is used. We validate the mesh insensitivity of the models using different mesh densities, which give surface deformation with the same critical mismatch strain and instability patterns. A random perturbation with small amplitude is applied to the mesh of film surface to trigger the instability during loading. We also validate that the instability patterns are perturbation insensitive. A self-contact interaction with the type of frictionless and “hard contact” is applied to the film surface to avoid element penetration during creasing and folding deformations.

The incompressible neo-Hookean constitutive model is adopted for both the film and substrate with different shear modulus $μf$ and $μs$⁠, respectively. The modulus ratio $μf/μs$ ranges from 10⁻¹ to 10⁴ in our simulation. The mismatch strain $εM$ can be induced by prestretching the substrate, shrinking the substrate, or expanding the film. The path independency of the loading method has been discussed in literature [40]. To induce the film–substrate mismatch strain in our simulation, for the negative-curvature system, we prestretch the substrate along the radial direction, bond the film on the substrate, and then relax the prestretch in the substrate. For the positive-curvature system, we directly shrink the substrate (without prestretch) by assigning thermal expansion coefficient to the substrate and applying thermal loading to control the shrinkage. The instability patterns are evaluated and categorized by the FFT method introduced in Sec. 3. In the simulation, creasing is usually hard to be triggered on smooth surface without sufficient imperfection or undulation. Following previous studies, the creasing instability is predicted to occur, when the compressive strain on the film surface reaches a critical value, 0.35 [17,46].

We summarize the results from the finite element simulations in phase diagrams given in Figs. 5(a)–5(e), for five normalized curvatures  $ρ=−0.1,−0.04, 0,0.04,and 0.1$⁠, respectively. Each phase on the phase diagrams represents a type of surface instability pattern as categorized by the FFT method introduced in Sec. 3, and the governing parameters for phase diagrams are normalized curvature, modulus ratio, and mismatch strain of the film–substrate system. From Fig. 5, it can be seen that the phase diagrams of surface instability patterns share a number of common features. Patternless and wrinkling phases occur at relatively small mismatch strains. Compared with patternless and wrinkling, other modes of instabilities take place when the mismatch strains are relatively large. Creasing and folding states develop in film–substrate structures with relatively low modulus ratios, while doubling and ridging states require higher modulus ratios. Despite these common features, the curvatures of film–substrate structures can affect the instability patterns significantly (Fig. 5). To understand the film–substrate curvature’s effects, we compare the surface instability patterns in film–substrate structures with different curvatures but the same mismatch strains and modulus ratios, as demonstrated in Figs. 5(a)–5(e).

Evidently, the critical mismatch strain for creasing increases with the normalized curvature of the film–substrate structure. For curved film–substrate structure with $ρ<0$⁠, the critical mismatch strain $εMCrease$ for creasing is smaller than 0.35. On the contrary, $εMCrease$ is larger than 0.35 for system with $ρ>0$⁠.

Folding is promoted by negative curvature and suppressed by positive curvature. The curvature of the film–substrate structure has significant effects on the phase boundaries between folding and adjacent phases including wrinkling, doubling, and ridging. From Figs. 5(a)–5(e), it reveals that the area of the folding region (marked by yellow) on the phase diagram decreases significantly when the surface of the film–substrate structure goes from concave to flat, then to convex (⁠$ρ$ goes from negative to zero, then to positive), which means that folding is promoted by concave surface. Folding can take place right after wrinkling instability when the modulus ratio is relatively low (near $μf/μs=1.3$⁠). In this case, the local film strain at the wrinkling valley reaches the critical mismatch strain for creasing (⁠$εf=0.35$⁠) immediately after wrinkling, and the surface initiates a crease tip at the valley of the wrinkle. Folding can also manifest upon well-developed wrinkles, when $μf/μs$ is larger. The phase diagram in Fig. 5(a) for concave structure $ρ=−0.1$ shows that the critical mismatch strain to trigger folding is within a range of $0.29≤εM≤0.40$⁠, for modulus ratio $1.3≤μf/μs≤80$⁠. For concave surface with normalized curvature $ρ=−0.04$ as shown in Fig. 5(b), folding initiates on film surface under mismatch strain $0.33≤εM≤0.35$⁠, for modulus ratio $1.3≤μf/μs≤15$⁠. Figure 5(c) illustrates the phase transitions for flat film–substrate structure (⁠$ρ=0$⁠). From the phase diagram, we can see that the folding transition on film surface starts with mismatch strain $0.35≤εM≤0.45$⁠, for modulus ratio $1.3≤μf/μs≤12$⁠. For convex film–substrate structures, the folding region is even more significantly reduced, as shown in the phase diagrams in Figs. 5(d) and 5(e). The critical mismatch strain for folding is postponed to $0.37≤εM≤0.40$ and $0.41≤εM≤0.45$⁠, for normalized film–substrate curvatures $ρ=0.04 and 0.1$⁠, respectively. Moreover, the range of modulus ratio for folding is shrunken to $1.3≤μf/μs≤6$ for both cases.

In general, negative curvature promotes folding by requiring relatively small mismatch strains for patternless-folding transition, and allowing relatively wide ranges of film–substrate modulus ratio for the folding phase. This trend can be understood as follows. The inhomogeneous distribution of the circumferential compressive strain throughout the thickness of the film leads to the curvature effect on surface instabilities. With the same interfacial mismatch strain $εM$ between the film and substrate, concave structure gives higher compressive strain on the film surface than flat structure, and thus promotes folding. On the contrary, convex structure gives lower compressive strain on the film surface compared to the flat structure, and thus suppresses folding.

We further use a few points with the same modulus ratio and mismatch strain (i.e., $μf/μs=30,εM=0.5$ for points (A–E) and $μf/μs=30,εM=0.6$ for points (F–J)) on the phase diagrams with different normalized curvature (i.e., $ρ=−0.1,−0.04,0, 0.04, and 0.1$⁠) to discuss curvature’s effect on folding, as shown in Figs. 6(a) and 6(b). Points (A) and (F) show folding instability on the film surface, for normalized film–substrate curvature $ρ=−0.1$⁠. Folding disappears on film surface when the normalized film–substrate curvature is increased to $ρ=−0.04$ or larger. In these cases, other instability phases such as doubling and ridging are observed and will be discussed in the following parts.

Ridging is promoted by positive curvature and suppressed by negative curvature. On flat film–substrate structure, ridging instability is observed on film surface as an outward localization for systems with relatively high mismatch strains and modulus ratio. As discussed in Sec. 3, ridging is characterized by FFT and its wave frequency components reveal the nonaxisymmetric geometry compared with wrinkling. Its FFT spectrum shows a dominant wave frequency, which is accompanied by massive small-amplitude frequencies as shown in Fig. 3(c). Here, our simulation results further show that for curved film–substrate structures, the curvature plays an important role in determining the ridging phase.

The phase diagram in Fig. 5 indicates that the area of ridging zone increases with the normalized film–substrate curvature. Flat film–substrate structure (⁠$ρ=0$⁠) shows that the range of the critical mismatch strain for wrinkling–ridging transition is within $0.30≤εM≤0.40$⁠, for modulus ratio $μf/μs≥150$⁠. For convex system with film–substrate curvature $ρ=0.04$ shown in Fig. 5(d), the critical mismatch strain for ridging is in the region of $0.35≤εM≤0.45$⁠, for a wide range of modulus ratio $μf/μs≥10$⁠. By further increasing the normalized film–substrate curvature to $ρ=0.1$⁠, ridging can be triggered at a small mismatch strain around $εM=0.30$ for very small modulus ratio $μf/μs≥1.3$ (illustrated in Fig. 5(e)). Ridge-folds or ridge-doubles will develop on film surface under continuous loading. On the contrary, for concave film–substrate structure $(ρ<0)$⁠, ridging is significantly suppressed. As illustrated in Fig. 5(b) for system with normalized curvature $ρ=−0.04$⁠, the triggering strain for ridging is postponed to $εM=0.50$ compared with flat and convex film–substrate structure, and the corresponding range of modulus ratio for ridging is $μf/μs≥60$⁠. If we further increase the normalized curvature of the system to $ρ=−0.1$⁠, no ridging transition is observed on film surface, as shown in Fig. 5(a).

The contour plots for selected cases (A–E), (F–J), and (K–O) in Figs. 6(a)–6(c) illustrate how positive curvature (convex structure) promotes ridging and suppresses folding. For system with mismatch strains and modulus ratios $εM=0.5,μf/μs=30$⁠, points (D) and (E) in column (a) show ridging instability for normalized film–substrate curvature $ρ=0.04 and 0.1$⁠, respectively. In column (b) with $εM=0.6,μf/μs=30$⁠, ridge is observed for point (I) and ridge-double is seen for point (J). For system with $εM=0.6,μf/μs=300$ as shown in column (c), points (L) and (N) reveal ridging instability and point (O) develops ridge-double on film surface.

In this paper, we studied the curvature’s effects on multimodal surface instabilities in tubular film–substrate structures, which are commonly found in nature and engineering applications. We proved that the curvature along with the mismatch strain and modulus ratio of the film–substrate system determines the occurrence of various instabilities including wrinkling, creasing, folding, doubling, and ridging. We applied the FFT method to characterize the surface instabilities quantitatively, by decomposing the wave frequencies of each instability mode. We then provided a systematic study on the formation of multimodal surface instabilities on curved film–substrate structures with different curvatures through combined theoretical analysis and numerical simulation. We used the incremental theory to predict the critical mismatch strain for wrinkling and wavenumber for systems with different normalized film–substrate curvatures, under a wide range of film–substrate modulus ratio (⁠$1<μf/μs≤104$⁠). The results showed that large curvature magnitudes can delay wrinkling significantly when the modulus ratio is larger than 100. Our finite element simulations on post-buckling of the curved film–substrate system showed that curvature plays an important role in determining the phase boundaries for various modes of instabilities. Concave structure (negative curvature) has the trend of promoting folding and suppressing ridging. Convex structure (positive curvature) promotes ridging and suppresses folding. A set of phase diagrams are calculated to systematically understand the curvature’s effect on multimodal instabilities in film–substrate structures.

This work was supported by ONR (No. N00014-14-1-0528), MIT Institute for Soldier Nanotechnologies, and NSF (No. CMMI-1253495).

The eigenvalue problem in Eqs. (A6)–(A9) can then be solved numerically to search for the critical film–substrate mismatch strain for wrinkling.