Numerical Investigation of Dynamic Wrinkling Behaviors in Stiff-Film/PDMS-Substrate Structure

TL;DR

This paper studies how stiff films on flexible substrates behave dynamically, aiming to improve the design of flexible electronics.

Contribution

A novel numerical approach is introduced to analyze dynamic wrinkling behaviors in stiff-film/PDMS-substrate structures.

Findings

Dynamic equations were derived using energy principles for wrinkled thin film/substrate systems.

Pre-strain and substrate stiffness significantly affect system stability and energy.

Numerical simulations reveal the impact on potential energy and chaos critical values.

Abstract

Thin film/substrate structures based on the principle of buckling mechanics exhibit both excellent stretchability and mechanical stability, and they have been recognized as a critical configuration in the design of flexible electronic devices. During application, flexible electronic devices are usually subjected to complex dynamic environments. Therefore, it is of great significance to investigate the dynamic behavior of thin film/substrate structures for the design of flexible electronic devices. The bending energy, membrane energy, and kinetic energy of the thin film and the elastic energy of the substrate were calculated. On this basis, the dynamic equation of the thin film/substrate structure with a checkerboard wrinkled pattern was derived by applying the principle of minimum energy combined with the Lagrangian function. Numerical simulations were conducted on the system to analyze the effect of pre-strain and the Young’s modulus of substrate on the system’s potential energy function, simulate the temporal response of the system’s dynamic behavior, and investigate the influences of pre-strain and the Young’s modulus of substrate on system stability and the chaos critical value. Theoretical support is expected to be provided for the design of two-dimensional (2D) thin film/substrate structures through this research.

Figures

Full text

1. Introduction

Based on the principle of buckling mechanics, the wrinkled structure of a thin, stiff film and polydimethylsiloxane(PDMS) substrate possesses both excellent stretchability and mechanical stability. This thin film/substrate structure has become a key configuration in the design of flexible electronic devices, demonstrating broad application potential in materials science and emerging functional device fields. Therefore, it is necessary to study the dynamic response of this structure under dynamic loading, as it not only reveals the instability evolution of this structure but also provides theoretical support for improving the structural reliability and service performance of flexible electronics. This is of great significance for promoting the design and engineering application of such devices [1,2,3].

With respect to the static buckling problem of film/substrate structures, extensive investigations have been carried out by numerous scholars in the field. In the research on the formation and evolution of film wrinkles, remarkable advancements have been achieved in multiple studies. Based on the principle of energy minimization, Huang et al. [4] investigated the wrinkling behavior in a multilayer structure consisting of a rigid film, a soft substrate, and a rigid support. To accurately predict the instability behavior of practical devices under tension, Li et al. [5] developed an analytical model capable of distinguishing between global buckling and local wrinkling for substrates with finite thickness. To reveal the critical influence of substrate nonlinearity on advanced wrinkling modes, Zang et al. [6] analyzed the buckling behavior of rigid films on pre-stretched elastic substrates. A theoretical framework for film buckling on soft substrates has been developed by Pan et al. [7], which clarifies the transition path and geometric laws from a flat state to a blister state under uniaxial and equibiaxial conditions. For the first time, the crystallographic anisotropy of films and the viscoelasticity of substrates have been coupled by Im et al. [8], and the entire process from the initiation of instability to the formation of equilibrium patterns has been elaborated. To investigate the factors influencing the evolution of wrinkles, Jin et al. [9] examined how the relative magnitude between the critical strain for wrinkling and the critical strain for wrinkling initiation affects the formation and evolution of wrinkles.

Further investigations into two-dimensional buckling have been conducted by numerous scholars. In terms of the buckling behavior in the hard film/soft substrate system driven by thermal stress, it has been investigated by Song et al. [10] that three buckling modes, namely one-dimensional, checkerboard, ordered herringbone modes, and a prediction model based on the film and various parameters, have been established. For the biaxial structure of silicon nanomembranes with two-dimensional buckling, Choi et al. [11] have studied the geometric shape and strain response of the material along different directions. By applying orthogonal strain to the interlayer of an elastomer, Pellegrino et al. [12] developed a method for generating two-dimensional dual-frequency patterns through the superposition of one-dimensional single-frequency wrinkles. Based on the buckling behavior of elastic films under residual compression, the interface boundary conditions of film/substrate have been analyzed, and linear stability analysis has been performed by Audoly et al. [13]. To optimize the design of functional surfaces in flexible electronics, Huang et al. [14] identified that the tension–compression asymmetry of soft substrates serves as the key factor governing pattern evolution. To investigate the difference in buckling behavior under uniaxial and biaxial compression, Peterson et al. [15] conducted a quantitative comparison of the buckling behavior in silicon-germanium thin films under these two loading conditions and contrasted their respective buckling propagation rates.

