Direct numerical simulation of buckling instability of thin films on a compliant substrate

Abstract

Surface buckling (wrinkling) driven by mechanical instability is commonly observed in thin-film structures with a compliant substrate. The resulting undulation, while sometimes undesirable, has been increasingly exploited to enhance mechanical and/or functional performances of many thin film devices. In this study a practical finite element modeling approach is introduced to simulate wrinkle formation in thin films atop a compliant substrate. The proposed technique is robust and easy to implement, and it overcomes typical challenges in computationally modeling the buckling instability. Using a two-dimensional geometry under the plane strain or generalized plane strain conditions, with randomly distributed imperfections bearing different material properties at the film/substrate interface, we demonstrate the model’s capability in triggering surface instability during direct compression and out-of-plane tensile loading. With sufficient mesh refinement, the predicted wrinkling wavelength, amplitude, and critical strain to activate wrinkle formation are shown to be close to analytical solutions. The effect of imperfection distribution is systematically studied, and a valid range of imperfection spacing is identified. The present numerical approach can be applied to predicting buckling instability in the design and analysis of thin film/compliant substrate systems over a wide range of material and geometric conditions. Directions for future studies are also discussed.

Keywords

Buckling wrinkling thin film finite element analysis instability

Introduction

Formation of buckles (wrinkles) on a thin film attached to a compliant substrate has received considerable attention due to the advancement of flexible electronics,^1–5 optoelectronics and photovoltaics,^6–11 and various functional devices,^12–16 as well as the development of self-organizing mechanisms.^17–21 Buckling can be triggered by direct application of an in-plane compressive strain, as schematically shown in Figure 1(a). It may also occur in response to applied tension, as shown in Figure 1(b), where a mismatch in lateral contraction (caused by the Poisson’s ratio effect) between the thin film and thick substrate forces the film into compression. Wrinkling relaxes the compressive strain in the coating, therefore reducing the elastic strain energy. Assuming elastic deformation with a semi-infinite substrate, the critical strain for buckling instability and the wavelength of the wrinkling pattern can be derived from linear buckling theories^22–25 (key equations to be presented in section “Brief overview of theory”). Thin-film metrology has also been developed, where the elastic modulus of a film, if initially unknown, can be calculated after measuring the wrinkling wavelength.^25–27

Figure 1.

Schematics showing the formation of buckles (wrinkles) during (a) compression along the x-direction and (b) applied tension along the z-direction. The scenarios in (a) and (b) are also termed direct compression and indirect compression, respectively.

While the analytical formulation has been widely employed in analyzing experimental data, relatively few studies have attempted to use numerical simulation to address wrinkle formation.^28–34 There are inherent challenges associated with instabilities in finite element modeling. To overcome the existing numerical complexities, in some studies the surface wrinkling phenomenon was simulated in two steps of pre-buckling and post-buckling analysis.^31,32,34 The linear modal analysis was mainly performed for the pre-buckling step, and the post-buckling analysis was done via the combination of nonlinear material/geometry models and the imperfection techniques. For the post-buckling step, the incompressible neo-Hookean and Arruda–Boyce constitutive models for elastomers were exploited that require the use of incompressible/hybrid elements. Although applying such elements may have some benefits within the reported numerical frameworks as to simulate the higher-order deformation modes at large strains, incompressible/hybrid elements are available only in limited finite element packages such as ABAQUS.³⁵ These types of elements are based on relatively complex mathematical theories, and are computationally more expensive compared with regular elements, thus significantly lowering the efficiency for three-dimensional analysis. In addition, they are applicable only for the static analysis and can lead to convergence issues when used with material models that exhibit volumetric plasticity.³³ Aside from the type of elements, other special treatments are required to capture the bifurcation and post-bifurcation wrinkling shapes.

The imperfection methods were categorized into two general types of material and geometrical imperfections, and were widely used in the simulation of instabilities in plates, shells, and cylinders. The geometrical imperfection method, employed more frequently for studying surface instability, may include the mesh, geometry, and boundary condition perturbation techniques (some were directly implemented in the commercial software packages such as ABAQUS). In some studies, a sinusoidal perturbation was directly imposed on the geometry of the film layer which dictates the formation/propagation of sinusoidal-form waves.^32,33 In another study,²⁸ two different types of sinusoidal imperfection and a periodic array of non-interacting exponential surface depressions were specified within the model and studied separately. In a similar fashion, a mesh perturbation technique was applied in that the mesh nodal points were displaced by the sinusoidal perturbation displacement field function.³¹ All such imperfection approaches are relatively complex, and they essentially build the formation/propagation of sinusoidal waves into the model. The implementation can be laborious and requires calibration by setting many free parameters. Furthermore, the interpretation and verification of results are not straightforward, which makes it less practical for common users. A more straightforward approach for modeling the wrinkle formation was reported,³⁶ which however involved special geometric features built into the film surface. Aside from the aforementioned geometrical approach, the material imperfection method is typically based on perturbation via the damage and plasticity techniques. For instance, a non-uniform plastic behavior was built into the model containing an elastoplastic film.³⁷ Some related studies regarding the material imperfection approach can be found in Andrews and Massabò.³⁸

The present work aims at developing a simple computational modeling approach to simulate temporal evolution of buckling. Having this capability enables increased accuracy and versatility in predicting wrinkle formation of a relatively stiff thin film on a thick and compliant substrate. The numerical approach may also be employed as a design tool and to study the effect of finite substrate thickness over a wide range of elastic property relationship between the constituents. Moreover, the effects of multilayer thin films and possible viscoelastic and plastic behavior, as well as delamination damage, may all be explored. Our current work seeks to establish a basic modeling framework and assess its feasibility. The proposed methodology is direct, in that the pre/post-buckling simulations can be performed successfully in only one analysis step. The approach is also robust and relatively easy to perform, and it can be implemented with any common finite element code and analysis platform.