In the research on the nonlinear behavior of film buckling, significant contributions have been made by numerous scholars to the exploration of its analytical or numerical solutions. Regarding the formation mechanism of film wrinkle patterns on planar and three-dimensional structured substrates, a hybrid analytical-numerical method has been proposed by Cutolo et al. [16] to simulate the wrinkle evolution process. Based on three-dimensional nonlinear finite element analysis, Xu et al. [17] explored the pattern formation process of rigid films on flexible substrates, and they adopted the asymptotic numerical method (ANM) to achieve robust path following and bifurcation detection. Based on a neo-Hookean substrate under uniaxial compression, the buckling behavior of neo-Hookean films has been studied by Fu et al. [18], and it has been revealed that there exists a class of local solutions with asymptotically decaying tails in the subcritical region. Under two-dimensional strain conditions, Nikravesh et al. [19] proposed a finite element method to simulate the wrinkle formation in thin films on flexible substrates, where random initial imperfections were introduced to trigger instabilities. Based on the finite deformation theory, the analytical solution for the highly nonlinear behavior of the mechanical model of buckled films on soft substrates has been obtained by Song et al. [20] through perturbation analysis. In the field of flexible electronic engineering, the precise regulation of surface patterns of film/substrate structures serves as a crucial design basis for the realization of device functionalities. To explore the variation in buckling patterns, Audoly et al. [21] investigated the buckling behavior of compressed elastic films and revealed the evolutionary path of buckling patterns during the accumulation of residual stresses. To reveal the energy characteristics of different wrinkles, Cai et al. [22] compared the energetic properties of various periodic wrinkle patterns formed in thin rigid films on biaxially compressed compliant substrates over a wide range of overstress. With respect to the influence of surface patterns on the energy of film/substrate structures, the effect of equibiaxial compressive stress on modes has been studied by Chen et al. [23], and it has been revealed that the herringbone mode is the configuration with the lowest energy. Based on the multi-scale method with slowly varying Fourier coefficients, Xu et al. [24] investigated the pattern formation of rigid thin films on soft substrates and realized the quantitative description and bifurcation analysis of sinusoidal and checkerboard patterns through a multiscale finite element framework. To predict buckling orientation, Song et al. [25] developed a theoretical model for chevron buckling patterns in anisotropic thin films on elastic substrates.

When film/substrate electronic devices operate in complex environments, it is of great significance to investigate their nonlinear dynamic behaviors to ensure the performance and reliability of the devices. Based on the analytical solution of Jacobi elliptic functions, Ou et al. [26] investigated the dynamic buckling behavior of thin films on soft substrates and determined the dynamic critical load by incorporating the Budiansky–Roth criterion. By employing the symplectic method, the nonlinear dynamic behaviors of piezoelectric film/soft substrate structures, as well as the random vibration behaviors under white noise excitation, have been analyzed by Wang et al. [27] and Bi et al. [28]. Considering the influence of the damping effect on film/substrate structures, a multi-level coupled time-varying parameter dynamic model has been developed by Zhou et al. [29] to analyze the nonlinear dynamic buckling of film/viscoelastic substrate structures under time-varying excitation. By introducing the dynamic boundary effect of sub-nanoscale van der Waals interactions, a dynamic analysis model for layered nanofilms on substrates has been established by Dong et al. [30]. Based on the small deformation theory, the governing equations have been derived by Zhang et al. [31] through the Lagrangian function and Euler-Lagrange equations, and the analytical solution of the response under linear loads has been obtained. To investigate the dynamic characteristics of buckling piezoelectric energy harvesters on pre-strained compliant substrates, Bi et al. [32] established governing equations and solved for steady-state responses, revealing that pre-strain and substrate modulus can effectively modulate high-energy orbits. Based on the nonlinear free vibration of beams on Winkler foundations, Ma et al. [33] established nonlinear equations of motion accounting for soil-structure dynamic coupling, and they obtained both linear and nonlinear natural frequencies and mode shapes of the beams through eigenvalue analysis and the multiple scales method. Based on the buckling behavior of thin films under nanoscale surface effects, Wang et al. [34] analyzed the influence of surface elasticity and residual surface tension on static and dynamic characteristics, and proposed design strategies for buckled film interconnects in flexible electronics. For the three-layer structure, Wang et al. [35] established static and dynamic theoretical models for the film/interlayer/substrate system, and the buckling amplitude and nonlinear vibration frequency have been analyzed. To elucidate the role of interlayers in regulating wrinkle stability, Bi et al. [36] investigated their influence on the system’s buckling behavior, distinguishing between global buckling and local wrinkling modes.