In this study, two different film-substrate systems are considered. The first model system consists of a conjugated polymer blend of regioregular poly (3-hexylthiophene) and phenyl-C61-butyric acid methyl ester (P3HT:PCBM) as the thin film, and the intrinsically conductive polymer poly(3,4-ethylenedioxythiophene) and polystyrene sulfonate acid (PEDOT:PSS) is the thin film for second model. The elastomeric polydimethylsiloxane (PDMS) is considered as the substrate for both systems. Note that this is a part of the mechano-optoelectronic composite structure being developed as a functional component of the autonomous structural composites capable of self-powered strain sensing.^39–41 The P3HT:PCBM or PEDOT:PSS film was prepared on pre-tension strained PDMS. Upon release of the pre-strain, the film experiences compression leading to buckling. While corrugated films showed improved mechanical resiliency, their optoelectronic performance may be affected.^39,41 Developing a robust numerical methodology to predict wrinkling is thus essential in designing the composite layout along with its fabrication.

In the remainder of this article, a brief overview of the buckling theory is first given in section “Brief overview of theory.” Section “Numerical model description” provides the model description. Simulation results are presented in section “Results and discussion,” with attention devoted to comparing the model predictions with the theory on buckling wavelength, critical strain, and amplitude. Details on the analysis procedure and the requirements for the placement of imperfections are also discussed in this section. Highlights of the numerical approach and salient findings are summarized in section “Conclusion.” As evidenced from the current Introduction section, the terms buckling and wrinkling are used interchangeably throughout this article.

Brief overview of theory

Various theoretical analyses on surface buckling instability of a thin film on a compliant substrate have been reported. Although there are no full-field closed-form analytical elastic solutions, some approximated planar solutions are available in the literature.^22–25 Here, the force balance approach^22,25 is followed and briefly discussed, and the key parameters related to surface buckling are introduced.

The problem considers a semi-infinite substrate under plane strain deformation, with negligible shear stress at the film-substrate interface. It is assumed that the buckling waves are sinusoidally formed throughout the entire width of the film. The classical ordinary differential equation for bending of an elastic film on an elastic substrate subjected to direct compression is

(\frac{E_{f} I}{1 - ν_{f}^{2}}) \frac{d^{4} y}{d x^{4}} + F \frac{d^{2} y}{d x^{2}} + (\frac{E_{s} π L}{(1 - ν_{s}^{2}) λ}) y = 0

(1)

where F is the uniaxially applied compressive force in the film, $E_{f}$ and $ν_{f}$ are, respectively, Young’s modulus and Poisson’s ratio of the film, $E_{s}$ and $ν_{s}$ are, respectively, Young’s modulus and Poisson’s ratio of the substrate, L is the out-of-plane length (here in z direction), $λ$ is the wavelength, $I = L {t_{f}}^{3} / 12$ is the moment of inertia of the film, and $t_{f}$ is the thickness of the film. The solution of equation (1) in a sinusoidal form is expressed as

y = A \sin \frac{2 π x}{λ}

(2)

where A is the amplitude of the waves. Without going through the details, the compressive force in the film, $F = F (λ)$ , can be obtained by substituting equation (2) into equation (1) and solving for F. Subsequently, by minimizing F with respect to $λ$ , the critical buckling force in the film, $F_{cr}$ , and the buckling wavelength can be obtained. Therefore, the film buckling wavelength is

λ = 2 π t_{f} {(\frac{(1 - ν_{s}^{2}) E_{f}}{3 (1 - ν_{f}^{2}) E_{s}})}^{\frac{1}{3}}

(3)

Note that the wavelength shown in equation (3) is only a function of the thickness of the film and the elastic properties of the film and substrate. The critical buckling stress and strain can also be determined, using $F_{cr}$ , as

σ_{cr} = {(σ_{xx})}_{cr} = \frac{F_{cr}}{t_{f} L} = {(\frac{9 E_{f} E_{s}^{2}}{64 (1 - ν_{f}^{2}) {(1 - ν_{s}^{2})}^{2}})}^{\frac{1}{3}}

(4)

and

e_{cr} = {(e_{xx})}_{cr} = \frac{σ_{cr}}{(E_{f} / (1 - ν_{f}^{2}))} = \frac{1}{4} {(\frac{3 (1 - ν_{f}^{2}) E_{s}}{(1 - ν_{s}^{2}) E_{f}})}^{\frac{2}{3}}

(5)

Also, the amplitude of the waves, A is derived with the assumption that once the applied compressive strain, e, exceeds $e_{cr}$ , the wavelength remains constant ( $λ$ is independent of e). It is also assumed that $A << λ$ . Hence, the wave’s amplitude is derived via the tensile strain release approach²⁵ and expressed as

A = (\frac{λ}{π}) {(e - e_{cr})}^{\frac{1}{2}} = t_{f} {(\frac{e}{e_{cr}} - 1)}^{\frac{1}{2}}

(6)

The parameters $λ$ , $e_{cr}$ , and A are the key variables of thin film buckling instability. In addition to direct compression along the x direction, this study also considers applied tension in the z direction. If $ν_{s} > ν_{f}$ , the substrate will undergo greater lateral contraction (along x, see Figure 1(b)) than the film, thus forcing the film into compression. Under this circumstance, the critical compressive stress for triggering film buckling is caused by the difference in lateral contraction between substrate and film, namely, $(ν_{s} - ν_{f}) e_{zz}$ . As a consequence, the critical tensile strain $(e_{zz})_{cr}$ is related to the critical compressive strain $(e_{xx})_{cr}$ by

{(e_{zz})}_{cr} = \frac{{(e_{xx})}_{cr}}{ν_{s} - ν_{f}}

(7)

Note that since the critical strain $(e_{xx})_{cr}$ in equation (5) was derived based on the plane strain condition, the resulting expression of $(e_{zz})_{cr}$ in equation (7), with the applied tensile strain in z, should be viewed as an approximation. In section “Results and discussion,”equations (3), (5)–(7) are compared with the numerical results.

Numerical model description

The finite element simulations are performed on the two-dimensional film-substrate system as shown in Figure 2. Buckling instability under both cases of direct and indirect compression are studied. For the case of direct compression, the surface wrinkles are simulated via both generalized plane strain and plane strain models; however, for the case of indirect compression, plane strain is fundamentally invalid so only the generalized plane strain condition is considered. The generalized plane strain condition takes into account the uniform deformation in the out-of-plane (z) direction and thus the realistic three-dimensional effect. Note that the plane strain assumption is consistent with the theoretical formulation, equations (3)–(6). The problem domain is defined with equal width and depth, $w = d = 1000 μ m$ , film thickness $t_{f} = 0.1 μ m$ , and the out-of-plane thickness used in generalized plane strain, L, of unity. Note that once wrinkles have developed, their wavelength and amplitude scale linearly with the film thickness according to the theories (equations (3) and (6)). In the current numerical study, we consider a single film thickness to study the feasibility of our numerical approach and to compare with theoretical values. Both the thin film and substrate materials are assumed to be linear-elastic. For the first system, the film is P3HT:PCBM with the Young’s modulus of $E_{f} = 7300 MPa$ and Poisson’s ratio of $ν_{f} = 0.35$ ; for the second system, the film is PEDOT:PSS with the Young’s modulus of $E_{f} = 2000 MPa$ and Poisson’s ratio of $ν_{f} = 0.35$ .⁴² For the compliant PDMS substrate in both cases, $E_{s} = 2.97 MPa$ and $ν_{s} = 0.495$ .⁴³ Note that the Poisson’s ratio of PDMS is defined slightly smaller than 0.5 to avoid potential convergence problems. The simulation domain is treated as a representative unit segment of a large periodic structure along x. The left-hand edge is a symmetry boundary, with displacement allowed only in the y direction while the bottom corner point is fixed in both x and y. The right-hand boundary can move but is constrained to remain vertical during deformation. For direct compression, the negative x-direction displacement is prescribed on the right-hand boundary. For indirect compression, a pulling action is imposed along the z direction resulting in lateral contraction.⁴⁴

Figure 2.

The film-substrate system and boundary conditions used for numerical simulations.

The model generation and simulations were all carried out using the commercial finite element software package ABAQUS (Version 2017, Dassault Systems Simulia Corp., Johnston, RI, USA). Four-noded quadrilateral continuum elements are utilized throughout the model. The top layer is assigned to be the thin film, with the element size increasing gradually from top to bottom. To help initiate buckling instability in the present static linear-elastic finite element model, an “embedded imperfections” approach is adopted here. Some elements were randomly chosen to be endowed with perturbed material properties. In this study, the imperfections, as schematically shown in Figure 2, are distributed in the substrate immediately below the film/substrate interface, with their elastic properties following those of the film material. This proposed methodology is somewhat similar to the material imperfection approach mentioned previously. However, one may interpret this practice as having an imperfect interface with occasional geometric irregularities, which is also physically plausible. Therefore, the imperfections here may be viewed as bearing both the material and geometric natures, but the methodology is straightforward to implement and a fully elastic system may be simulated without any other artificial treatment. It is illustrated in the following section that the methodology is able to generate well-defined wrinkles. The details regarding the desired number and distribution of the embedded imperfections, and the associated wrinkle configurations, are discussed in the following section.

Once a wavy pattern forms, the wavelength $(λ)$ and amplitude (A) can be directly obtained from the simulation output. The critical buckling strain $(e_{cr})$ may be identified by monitoring the evolutions of overall load and displacement. A deviation from the pre-buckling path signifies the onset of instability. The details are presented in the following section.

Results and discussion

In this section, the simulation results are presented, with emphasis given to the wrinkling wavelength, amplitude, and critical strain. For illustration purposes, a mesh consisting of 500 columns and 100 rows of elements (total of 50,000) is initially considered for both the direct and indirect compression problems (sections “Direct compression” and “Indirect compression”). The effect of imperfection placement is discussed in section “Effect of imperfection distribution.” Mesh convergence analyses are presented in section “Convergence analysis and comparisons with theory,” along with the comparisons between simulation results and analytical solutions described in section “Brief overview of theory.” Section “Further discussion” provides additional discussion of the numerical technique.

Direct compression

Simulations based on both the plane strain and the more realistic generalized plane strain conditions were performed, but this section only includes the results of generalized plane strain (for P3HT:PCBM film). The plane strain results are presented in section “Convergence analysis and comparisons with theory,” in the context of making comparisons with theories. The wavelengths, critical strains, and amplitudes obtained from all simulations for the two film-substrate systems are also listed in Tables 1 and 2.

Table 1.

Convergence study results (P3HT:PCBM film).

Number of elements	Direct compression						Indirect compression
	Wavelength $(μ m)$		Amplitude $(μ m)$		Critical buckling strain		Wavelength $(μ m)$	Amplitude $(μ m)$	Critical buckling strain
	GPE	PE	GPE	PE	GPE	PE	GPE
50,000	28.40	28.50	0.305	0.286	0.0184	0.0156	28.00	0.271	0.1440
100,000	18.32	18.80	0.155	0.142	0.0120	0.0098	17.72	0.158	0.0893
200,000	12.64	11.90	0.097	0.090	0.0077	0.0065	12.16	0.099	0.0564
400,000	6.96	6.97	0.043	0.046	0.0051	0.0043	6.88	0.044	0.0365
800,000	5.23	4.88	0.029	0.029	0.0036	0.0030	4.85	0.028	0.0252
1,600,000	4.30	4.30	0.014	0.012	0.0029	0.0024	4.20	0.016	0.0207
Theory		5.59		0.020		0.0032	5.59	0.020	0.0218