It has been widely documented that in film/substrate structure, two-dimensional buckling modes, including herringbone and checkerboard patterns, have been extensively investigated under static conditions. Nevertheless, it should be emphasized that their mechanical response under dynamic loading remains inadequately explored. Consequently, the checkerboard buckling mode is selected as the research focus in this study, and its vibration characteristics are systematically examined. In Section 2 of the paper, the dynamic equations of the thin film/substrate structure with a checkerboard wrinkle pattern are derived. In Section 3, numerical simulations are performed on the dynamic system, and the effects of parameters such as pre-strain on the dynamic behaviors of the system are discussed. Finally, in Section 4, conclusions are drawn.

2. Materials and Methods

As shown in Figure 1, for the checkerboard film/substrate structure, $[eqn]$ denotes the thickness of the film, $[eqn]$ and $[eqn]$ represent the Poisson’s ratios of the film and the substrate, respectively, and $[eqn]$ and $[eqn]$ denote the Young’s moduli of the film and the substrate, respectively. In this paper, it is assumed that $[eqn]$ and $[eqn]$ are the equivalent Young’s modulus of the film and the equivalent Young’s modulus of the substrate, respectively. The $[eqn]$ and $[eqn]$ axes point to the principal strain directions of the biaxial load, respectively, and the $[eqn]$ direction represents the thickness direction of the film. The symbol explanations in the text are also uniformly presented in Table A1 of Appendix A. The thickness of the substrate is sufficiently larger than that of the film, and it is regarded as a semi-infinite space in the model.

According to Huang et al. [4] and Song et al. [10], assuming the out-of-plane displacement function is given by

$$[eqn]$$

where $[eqn]$ denotes the out-of-plane displacement of the checkerboard membrane structure, $[eqn]$ represents the time-varying wrinkle amplitude of this structure, and $[eqn]$ and $[eqn]$ are the characteristic wave numbers along the $[eqn]$ and $[eqn]$ directions, respectively. The bending energy density of the film is expressed as Huang et al. [4] and Song et al. [10]

$$[eqn]$$

By computing the second-order partial derivatives of $[eqn]$ with respect to $[eqn]$ and $[eqn]$ , along with the mixed partial derivative with respect to $[eqn]$ , and substituting these into Equation (2), the bending energy density can be expressed in terms of $[eqn]$ , $[eqn]$ , and $[eqn]$ as

$$[eqn]$$

Then, the bending energy of the film is

$$[eqn]$$

Substituting Equation (3) into Equation (4) yields

$$[eqn]$$

Since the shear stress at the film/substrate interface has no effect on the wavelength and amplitude of the structure, i.e., the shear stress at the interface is zero, it follows that the shear stress in different directions can be obtained as

$$[eqn]$$

where $[eqn]$ are the membrane forces to the membrane strains, and can be obtained as

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ , $[eqn]$ denotes the strain tensor of the film, $[eqn]$ , and $[eqn]$ represents the pre-strain of the thin film in each direction.

The relationships between the in-plane displacement of the film $[eqn]$ , $[eqn]$ and $[eqn]$ , $[eqn]$ can be obtained, leading finally to

$$[eqn]$$ $$[eqn]$$

The membrane energy density of the film is