P3HT: poly (3-hexylthiophene); PCBM: phenyl-C61-butyric acid methyl ester; GPE: generalized plane strain; PE: plane strain.

Table 2.

Convergence study results (PEDOT:PSS film).

Number of elements	Direct compression						Indirect compression
	Wavelength $(μ m)$		Amplitude $(μ m)$		Critical buckling strain		Wavelength $(μ m)$	Amplitude $(μ m)$	Critical buckling strain
	GPE	PE	GPE	PE	GPE	PE	GPE
50,000	17.20	17.32	0.320	0.330	0.0430	0.0360	15.60	0.300	0.3150
100,000	11.67	11.72	0.140	0.150	0.0280	0.0230	10.89	0.130	0.1870
200,000	7.84	7.86	0.100	0.090	0.0180	0.0150	7.51	0.090	0.1170
400,000	5.17	5.18	0.060	0.060	0.0120	0.0100	5.03	0.050	0.0770
800,000	3.47	3.47	0.030	0.027	0.0090	0.0070	3.41	0.028	0.0550
1,600,000	2.79	2.79	0.016	0.014	0.0070	0.0060	2.75	0.014	0.0430
Theory		3.63		0.020		0.0075	3.63	0.020	0.0516

PEDOT: poly(3,4-ethylenedioxythiophene); PSS: polystyrene sulfonate acid; GPE: generalized plane strain; PE: plane strain.

A linearly ramped negative x-displacement is applied to the right-hand boundary in Figure 2 to simulate direct compression. Both the cases with and without the embedded imperfections are considered. Figure 3 shows the deformed configurations for the generalized plane strain model when the applied displacement reaches 20 µm (or equivalently, 0.02 compressive engineering strain), for the cases with (Figure 3(a)) and without (Figure 3(b)) the embedded imperfections. As can be seen in Figure 3(b), no buckling occurred in the perfect model. With the imperfections, Figure 3(a), however, a uniform wave form was generated throughout the entire width of the specimen in spite of the relatively coarse mesh used. The simulated wavelength and amplitude for this mesh size (50,000 elements) are listed in Table 1.

Figure 3.

Deformed configurations for the direct compression problem (under the generalized plane strain condition) when the applied compressive displacement reaches 20 µm, for the case (a) with embedded imperfections and (b) without imperfection. A displacement scaling factor of 10 is used here for better visualization.

The onset of buckling can be identified in several ways, by comparing the simulation history results with and without imperfections. Figure 4(a) shows how the induced vertical (y) displacement at the upper-right corner node evolves with the applied overall compressive x-displacement at the right-hand boundary. The zoomed-in plot is used to identify where the curve for the case with imperfections starts to deviate from that of the perfect model. Figure 4(b) shows the overall reaction force as a function of the applied compressive displacement. Again, a clear deviation point can be identified. In both cases, the critical applied displacement is found to be −18.4 µm, which is equivalent to a critical strain e_cr of −0.0184. By monitoring the deformed shape of the model, it was confirmed that the onset of buckling corresponds to this critical strain.

Figure 4.

(a) Variation of the induced vertical displacement at the upper-right corner node with the applied displacement during direct compression and (b) Variation of the reaction force with the applied displacement during direct compression.

Indirect compression

A linearly ramped tensile displacement in the z direction is applied to the generalized plane strain model. Figure 5 shows the deformed configurations when the applied tensile strain e_zz reaches 0.16, for the cases with (Figure 5(a)) and without (Figure 5(b)) the embedded imperfections. No buckling occurred in the perfect model as seen in Figure 5(b). With imperfections, however, lateral contraction due to the Poisson’s ratio effect instigated buckling and a uniform wave form was generated throughout the entire width of the specimen (Figure 5(a)) despite the relatively coarse mesh used. The simulated buckling wavelength and amplitude for this mesh size (50,000 elements) are also included in Table. 1.

Figure 5.

Deformed configurations for the indirect compression problem when the applied tensile strain e_zz reaches 0.16, for the case (a) with embedded imperfections and (b) without imperfection. A scaling factor of 5 is used here for better visualization.

A similar procedure as outlined in section “Direct compression” was used to obtain the critical strain. Figure 6(a) and (b) shows the evolutions of the vertical displacement u_y and horizontal displacement u_x, respectively, at the upper-right corner node as the applied tensile strain e_zz increases. In this case of indirect compression, the deviation of the curves with imperfections from their perfect counterparts can be easily detected in only the displacement response as shown in Figure 6. No noticeable deviation can be seen from the load-displacement curve. A consistent critical tensile strain value $(e_{zz})_{cr}$ of 0.144 is obtained from both Figure 6(a) and (b). Again, it was confirmed from the deformed shape that the onset of buckling corresponds to this critical strain.

Figure 6.

Variation of the induced: (a) vertical displacement and (b) horizontal displacement at the upper-right corner node with the applied displacement u_z during indirect compression (tension in the z direction).

Effect of imperfection distribution

The embedded imperfections used in the current approach are elements immediately below the interface, with the elastic properties identical to those of the film. A large number of simulations were conducted to explore how the distribution of imperfections can lead to uniform and reliable wrinkling characteristic of a given film-substrate system. It was found that assigning only one imperfect element in the model is sufficient to trigger instabilities; however, different distributions of such defects in the model may alter the induced wrinkling pattern. Here, a systematic study on the effects of uniform and non-uniform distributions of imperfections on the wrinkling patterns is presented.

Consider the film-substrate system shown in Figure 2 (with the dimensions and boundary conditions defined in section “Numerical model description”). Starting with the case of uniformly distributed imperfections, models with various numbers of imperfect elements (or equivalently, with different values of imperfection spacing, $S_{imp}$ ) are analyzed for the two film materials. In each model, a mesh with the total number of 1,600,000 elements is considered (the finest mesh size yielding converged results in this study, section “Convergence analysis and comparisons with theory”). The simulated wrinkle wavelength versus the number of imperfections are plotted in Figure 7(a), and the wavelength versus the imperfection spacing, $S_{imp}$ , is shown in Figure 7(b). It should be noted that all models in this section were generated with the same number, size, and configuration of elements. As can be seen from Figure 7, the obtained wavelength is sensitive to the number of embedded imperfections; however, for a wide range of imperfection spacing, the wavelength tends to remain constant. This stable value of wavelength is henceforth denoted as $λ_{s}$ .

Figure 7.

Variation of the wrinkle wavelength, λ, with (a) the total number of uniformly embedded imperfections and (b) the uniform imperfection spacing, S_imp. The results pertain to the generalized plane strain (GPE) condition, at e/e_cr = 1.20.

The results presented in Figure 7 for the two film materials can be combined by normalizing both the wavelength and imperfection spacing with $λ_{s}$ (which is different for the two film materials), as shown in Figure 8. It is evident that the range of imperfection spacing which generates the constant wavelength coincides for the two film materials, suggesting the generality of the current approach. The curves in Figure 8 can be divided into three regions as follows. Region (II) $(6 ⩽ S_{imp} / λ_{s} ⩽ 55)$ is associated with the constant/stable wavelength $(λ / λ_{s} = 1)$ , in which the wrinkling wavelength is invariant with respect to the imperfection distribution. The wrinkling configuration within region (II) is sinusoidal and uniform throughout the film, a typical case of which is shown in Figure 9(b). In the other two regions defined in Figure 8, regions (I) and (III), the value of $λ / λ_{s}$ deviates from unity. The simulated wave configurations differ from the conventional sinusoidal form. Note that in region (III), the imperfection density is very low $(S_{imp} / λ_{s} > 55)$ . The resulting wrinkle pattern is not uniform; the instability was initiated locally at the imperfection site and the wrinkles do not propagate uniformly throughout the entire width (a typical form is shown in Figure 9(a)). Finally, region (I) in Figure 8 corresponds to very high imperfection densities $(S_{imp} / λ_{s} < 6)$ ; in this case, depending on the value of $S_{imp} / λ_{s}$ , various forms of wavy patterns were seen which were dictated by the imperfection distribution and their interactions. Within the span of $1 < S_{imp} / λ_{s} < 4$ , uneven wave patterns were observed (representative forms are shown in Figure 9(c) and (d)). For the case of $S_{imp} / λ_{s} ⩽ 1$ , a sinusoidal wave with the wavelength artificially controlled by the imperfection spacing was seen (a typical form is shown in Figure 9(e)). Note that for the cases of non-uniform wrinkle configurations $(λ / λ_{s} \neq 1)$ , the maximum measurable wavelength values are reported here.

Figure 8.

The curves in Figure 7(b) normalized by the value of stable wavelength. Three regions with different ranges of imperfection spacing are identified: (a) Imperfection controlled waves (S_imp/λ_s < 6), (b) Fully developed uniform waves independent of the imperfections (6 ⩽S_imp/λ_s⩽ 55) and (c) Non-uniform/localized waves (S_imp/λ_s > 55).

Figure 9.

Typical wrinkling configurations obtained from the simulations for various uniform imperfection distributions, for the models with PEDOT:PSS film (the green dots are the embedded imperfections barely visible in these images). (a) Non-uniform/localized wave pattern (S_imp/λ_s > 55) (b) Fully sinusoidal wave form independent of the imperfection distribution (6 ⩽S_imp/λ_s⩽ 55) (c) and (d) Complex wave patterns as the imperfection spacing becoming smaller (S_imp/λ_s = 3.60 in (c) and S_imp/λ_s = 2.40 in (d)) (e) Fully imperfection controlled sinusoidal wave (S_imp/λ_s = 0.89).

In addition to uniform imperfection distributions, a large number of simulations were performed involving random (non-uniform) imperfection distributions. Here, only the salient findings are reported. Depending on the maximum and minimum values of the imperfection spacing, $S_{imp}$ , in a given model, the buckling pattern for the random imperfection distribution can be both uniform or non-uniform. The sinusoidal and fully uniform waveform, independent of imperfection locations, can only be generated when both the conditions $6 ⩽ (S_{imp})_{\max} / λ_{s} ⩽ 55$ and $(6 ⩽ (S_{imp})_{\min} / λ_{s} ⩽ 55)$ are met. These conditions correlate well with the satisfactory range identified in the case of uniform imperfection distributions, that is, region (II) as shown in Figure 8. This consistency further verifies the generality of the valid range of $S_{imp} / λ_{s}$ in generating uniform surface wrinkles. Deviation from the conditions leads to various non-uniform and complex wave patterns. Figure 10 shows representative wrinkle patterns induced by random imperfection distributions. Note that in Figure 10(a), a non-uniform wavy configuration forms locally, but in regions far away from densely populated imperfections, a common sinusoidal pattern may still develop.

Figure 10.

Typical wrinkling patterns seen for nonuniform distribution of imperfections (for P3HT:PCBM film). (a) Nonuniform and complex wrinkles for a case that does not satisfy the condition of 6 ⩽S_imp/λ_s⩽ 55. (b) Uniform wrinkles obtained for a case that satisfies the condition of 6 ⩽S_imp/λ_s⩽ 55, independent of the imperfection locations. The images cover only a portion of the surface region. The barely visible green dots are the embedded imperfections, and the scaling factor is 40 for both (a) and (b).