$$[eqn]$$

Substituting Equations (1), (9) and (10) into Equation (8), $[eqn]$ can then be determined.

$$[eqn]$$

The membrane energy can then be calculated as

$$[eqn]$$

The strain energy density of the substrate is

$$[eqn]$$

where $[eqn]$ denotes the stress of the substrate.

And the strain energy of the substrate is obtained as

$$[eqn]$$

This study employs the principle of energy minimization. By taking the partial derivatives of the total energy $[eqn]$ with respect to $[eqn]$ , $[eqn]$ , $[eqn]$ , and setting them to zero, the wave numbers ( $[eqn]$ ) are finally obtained as [10]

$$[eqn]$$

The total kinetic energy of the film/substrate structure is

$$[eqn]$$

where $[eqn]$ is the density of the film. The work performed by external forces is

$$[eqn]$$

where $[eqn]$

The Lagrangian function of this structure is $[eqn]$ , and the Euler equation is given by

$$[eqn]$$

Taking linear damping into account for this system, the dynamic equation of this structure can be derived as

$$[eqn]$$

where $[eqn]$ is the damping term. $[eqn]$ is the external excitation, $[eqn]$ , $[eqn]$ is the amplitude of the external excitation, $[eqn]$ is the frequency of the external excitation, and

$$[eqn]$$

For convenience, the following nondimensional parameters of this structure are introduced,

$$[eqn]$$

Substitution of the above equations into the system (20) yields

$$[eqn]$$

where $[eqn]$ .

The potential energy function of the system (21) is

$$[eqn]$$

3. Results and Discussion

3.1. Comparison with Experiment

To enable a more specific comparison with previous experiments, this study presents plots of both wavelength and amplitude versus pre-strain at an evolution time of 50 s, which are compared with the experimental results reported in Pellegrino et al. [12]. The comparison is illustrated in Figure 2, and the parameters used in this work are listed as follows: $[eqn]$ , $[eqn]$ , $[eqn]$ .

Figure 2 is compared with the amplitude and wavelength at an evolution time of 50 s in Pellegrino et al. [12]. As shown in Figure 2, when the pre-strain is 10%, the theoretical wavelength is 2.36 $[eqn]$ , while the actual value reported in Pellegrino et al. [12] is 2.3 $[eqn]$ , resulting in a relative error of 2.6%. When the pre-strain is 20%, the theoretical wavelength is 2.15 $[eqn]$ , compared to the actual value of 2 $[eqn]$ in Pellegrino et al. [12], yielding a relative error of 7.5%. Both cases show relatively small discrepancies. These results further corroborate the accuracy of the theoretical derivations presented in this paper.

3.2. The Effect of the Pre-Strain

To illustrate the dynamic energy characteristics of the buckling system, we examine the potential energy function of a structure exhibiting a checkerboard buckling mode. The influence of the pre-strain parameters on this potential energy function is explored, with the results visualized in Figure 3. The other physical parameters of the structure with a checkerboard buckling mode are defined as Song et al. [10] $[eqn]$ , $[eqn]$ , $[eqn]$ $[eqn]$ , $[eqn]$ .

As shown in Figure 3, the system’s potential energy function possesses a single saddle point and two stable equilibrium points. In Figure 3a, with $[eqn]$ , an increase in $[eqn]$ causes the positions of the stable equilibria to shift away from the initial dimensionless displacement. Concurrently, the potential barrier height increases monotonically. Furthermore, the distance between the two stable equilibria and the separation distance between the wells would increase with increasing $[eqn]$ . Similarly, Figure 3b demonstrates that for $[eqn]$ , increasing $[eqn]$ also shifts the stable equilibria away from the initial displacement, elevates the potential barrier, and increases the separation between the two potential wells. These observations indicate that the pre-strain parameters $[eqn]$ and $[eqn]$ exert a symmetrical impact on the potential energy landscape. Therefore, we only need to discuss the influence of the value of $[eqn]$ on the dynamic behavior.

Driven by the aim to preliminarily investigate the influence of the parameter $[eqn]$ on the system’s dynamic behavior, we present the numerical results in Figure 4. The parameter $[eqn]$ was varied across the values 0.03, 0.05, and 0.07.