In summary, in the current numerical approach the embedded imperfections are required to trigger buckling instabilities, but a range of $S_{imp} / λ_{s}$ was identified such that the wrinkle pattern is not controlled by the number and locations of the imperfections. For any arbitrary distribution of imperfections, an “objective” wrinkle pattern can be achieved provided that the extreme values of imperfection spacing in a given model fall within the range of $6 ⩽ S_{imp} / λ_{s} ⩽ 55$ . The presented upper and lower bounds for $S_{imp} / λ_{s}$ are specified conservatively based on the numerical results. It should be noted that, in addition to the wavelength results presented here, the critical buckling strain and amplitude were also found to remain invariant for any random imperfection distribution within this valid range. Moreover, our sensitivity studies showed that the desirable spacing of imperfections is simply proportional to the variation of the numerically obtained wavelength. In other words, although making any changes in the model (such as mesh size and material properties) can alter the wavelength, the objective range of $S_{imp} / λ_{s}$ remains invariant. Therefore, the range of $6 ⩽ S_{imp} / λ_{s} ⩽ 55$ is valid for any combinations of mesh size and material properties.

Convergence analysis and comparisons with theory

The results presented in sections “Direct compression” and “Indirect compression” illustrated that the incorporation of imperfections is able to trigger wrinkle formation even with relatively coarse finite element meshes. In this section a comprehensive convergence analysis is presented to assess the numerical approach and to compare the numerical and analytical solutions. The range of mesh refinement spans from a total of 50,000–1,600,000 elements, with the number increasing by a factor of two at each level. We have conducted extensive preliminary simulations and observed that the wavelength tends to become smaller as the mesh is refined, and as a result, more imperfections are required to obtain an objective wavy pattern. Considering the valid range of $S_{imp} / λ_{s}$ identified in the previous section, in the case of the coarsest mesh (50,000 elements) placing only one imperfection was sufficient to initiate uniform and stable wrinkles. More imperfections are required for the finer meshes. All numerical results on wavelength, critical buckling strain, and amplitude are listed in Tables 1 and 2.

The buckling wavelength is considered first. Figure 11 shows the variation of wavelength with the total number of elements. Three sets of numerical results are included—direct compression under generalized plane strain (GPE), direct compression under plane strain (PE), and indirect compression under generalized plane strain. The theoretical wavelength, based on equation (3), is also included in the figure for comparison. It is evident that, in both cases of direct and indirect compression, the wavelength converges to a value slightly lower than the theoretical solution as the mesh becomes sufficiently fine. As described in section “Brief overview of theory,” the theoretical value was derived under idealized assumptions. However, its difference with the converged numerical solution is small. Moreover, no significant difference exists between the numerical results of generalized and regular plane strain models in the case of direct compression.

Figure 11.

Variation of the simulated wrinkling wavelength with the number of elements for the cases of generalized plane strain (GPE) direct compression, plane strain (PE) direct compression, and indirect compression. The theoretical value is also included for comparison. (a) P3HT:PCBM film, (b) PEDOT:PSS film. The finite element meshes are all composed of four-noded quadrilateral elements in the present work.

Figure 12(a) and (c) shows the variation of simulated critical buckling strain with the number of elements, in the case of direct compression for both the generalized plane strain (GPE) and plane strain (PE) conditions. The theoretical value of $(e_{xx})_{cr}$ , based on equation (5), is included. The corresponding result for indirect compression is shown in Figure 12(b) and (d). Here, the critical strain is taken as the tensile strain $(e_{zz})_{cr}$ instead of the compressive strain $(e_{xx})_{cr}$ , with the theoretical value calculated from equation (5) and then substituted into equation (7). It is apparent that, in both cases of loading, the numerical solution is only slightly lower than the theoretical value as the mesh is refined. It is also noticed from Figure 12(a) and (c) that the PE numerical model results in a lower critical strain than the GPE model, suggesting that relaxing the out-of-plane (z) constraint caused by PE somewhat delays the onset of instability. In general, the results shown in Figures 11 and 12 have demonstrated the credibility of the present numerical approach in predicting mechanically induced buckling instability for thin film-thick substrate systems.

Figure 12.

Variations of the simulated critical buckling strain with the number of elements, for the cases of direct and indirect compression. (a, b) P3HT:PCBM film, (c, d) PEDOT:PSS film. Theoretical values are also included for comparison.

Attention is now turned to the wave amplitude displayed by the wrinkles. Note that the theoretical amplitude depends on both the applied strain and the critical buckling strain, equation (6). As shown in Figure 12, the critical buckling strain is different for each analysis. To ensure an objective comparison, a constant ratio of $e / e_{cr}$ is considered for all the simulations as well as the theoretical solution. (Note that according to equation (6), for any given geometry, the wave amplitude is uniquely determined by the $e / e_{cr}$ ratio.) Figure 13 shows the simulated amplitudes plotted against the number of elements, for the cases of direct and indirect compression. A constant ratio of $e / e_{cr} = 1.04$ , shortly after the onset of instability, was chosen for all cases in the figure. The amplitudes are also listed in Tables 1 and 2. As can be seen from Figure 13, the numerical results are close to one another and they converged to slightly below the theoretical value.

Figure 13.

Variations of simulated amplitudes with the number of elements, for the cases of generalized plane strain (GPE) under direct compression, plane strain (PE) under direct compression, and indirect compression. Note that the numerical results pertain to e/e_cr = 1.04 in all cases. The theoretical value is also included for comparison. (a) P3HT:PCBM film, (b) PEDOT:PSS film.

Further discussion

A numerical simulation framework was developed for predicting buckling instability of thin films attached to a thick compliant substrate. The incorporation of imperfections in the model enables simulations on a regular static finite element platform, with no additional treatment or special elements needed. The imperfections perturb the deformation field and trigger the buckling deformation modes. Since under normal circumstances the wrinkling wavelength is much smaller than typical specimen dimensions, a finite-size model is employed with periodic boundary conditions, which proves to be adequate in bringing out a uniform pattern of wrinkles once the critical buckling strain is reached. While the plane strain and generalized plane strain models are used in this study, the approach is readily applicable to any three-dimensional geometry.

In section “Brief overview of theory,” the theoretical basis for wrinkling was summarized. These analytical expressions were derived under the premise of a uniform sinusoidal surface contour. The present numerical simulations, however, can be used to visualize the progression of wrinkle formation. An example is given in Figure 14, where the generalized plane strain model during direct compression is shown under the applied compressive strain of (a) $e / e_{cr} = 0$ , (b) $e / e_{cr} = 0.9$ , (c) $e / e_{cr} = 1.0$ , (d) $e / e_{cr} = 1.04$ , and (e) $e / e_{cr} = 1.30$ . It can be seen that, shortly before the critical strain was reached, Figure 14(b), a small extent of disturbance has commenced, although the macroscopic load-displacement response has not displayed any deviation from that of the perfect model (e.g. Figure 4). At the instant of the critical buckling strain (Figure 14(c)), a localized (non-uniform) form of wavy pattern is evident. Note that this is only temporary, unlike the sustained irregular waveform due to the valid range of $S_{imp} / λ_{s}$ not being satisfied (as in Figure 10(a)). Despite the varying amplitude in Figure 14(c), the wavelength is a constant throughout the width of the specimen. The same wavelength persists with further straining, and in Figure 14(d) well developed uniform wrinkles have taken shape. This type of evolution process cannot be captured by the theoretical approach, hence illustrating the importance of numerical simulations in predicting surface instability phenomena. Moreover, it shows that depending on $e / e_{cr}$ , the wavy pattern can be both localized or uniform during the temporal evolution of surface contours irrespective of the imperfection distribution effect discussed in section “Effect of imperfection distribution” (note that theory predicts uniform buckling pattern regardless of the value of $e / e_{cr}$ ).

Figure 14.

Development of surface wrinkles for the case of generalized plane strain under direct compression (the model consists of 100,000 elements), when the applied compressive strain is (a) e/e_cr = 0, (b) e/e_cr = 0.9, (c) e/e_cr = 1.0, (d) e/e_cr = 1.04, (e) e/e_cr = 1.30. Images cover only a portion of the surface region, with a displacement scaling factor of 10.

The present technique is not only of significance in numerical studies but also important from the physical point of view. One may draw an analogy between the embedded imperfections and interface defects in an actual specimen. Deviations from a uniform wavy pattern observed in real physical experiments may be rationalized by the numerical simulations with different imperfection placements. Uncommon forms and irregularities of wrinkles can potentially be captured even with a fully linear-elastic film/substrate system. The current approach adds to the existing numerical capability and can be further employed to study how buckling instability is influenced by the thickness ratio and elastic property relationship between the film and substrate materials. It is also well suited as a design tool for the case of multilayer thin films on a substrate. In addition, the approach can potentially be combined with various constitutive behaviors such as viscoelasticity, plasticity, damage, and delamination.

Conclusion

A computational modeling approach was developed to predict the evolution of buckling (wrinkling) of relatively stiff thin coatings above a thick compliant substrate. Buckling instability can be triggered by locally perturbed deformation field through the incorporation of imperfections. Wrinkles form under directly applied compressive strain. When the Poisson’s ratio of the substrate is greater than that of the film, applied tension also induces lateral compressive strain in the film, leading to wrinkling. Instabilities during both of these direct and indirect compression are successfully predicted by the simulations. Through the systematic analyses, a valid range of imperfection density, being able to generate an “objective” wrinkle pattern, was identified. With a sufficiently fine mesh, the simulated wavelength, amplitude, and critical buckling strain converge to near the respective analytical solutions. The current approach does not involve any multi-step or iterative numerical procedure, and is straightforward, computationally efficient, and can be implemented with any finite element analysis platform using regular continuum elements. It is well suited for predicting buckling instability in the design and analysis of thin film systems over a wide range of material and geometric conditions.

Footnotes