Figure 4 shows three different oscillating modes with $[eqn]$ , $[eqn]$ , $[eqn]$ . For $[eqn]$ , the system exhibits a periodic response with one dominant cycle. At this point, combined with Figure 3, the potential well depth is relatively small, the energy required to cross the potential well is relatively low, and inter-well vibration is prone to occur. When $[eqn]$ increases to 0.05, the system transitions into a chaotic state, characterized by an aperiodic and irregular time-history response of the dimensionless displacement. In this case, the system can vibrate in cross-well. Further increasing $[eqn]$ to 0.07 restores a periodic oscillation. As $[eqn]$ increases, the potential well depth increases, and the energy required to cross the potential well becomes larger; at this point, the system undergoes intra-well vibration.

In order to investigate the influence of the pre-strain parameters $[eqn]$ and $[eqn]$ on the critical range for chaotic vibration in the checkerboard buckling mode system, Figure 5 is presented. The parameters for the system are set as, $[eqn]$ , $[eqn]$ , $[eqn]$ .

As shown in Figure 5a, the system exhibits a stable 1-periodic vibration when it lies between 0.01 and 0.036. At this point, since $[eqn]$ is small, the potential well depth is relatively low (as shown in Figure 3). As the energy can easily overcome the potential barrier, the structure undergoes inter-well vibration. When $[eqn]$ increases to 0.037, the structure undergoes chaotic vibration, and this process terminates when $[eqn]$ = 0.059. As $[eqn]$ further increases, the system recovers to a stable 1-periodic motion. At this point, since $[eqn]$ is large, the potential well depth is high (as shown in Figure 3). Because the energy is insufficient to overcome the potential barrier, the structure undergoes intra-well vibration. Similarly, system response to variations in $[eqn]$ are presented in Figure 5b. A comparison between Figure 5a,b reveals that the influences of the pre-strain parameters $[eqn]$ and $[eqn]$ on the dynamic stability of the system exhibit a clear symmetric characteristic.

By investigating the effect of pre-strain on the system, it can be concluded that as the pre-strain increases, the potential well depth increases, and the system undergoes a transition from period-1 to chaos and then back to period-1. When the pre-strain is small, the structure is prone to inter-well vibration because the energy is sufficient to cross the potential barrier; when chaos occurs, the structure vibrates in the cross-well. When the pre-strain is large, the structure is prone to intra-well vibration due to insufficient energy to cross the potential barrier. The results are consistent with the findings on one-dimensional buckling structures reported by Bi et al. [32]. While quantitative discrepancies exist in numerical simulations due to differences in parameters, qualitative similarities remain evident. For instance, both studies indicate that as pre-strain increases, dynamic behavior evolves from period-1 inter-well oscillations to chaos, eventually returning to period-1 intra-well oscillations. This observation provides a feasible strategy for modulating dynamic behaviors through pre-strain adjustment.

3.3. The Effect of the Young’s Moduli of the Substrate

To investigate the influence of the Young’s modulus of the substrate on the potential energy function of the system. The influence of the Young’s modulus of substrate parameters on this potential energy function is explored, with the results visualized in Figure 6. The physical parameters of the structure are defined as Song et al. [10]: $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ .

It can be observed from Figure 6 that with the increase in the Young’s modulus of the substrate, the position of the stable equilibrium point is shifted closer to the initial dimensionless displacement, and the height of the potential barrier is correspondingly reduced. Meanwhile, the distance between the two stable equilibrium points is decreased, and the inter-well distance is also diminished. Consequently, the energy required for the system to perform inter-well vibration is lowered, and it is more likely that inter-well vibration occurs.

To preliminarily investigate the influence of different $[eqn]$ values on the system’s dynamic behavior, we performed numerical simulations with the results presented in Figure 7. The parameters of the checkerboard buckling mode system were set as $[eqn]$ , $[eqn]$ , $[eqn]$ . The parameter $[eqn]$ takes the values 1 MPa, 2 MPa, 3 MPa.

Figure 7 reveals three distinct dynamic regimes. For $[eqn]$ , the system exhibits periodic-1 oscillations. Combined with Figure 6, due to the small Young’s modulus of the substrate, the potential well depth is relatively large. The energy of the system is insufficient to overcome the potential barrier, and the structure undergoes intra-well vibration at this point. As $[eqn]$ increases to 2 MPa, the system bifurcates into a chaotic state, characterized by an aperiodic and irregular displacement response. The structure vibrates in the cross-well. Finally, at $[eqn]$ , the dynamics revert to periodic-1 vibrations. Combined with Figure 6, at this point, since the Young’s modulus of the substrate is relatively large, the potential well depth is relatively small. The system’s energy can easily overcome the potential barrier; thus, the structure undergoes intra-well vibration.

It aims to explore the influence of the Young’s modulus of the substrate $[eqn]$ of the system on the critical value range of chaotic vibration in the checkerboard buckling mode vibration system. Accordingly, Figure 8 was plotted in this study, where the physical parameters of the checkerboard buckling mode system are $[eqn]$ , $[eqn]$ , $[eqn]$ .

It can be observed from Figure 8 that a bifurcation diagram corresponding to the Young’s modulus of the substrate $[eqn]$ in the range of 1~4 has been constructed in this study. When $[eqn]$ increases from 1 MPa to 1.56 MPa, the system is maintained in a stable period-1 state and performs inter-well vibration during this process. As $[eqn]$ continues to increase, the system undergoes chaotic vibration, and the chaotic behavior terminates when $[eqn]$ rises to 3.28 MPa. During this process, the system can vibrate in cross-well. With the further increase of $[eqn]$ , the system recovers to a period-1 vibration. Meanwhile, the energy is sufficient to overcome the potential barrier, and the structure undergoes inter-well vibration.

By investigating the effect of the Young’s modulus of the substrate on the system, it can be concluded that as the Young’s modulus increases, the potential well depth decreases, and the system undergoes a transition from period-1 to chaos and then back to period-1. When the Young’s modulus of the substrate is small, the structure is prone to intra-well vibration because the energy is insufficient to overcome the potential barrier; when chaos occurs, the structure vibrates in the cross-well. When the pre-strain is large, the structure is prone to inter-well vibration due to sufficient energy to overcome the potential barrier. The present results are qualitatively consistent with those reported by Bi et al. [32] for one-dimensional buckling structures, though quantitative deviations exist due to parametric differences. Specifically, both studies demonstrate that as the substrate Young’s modulus increases, the dynamic behavior transitions from period-1 intra-well oscillation to chaos and eventually returns to period-1 inter-well oscillation. This agreement suggests a viable strategy for tailoring dynamic response regimes through changing the Young’s modulus of the substrate.

4. Conclusions

In this study, the governing equations for the nonlinear vibration of rigid-film structures with a checkerboard wrinkle pattern are established. The fourth-order Runge–Kutta method is employed to solve these governing equations, and the dynamic responses of the system over time under different physical parameters are investigated. Additionally, the variation in the potential energy function under different pre-strains and Young’s modulus of the substrate is analyzed, along with the changes in the critical values and stability of the system’s chaotic vibration under various parameters. Numerical simulations reveal that the excitation frequency and the system’s pre-strain can alter the dynamic behavior of the system.

The pre-strain and its influence on the potential energy function exhibit symmetry; as the pre-strain increases, the distance between the two stable equilibrium points of the potential energy function increases, and the inter-well distance and potential barrier height also increase in accordance.As the Young’s modulus of the substrate increases, the distance between the two stable equilibrium points of the potential energy function decreases, and the inter-well distance and potential barrier height also decrease accordingly.When the Young’s modulus of the substrate is small, the energy cannot overcome the potential barrier, and the system undergoes intra-well vibration; when the Young’s modulus of the substrate is large, the energy can overcome the potential barrier, and the system will exhibit inter-well vibration.Under the parameters adopted in this study, the stability of the system is affected by the selection of the Young’s modulus of the substrate and the system pre-strain. For instance, the system may transition from period-1 vibration to chaotic vibration.

This study on the nonlinear dynamics of checkerboard-patterned thin film/substrate structures can effectively prevent flexible electronics from entering unstable dynamic regions due to inappropriate selection of material dimensions and parameters during operation, thereby ensuring reliability and robustness in their performance. This study aims to provide theoretical support for the design of two-dimensional checkerboard rigid film structures.