Handling Editor: Guian Qian

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors would like to acknowledge NASA EPSCoR CAN (grant #: 80NSSC17M0050), New Mexico Space Grant Consortium, and NASA’s Space Grant College and Fellowship Program for supporting this study.

ORCID iDs

Siavash Nikravesh

Yu-Lin Shen

References

Khang

Jiang

Huang

et al . A stretchable form of single-crystal silicon for high-performance electronics on rubber substrates. Science 2006; 311: 208–212.

Kim

Ahn

Choi

et al . Stretchable and foldable silicon integrated circuits. Science 2008; 320: 507–511.

Shin

Wang

et al . Curvilinear electronics formed using silicon membrane circuits and elastomeric transfer elements. Small 2009; 5: 2703–2709.

Kaltenbrunner

Sekitani

Reeder

et al . An ultra-lightweight design for imperceptible plastic electronics. Nature 2013; 499: 458.

Zhu

Highly conductive and stretchable silver nanowire conductors. Adv Mater 2012; 24: 5117–5122.

Koo

Jeong

Araoka

et al . Light extraction from organic light-emitting diodes enhanced by spontaneously formed buckles. Nat Photonics 2010; 4: 222.

Lipomi

Tee

BCK

Vosgueritchian

et al . Stretchable organic solar cells. Adv Mater 2011; 23: 1771–1775.

Kim

Pégard

et al . Wrinkles and deep folds as photonic structures in photovoltaics. Nat Photonics 2012; 6: 327.

Wang

Zhang

Yang

et al . Enhanced light scattering effect of wrinkled transparent conductive ITO thin film. RSC Adv 2017; 7: 25483–25487.

10.

Hsieh

Chen

Shih

et al . Mechanically robust, stretchable organic solar cells via buckle-on-elastomer strategy. Org Electron 2018; 53: 339–345.

11.

Zhang

Zheng

Fang

et al . Enhancement of silicon-wafer solar cell efficiency with low-cost wrinkle antireflection coating of polydimethylsiloxane. Sol Energ Mat Sol C 2018; 181: 15–20.

12.

Tang

Guo

Chen

et al . Wrinkled Ag nanostructured gratings towards single molecule detection by ultrahigh surface Raman scattering enhancement. Sensor Actuat B: Chem 2015; 218: 145–151.

13.

Kim

Crosby

AJ.

Solvent-responsive surface via wrinkling instability. Adv Mater 2011; 23: 4188–4192.

14.

Zhao

et al . Highly stretchable carbon nanotube transistors enabled by buckled ion gel gate dielectrics. Appl Phys Lett 2015; 107: 751.

15.

Ohzono

Monobe

Yamaguchi

et al . Dynamics of surface memory effect in liquid crystal alignment on reconfigurable microwrinkles. Appl Phys Lett 2009; 95: 014101.

16.

Chen

Yang

Wrinkling instabilities in polymer films and their applications. Polym Int 2012; 61: 1041–1047.

17.

Moon

Lee

Sun

et al . Wrinkled hard skins on polymers created by focused ion beam. P Natl Acad Sci USA 2007; 104: 1130–1133.

18.

Sun

Xia

Moon

et al . Folding wrinkles of a thin stiff layer on a soft substrate. Proc R Soc A 2012; 468: 932–953.

19.

Yan

Zhang

Wrinkling of structured thin films via contrasted materials. Soft Matter 2016; 12: 3937–3942.

20.

Lam

Sim

Zhu

et al . The effect of continuous wavy micropatterns on silicone substrates on the alignment of skeletal muscle myoblasts and myotubes. Biomaterials 2006; 27: 4340–4347.

21.

Vandeparre

Gabriele

Brau

et al . Hierarchical wrinkling patterns. Soft Matter 2010; 6: 5751–5756.

22.

Volynskii

Bazhenov

Lebedeva

et al . Mechanical buckling instability of thin coatings deposited on soft polymer substrates. J Mater Sci 2000; 35: 547–554.

23.

Biot

MA.

Mechanics of incremental deformations. New York: John Wiley & Sons, 1965.

24.

Groenewold

Wrinkling of plates coupled with soft elastic media. Physica A 2001; 298: 32–45.

25.

Chung

Nolte

Stafford

CM.

Surface wrinkling: a versatile platform for measuring thin-film properties. Adv Mater 2011; 23: 349–368.

26.

Stafford

Harrison

Beers

et al . A buckling-based metrology for measuring the elastic moduli of polymeric thin films. Nat Mater 2004; 3: 545.

27.

Stafford

Vogt

Harrison

et al . Elastic moduli of ultrathin amorphous polymer films. Macromolecules 2006; 39: 5095–5099.

28.

Cao

Hutchinson

JW.

From wrinkles to creases in elastomers: the instability and imperfection-sensitivity of wrinkling. Proc R Soc A 2012; 468: 94–115.

29.

Zang

Zhao

Cao

et al . Localized ridge wrinkling of stiff films on compliant substrates. J Mech Phys Solid 2012; 60: 1265–1279.

30.

Zheng

Cao

et al . Wrinkling of a stiff film resting on a fiber-filled soft substrate and its potential application as tunable metamaterials. Extreme Mech Lett 2017; 11: 121–127.

31.

Saha

SK.

Sensitivity of the mode locking phenomenon to geometric imperfections during wrinkling of supported thin films. Int J Solids Struct 2017; 109: 166–179.

32.

Bayat

Gordaninejad

Switching band-gaps of a phononic crystal slab by surface instability. Smart Mater Struct 2015; 24: 075009.

33.

Huck

Bowden

Onck

et al . Ordering of spontaneously formed buckles on planar surfaces. Langmuir 2000; 16: 3497–3501.

34.

Mei

Landis

Huang

Concomitant wrinkling and buckle-delamination of elastic thin films on compliant substrates. Mech Mater 2011; 43: 627–642.

35.

ABAQUS. ABAQUS Analysis user’s manual 2016. Providence, RI: Dassault Systems Simulia Corp, 2016.

36.

Zheng

Cao

et al . Surface wrinkling of nanostructured thin films on a compliant substrate. Comput Mater Sci 2010; 49: 767–772.

37.

Cao

Zheng

Jia

et al . Wrinkling and creasing of a compressed elastoplastic film resting on a soft substrate. Comput Mater Sci 2012; 57: 111–117.

38.

Andrews

Massabò

Delamination in flat sheet geometries with material imperfections and thickness variations. Compos Part B: Eng 2008; 39: 139–150.

39.

Pulliam

Hoover

Tiparti

et al . Development of self-powered strain sensor using mechano-luminescent ZnS: Cu and mechano-optoelectronic P3HT. In: Sensors and smart structures technologies for civil, mechanical, and aerospace systems, vol. 101680, Bellingham, WA, 12 April 2017. Bellingham, WA: SPIE.

40.

Pulliam

Hoover

Ryu

. Multifunctional mechano-luminescent-optoelectronic composites for self-powered strain sensing. In: Proceedings of the ASME 2017 conference on smart materials, adaptive structures and intelligent systems, Snowbird, UT, 18–20 September 2017. New York: ASME.

41.

Ryu

Mongare

Corrugated photoactive thin films for flexible strain sensor. Materials 2018; 11: 1970.

42.

Lang

Naujoks

Dual

Mechanical characterization of PEDOT: PSS thin films. Synthetic Met 2009; 159: 473–479.

43.

Tahk

Lee

Khang

DY.

Elastic moduli of organic electronic materials by the buckling method. Macromolecules 2009; 42: 7079–7083.

44.

Shen

YL.

Constrained deformation of materials: devices, heterogeneous structures and thermo-mechanical modeling. New York: Springer, 2010.