Vibration Analysis of Post-Buckled Thin Film on Compliant Substrates

Abstract

Buckling stability of thin films on compliant substrates is universal and essential in stretchable electronics. The dynamic behaviors of this special system are unavoidable when the stretchable electronics are in real applications. In this paper, an analytical model is established to investigate the vibration of post-buckled thin films on a compliant substrate by accounting for the substrate as an elastic foundation. The analytical predictions of natural frequencies and vibration modes of the system are systematically investigated. The results may serve as guidance for the dynamic design of the thin film on compliant substrates to avoid resonance in the noise environment.

1. Introduction

As is well known, buckling instability appears in the film/substrate system when subjected to external loading, which is utilized by stretchable electronics to achieve stretchability. Recent advances in the buckling of thin films on a compliant substrate enable high stretchability for stretchable electronics [1,2,3,4,5,6,7,8], which can be used as human vital signs monitors [9,10,11,12,13], sensitive electronic skins [14,15,16], eye-like digital cameras [17,18,19], tunable phase optics [20,21], and tunable optical window [22,23].

In these systems, controlled buckling is generated in thin films deposited onto prestrained substrate after releasing the prestrain in the substrate [24,25,26,27,28,29,30,31,32]. The post-buckling morphology can be global buckling and local buckling with different modulus of substrate [24,25]. Bowden et al. [26] showed the buckling morphology surface within 100 nm–100 μm range in the film. Furthermore, Khang et al. [28] illustrated wave length and amplitude with specific external pre-strains through the experiments and the analytical model. A series of analytical models were developed to precisely predict the post-buckling shape with different material, geometric, and loading conditions [29,30].

However, in operation conditions, the system undergoes complex noise environment with electrical, mechanical, and thermal loadings. The dynamic behavior of the post-buckled system is one of the most important problems to solve for the applications of stretchable electronics. Many investigations on vibration analysis of post-buckling structures have been carried out in recent years [33,34,35,36,37,38,39,40,41,42]. Nayfeh et al. [33] presented an exact solution for the post-buckling configurations of beams, investigates the lowest natural frequencies of vibration around each of the first three buckled configurations beams with various boundary conditions, and discusses the stability of the beam. Emam et al. [34] derived approximate analytical solutions for the nonlinear free vibrations of laminated composite beams in prebuckling and postbuckling. Cristobal et al. [35] investigated dynamic and impact behaviors of glass fiber reinforced polymer composites. Neukirch et al. [36] investigated the vibration of post-buckled rods clamped at both ends. Emam and Nayfeh studied the non-linear response of the buckled beam considering internal resonances [37]. Furthermore, Ansari et al. [38] derived an analytical model to describe the dynamic behaviors of a postbuckled microscale functionally graded (FG) beam with modified coupled stress theory. The response of nonlinear vibration of postbuckled beams has also been investigated numerically [39]. Focusing on stretchable electronics, Wang et al. [41] established an analytical model to predict the dynamic behaviors of buckled thin films considering surface effects, and Wang et al. [42] derived an explicit analytical solution for the vibration mode and the linear natural frequency of a buckled ribbon, which was verified by finite element method (FEM) and experiments. Tseng et al. [43] derived a solution to the linear natural frequency of the first-order vibration mode for a buckled ribbon with fixed ends. The former papers have paid much attention to buckling and dynamic behaviors of thin film. However, they are not directly applicable in solving the vibration problem for flexible electronics, in which the compliant substrates beneath the devices vibrate together with thin films. Currently, there are very scant investigations of analytical solution on the vibration analysis of post-buckled thin films on a compliant substrate. With only the assistance of numerical simulations, it is difficult to understand the mechanical nature of the film/substrate system and is time-consuming in practical device design and optimization.

Therefore, this paper aims to establish an analytical model to illustrate the dynamic behaviors of this special system. However, dynamic loading, temperature change, vibration environment, and other factors affect the mechanical properties of the substrate, which can make the substrate viscoelastic instead of just being simplified as springs. This paper mainly focuses on hyperelastic substrate such as polydimethylsiloxane (PDMS) with small temperature change and low vibration frequency. In such cases, viscoelasticity of the substrate has negligible influence on the vibration analysis. The outline of this paper is as follows. Section 2 illustrates the analytical modeling for buckling analysis, and the vibration analysis is shown in Section 3. Results and discussion are given in Section 4. The main conclusions are presented in Section 5.

2. Buckling Analysis

Figure 1a illustrates the thin film with thickness of h_f on top of a compliant substrate with thickness of h_s subjected to a compressive load $\widehat{P}$ and non-conservative force $\widehat{q}$. Buckling occurs after loading because the thickness of the film is extremely small, which is very similar to a slender beam buckling. In the analytical model, the thin film is considered as Euler–Bernoulli beam with fixed condition at both ends, since the thickness of film h_f is far less than the length L [42]. Planes of the cross sections remain planes after deformation, and the plane of the cross section is still perpendicular to the axis after deformation. It is unnecessary to employ the von Karman plate theory to model the film because it leads to lengthy solutions not convenient for practical use due to the in-plane displacement and the shear traction [44]. A common approximation is to ignore the in-plane displacement and the shear traction [45,46]. The film deformed out of a plane only can be modeled into a beam, which greatly simplifies theoretical analysis [24]. The compliant substrate is considered to be a Winkler elastic foundation [47,48,49,50], whose reaction at any point is proportional to the deflection, with stiffness $\widehat{k}$ and deflection of the thin film $\widehat{w}$ shown in Figure 1b, where $\widehat{w}$ is the function of Cartesian coordinate $\widehat{x}$, which is along the axial direction of the thin film.

The kinetic energy of the thin film intra-domain is

$$T=\frac{1}{2}m{\displaystyle {\int }_0^L{(\frac{\partial \widehat{w}}{\partial \widehat{t}})}^2\mathrm{d}\widehat{x}},$$

(1)

where L is the undeformed length of the film, and m is the mass density. The total static potential energy in the film can be expressed as

$$V=\frac{1}{2}EI{\displaystyle {\int }_0^L{(\frac{{\partial }^2\widehat{w}}{\partial {\widehat{x}}^2})}^2}\mathrm{d}\widehat{x}-\frac{1}{2}\widehat{P}{\displaystyle {\int }_0^L{(\frac{\partial \widehat{w}}{\partial \widehat{x}})}^2}\mathrm{d}\widehat{x}+\frac{EA}{8L}{[{\displaystyle {\int }_0^L{(\frac{\partial \widehat{w}}{\partial \widehat{x}})}^2}\mathrm{d}\widehat{x}]}^2+\frac{1}{2}{\displaystyle {\int }_0^L\widehat{k}{\widehat{w}}^2}\mathrm{d}\widehat{x},$$

(2)

where E is Young’s modulus, and A and I are the area and the moment of inertia of the cross section, respectively. Here, the potential energy due to bending can be expressed as $\frac{1}{2}EI{\displaystyle {\int }_0^L{(\frac{{\partial }^2\widehat{w}}{\partial {\widehat{x}}^2})}^2}\mathrm{d}\widehat{x}$. The potential energy due to the axial force P is given by $-\frac{1}{2}\widehat{P}{\displaystyle {\int }_0^L{(\frac{\partial \widehat{w}}{\partial \widehat{x}})}^2}\mathrm{d}\widehat{x}$. The potential energy due to the midplane stretching is given by $\frac{EA}{8L}{[{\displaystyle {\int }_0^L{(\frac{\partial \widehat{w}}{\partial \widehat{x}})}^2}\mathrm{d}\widehat{x}]}^2$. The potential energy due to elastic foundation is given by $\frac{1}{2}{\displaystyle {\int }_0^L\widehat{k}{\widehat{w}}^2}\mathrm{d}\widehat{x}$. Then, the work done by non-conservative force in the film is

$$W_{nc}={\displaystyle {\int }_0^L\widehat{q}\widehat{w}\mathrm{d}\widehat{x}}$$

(3)

According to Hamilton’s principle ${\displaystyle {\int }_{{\widehat{t}}_0}^{{\widehat{t}}_f}(\delta T-\delta V+\delta W_{nc})\mathrm{d}\widehat{t}=0}$, it can be derived,

$$\begin{array}{l}{\displaystyle {\int }_{{\widehat{t}}_0}^{{\widehat{t}}_f}{\displaystyle {\int }_0^L[-EI\frac{{\partial }^4\widehat{w}}{\partial {\widehat{x}}^4}}-\widehat{P}\frac{{\partial }^2\widehat{w}}{\partial {\widehat{x}}^2}+\frac{EA}{2L}\frac{{\partial }^2\widehat{w}}{\partial {\widehat{x}}^2}{\displaystyle {\int }_0^L{(\frac{\partial \widehat{w}}{\partial \widehat{x}})}^2}\mathrm{d}\widehat{x}+\widehat{q}-\widehat{k}\widehat{w}-m\frac{{\partial }^2\widehat{w}}{\partial {\widehat{t}}^2}]\delta \widehat{w}\mathrm{d}\widehat{x}\mathrm{d}\widehat{t}}\\ -{\displaystyle {\int }_{{\widehat{t}}_0}^{{\widehat{t}}_f}{[EI\frac{{\partial }^2\widehat{w}}{\partial {\widehat{x}}^2}\delta (\frac{\partial \widehat{w}}{\partial \widehat{x}})]}_0^L\mathrm{d}\widehat{t}}+{\displaystyle {\int }_{{\widehat{t}}_0}^{{\widehat{t}}_f}{[(EI\frac{{\partial }^3\widehat{w}}{\partial {\widehat{x}}^3}+\widehat{P}\frac{\partial \widehat{w}}{\partial \widehat{x}}-\frac{EA}{2L}\frac{\partial \widehat{w}}{\partial \widehat{x}}{\displaystyle {\int }_0^L{(\frac{\partial \widehat{w}}{\partial \widehat{x}})}^2}\mathrm{d}\widehat{x})\delta \widehat{w}]}_0^L\mathrm{d}\widehat{t}}=0\end{array}$$

(4)

where $\widehat{t}$ is time. In consideration of the arbitrariness of variation and boundary conditions for fixed support at both ends, the governing equations are obtained:

$$m\frac{{\partial }^2\widehat{w}}{\partial {\widehat{t}}^2}+EI\frac{{\partial }^4\widehat{w}}{\partial {\widehat{x}}^4}+P\frac{{\partial }^2\widehat{w}}{\partial {\widehat{x}}^2}+\widehat{k}\widehat{w}-\frac{EA}{2L}\frac{{\partial }^2\widehat{w}}{\partial {\widehat{x}}^2}{\displaystyle {\int }_0^L{(\frac{\partial \widehat{w}}{\partial \widehat{x}})}^2}\mathrm{d}\widehat{x}-\widehat{q}=0$$

(5)

The boundary condition is $\partial \widehat{w}/\partial \widehat{x}=0$ and $\widehat{w}=0$ at $\widehat{x}=0$ and $\widehat{x}=L$, respectively. For the convenience of deduction, we use the non-dimensional variables as below

$$\begin{array}{l}x=\frac{\widehat{x}}{L},w=\frac{\widehat{w}}{r},t=\widehat{t}\sqrt{\frac{EI}{m{L}^4}},\\ q=\widehat{q}\frac{L^4}{rEI},P=\frac{\widehat{P}{L}^2}{EI},k=\widehat{k}\frac{L^4}{EI}\end{array}$$

(6)

where $r=\sqrt{I/A}$ is the radius of gyration of the cross section. Substituting Equation (6) into Equation (5) gives,

$$\ddot{w}+w^{\mathrm{IV}}+P{w}^{″}+kw-\frac{1}{2}{w}^{″}{\displaystyle {\int }_0^1{w^{\prime }}^2}\mathrm{d}x-q=0,$$

(7)

for fixed condition at both ends,

$$w^{\prime }=0\mathrm{and}w=0,\mathrm{at}(x=0,x=1),$$

(8)

where superscript “˙” and “′” stand for the derivation to t and x, respectively.

In order to solving buckling problem first, the buckling equilibrium equation can be obtained through omitting the time term and transverse load term from Equation (7) and denoting the buckled configuration by $\Phi (x)$,

$$\{\begin{array}{c}{\Phi }^{\mathrm{IV}}+P{\Phi }^{″}-\frac{1}{2}{\Phi }^{″}{\displaystyle {\int }_0^1{{\Phi }^{\prime }}^2}\mathrm{d}x+kw=0\\ {\Phi }^{\prime }=0,(x=0,x=1)\\ {\Phi }^{\prime }=0,(x=0,x=1)\end{array}$$

(9)

Here, assuming ${\sigma }^2=P-\kappa =P-\frac{1}{2}{\displaystyle {\int }_0^1{{\Phi }^{\prime }}^2}\mathrm{d}x$ is a constant, then Equation (9) becomes

$$\{\begin{array}{c}{\Phi }^{\mathrm{IV}}+{\sigma }^2{\Phi }^{″}+k\Phi =0\hfill \\ {\Phi }^{\prime }=0,(x=0,x=1)\hfill \\ \Phi =0,(x=0,x=1)\hfill \end{array}$$

(10)

When ${\sigma }^4>4k$, the general solution of Equation (10) can be assumed as

$$\Phi (x)=c_1\sin (k_{1}x)+c_2\cos (k_{1}x)+c_3\sin (k_{2}x)+c_4\cos (k_{2}x)$$

(11)

where $k_1=\sqrt{\left({\sigma }^2+\sqrt{{\sigma }^4-4k}\right)/2}$ and $k_2=\sqrt{\left({\sigma }^2-\sqrt{{\sigma }^4-4k}\right)/2}$; c_i (i = 1,2,3,4) is the coefficient to be determined by the boundary conditions. Substituting Equation (11) into the four boundary conditions of Equation (10) derives

$$Ac=0,$$

(12)

where,

$$A=\left[\begin{array}{cccc}0& 1& 0& 1\\ \sin (k_1)& \cos (k_1)& \sin (k_2)& \sin (k_2)\\ k_1\cos (k_1)& -k_1\sin (k_1)& k_2\cos (k_2)& -k_2\sin (k_2)\\ k_1& 0& k_2& 0\end{array}\right],c=\left[\begin{array}{c}{c}_1\\ c_2\\ c_3\\ c_4\end{array}\right]$$

(13)

Demanding that the determinant of the coefficient matrix equals zero, σ can be obtained via solving the following characteristic equation,

$$\begin{array}{c}{k}_1{k}_{2}cos{\left(k_1\right)}^2-2k_1{k}_{2}cos\left(k_1\right)cos\left(k_2\right)+k_1{k}_{2}sin{\left(k_1\right)}^2-{k_1}^{2}sin\left(k_1\right)sin\left(k_2\right)\hfill \\ -{k_2}^{2}sin\left(k_1\right)sin\left(k_2\right)+k_1{k}_{2}cos{\left(k_2\right)}^2+k_1{k}_{2}sin{\left(k_2\right)}^2=0\hfill \end{array}$$

(14)

It can be seen from the equation ${\sigma }^2=P-\kappa =P-\frac{1}{2}{\displaystyle {\int }_0^1{{\Phi }^{\prime }}^2}\mathrm{d}x$ that $P\ge {\sigma }^2$, which means that the critical buckling force is ${\sigma }^2$. We can obtain the value of c_i by solving the Equation (12),

$$\begin{array}{c}{c}_1=-c\frac{k_2}{k_1}\hfill \\ c_2=c\frac{k_{1}sin\left(k_2\right)-k_{2}sin\left(k_1\right)}{-k_{1}cos\left(k_1\right)+k_{1}cos\left(k_2\right)},\hfill \\ c_3=c\hfill \\ c_4=-c\frac{k_{1}sin\left(k_2\right)-k_{2}sin\left(k_1\right)}{-k_{1}cos\left(k_1\right)+k_{1}cos\left(k_2\right)}\hfill \end{array}$$

(15)

where c is constant. Substituting Equation (15) and Equation (11) into ${\sigma }^2=P-\kappa =P-\frac{1}{2}{\displaystyle {\int }_0^1{{\Phi }^{\prime }}^2}\mathrm{d}x$ gives the value of c, and the bucking mode $\Phi (x)$ can be obtained.

3. Vibration Analysis of Post-Buckled System

In order to investigate the vibration problem of the beam near the buckling configuration, we need to introduce a small dynamic displacement $\epsilon (x,t)$,

$$w(x,t)=\Phi (x)+\epsilon (x,t)$$

(16)

Substituting Equation (16) into Equation (7) gives,

$${\epsilon }^{\mathrm{IV}}+\ddot{\epsilon }+k\dot{\epsilon }+{\sigma }^2{\epsilon }^{″}=\frac{1}{2}{\Phi }^{″}{\int }_0^1(2{\Phi }^{\prime }{\epsilon }^{\prime }+{{\epsilon }^{\prime }}^2)\mathrm{d}x+\frac{1}{2}{\epsilon }^{″}{\int }_0^1(2{\Phi }^{\prime }{\epsilon }^{\prime }+{{\epsilon }^{\prime }}^2)\mathrm{d}x+q,$$

(17)

with boundary conditions: $\epsilon =0$ and ${\epsilon }^{\prime }=0$ at $x=0,x=1$. For the post-buckling linear free vibration problem, the post-buckling vibration equilibrium equation can be derived by omitting the nonlinear term and the transverse load term from Equation (17),

$$\ddot{\epsilon }+{\epsilon }^{\mathrm{IV}}+k\dot{\epsilon }+{\sigma }^2{\epsilon }^{″}={\Phi }^{″}{\int }_0^1{\Phi }^{\prime }{\epsilon }^{\prime }\mathrm{d}x$$

(18)

We assume that $\epsilon (x,t)$ has the following form,

$$\epsilon (x,t)=\xi (x)e^{i\omega t},$$

(19)

where $\xi (x)$ denotes the post-buckling vibrational mode, and ω denotes the vibration frequency. Substituting Equation (19) into Equation (18),

$$\begin{array}{c}{\xi }^{\mathrm{IV}}+{\sigma }^2{\xi }^{″}-{\omega }^2\xi +k\xi ={\Phi }^{″}{\int }_0^1{\Phi }^{Ȃ}{\xi }^{\prime }\mathrm{d}x\hfill \\ \xi =0(x=0,1)\hfill \\ {\xi }^{\prime }=0(x=0,1)\hfill \end{array}$$

(20)

The general solution of Equation (20) can be expressed as,

$$\xi (x)={\xi }_h(x)+{\xi }_p(x),$$

(21)

where ${\xi }_h(x)$ is the general solution of the equation ${\xi }^{\mathrm{IV}}+{\sigma }^2{\xi }^{″}-({\omega }^2-k)\xi =0$, which can be expressed as when ${\omega }^2>k$,

$${\xi }_h(x)=d_1\sin (s_{1}x)+d_2\cos (s_{1}x)+d_3\sinh (s_{2}x)+d_4\cosh (s_{2}x),$$

(22)

where d_i are constants, and s_i is

$$s_{1,2}={\left(\pm \frac{1}{2}{\sigma }^2+\frac{1}{2}\sqrt{{\sigma }^4+4({\omega }^2-k)}\right)}^{\frac{1}{2}}$$

(23)

The particular solution ${\xi }_p(x)$ can satisfy Equation (20),

$${{\xi }_p}^{\mathrm{IV}}+{\sigma }^2{{\xi }_p}^{″}-{\omega }^2{\xi }_p+k{\xi }_p={\Phi }^{″}{\int }_0^1{\Phi }^{\prime }{{\xi }_h}^{\prime }\mathrm{d}x+{\Phi }^{″}{\int }_0^1{\Phi }^{\prime }{{\xi }_p}^{\prime }\mathrm{d}x$$

(24)

Assuming that

$${\xi }_p(x)=d_5{\Phi }^{″}$$

(25)

Substituting Equation (25) into Equation (24) considering ${\Phi }^{\mathrm{IV}}+{\sigma }^2{\Phi }^{″}+k\Phi =0$ provides,

$${\int }_0^1{\Phi }^{\prime }{{\xi }_h}^{\prime }\mathbf{d}x+d_5({\omega }^2+{\int }_0^1{\Phi }^{\prime }{\Phi }^{‴}\mathbf{d}x)=0$$

(26)

As a result, the solution of Equation (20) can be expressed as,

$$\xi (x)=d_1\sin (s_{1}x)+d_2\cos (s_{1}x)+d_3\sinh (s_{2}x)+d_4\cosh (s_{2}x)+d_5{\Phi }^{″}$$

(27)

Then, substituting $\Phi (x)$ and $\xi (x)$ into Equation (26) and boundary conditions of Equation (20) gives homogeneous linear equations, which can be rewritten as,

$$\mathbf{B}\mathbf{d}=\mathbf{0},$$

(28)

where

$$B=\left[\begin{array}{ccccc}0& 1& 0& 1& -c_2{k}_1{}^2-c_4{k}_2{}^2\\ s_1& 0& s_2& 0& -c_1{k}_1{}^3-c_3{k}_2{}^3\\ \sin (s_1)& \cos (s_1)& \sinh (s_2)& \cosh (s_2)& {\alpha }_6\\ s_1\cos (s_1)& -s_1\sin (s_1)& s_2\cosh (s_2)& s_2\sinh (s_2)& {\alpha }_7\\ {\alpha }_1& {\alpha }_2& {\alpha }_3& {\alpha }_4& {\omega }^2+{\alpha }_5\end{array}\right],\mathbf{d}=\left[\begin{array}{c}{d}_1\\ d_2\\ d_3\\ d_4\\ d_5\end{array}\right]$$

(29)

where α_i are constants as below,

$$\begin{array}{c}{\alpha }_1={\int }_0^1{\Phi }^{\prime }{\sin }^{\prime }\left(s_{1}x\right)\mathrm{d}x\hfill \\ {\alpha }_2={\int }_0^1{\Phi }^{\prime }{\cos }^{\prime }\left(s_{1}x\right)\mathrm{d}x\hfill \\ {\alpha }_3={\int }_0^1{\Phi }^{\prime }{\sinh }^{\prime }\left(s_{2}x\right)\mathrm{d}x\hfill \\ {\alpha }_4={\int }_0^1{\Phi }^{\prime }{\cosh }^{\prime }\left(s_{2}x\right)\mathrm{d}x\hfill \\ {\alpha }_5={\int }_0^1{\Phi }^{\prime }{\Phi }^{‴}\mathrm{d}x\hfill \\ {\alpha }_6=-c_1{k_1}^{3}cos\left(k_1\right)+c_1{k_1}^{3}sin\left(k_1\right)-c_1{k_2}^{3}cos\left(k_2\right)+c_1{k_2}^{3}sin\left(k_2\right)\hfill \\ {\alpha }_7=-c_1{k_1}^{3}cos\left(k_1\right)+c_1{k_1}^{3}sin\left(k_1\right)-c_1{k_2}^{3}cos\left(k_2\right)+c_1{k_2}^{3}sin\left(k_2\right)\hfill \end{array}$$

(30)

Since d_i cannot simultaneously be zero, $\det (B)=0$, which provides the condition to obtain the value of ω. Furthermore, the vibrational mode $\xi (x)$ corresponding to natural frequency ω can be obtained.

4. Results and Discussion

The effective stiffness of substrate k can have significant effects on the buckling behaviors of the thin film. Figure 2 shows the deflection of the first-order buckling mode in the films with different substrate stiffness with normalized pre-stress of 10 × π². When the pre-stress is chosen as 0.001 very close to zero, the deflection can nearly equal the results obtained without substrate [29], which can verify the accuracy of the analytical model in the post-buckling analysis. As substrate stiffness increases, the Young’s modulus of compliant substrate increases, and the deflection of the thin films decreases. When substrate stiffness changes from 0 to 500, the deflection of the film decreases from 4.9 to 3.0.

The first-order critical buckling force and the buckling modes can also be affected by substrate stiffness, as shown in Figure 3. When k is lower than 877, the buckling mode is shown in Figure 3A, but with increase of k, the first-order buckling mode gradually varies from global buckling to local buckling. The values of k at the transition points between the first fourth buckling modes are 877, 6225, and 21,750, respectively. Figure 3D exhibits the first-order buckling mode with k of 5.5 × 10⁷ where the number of wrinkles reaches up to 27, which is the typical morphing of local buckling. For comparisons, the number of wrinkles is just two or three when k belongs to the range of 877~6225 or 6225~21,750, respectively. The buckling modes and the first-order critical buckling force obtained from theoretical calculation agree reasonably well with simulation results.

Figure 4 illustrates the deflections of first and second vibration modes of the post-buckled thin film in the first-order buckling mode with compressive load of 10 × π². When k is close to zero, the vibration mode can be validated by the results obtained by Ref. [29], and if k is larger, the deflections obviously decrease. When k increases from 0 to 500, the deflection of the first-order vibration mode decreases from 57.6 to 40.9.

The substrate stiffness can have negative effects on the natural frequencies, which can be found in Figure 5. Blue lines and red lines in the main figure stand for the first and the second order natural frequencies for the system. However, the corresponding buckling modes are shown in the graphs A–D, where A and B show the vibration mode in the first buckling mode; C and D show the vibration modes in the second buckling mode. With increase of k, the first linear natural frequency and the second linear natural frequency both decrease gradually, and the first vibration mode disappears when the first nature frequency decreases to 24.13 below $\sqrt{k}$ in first buckling mode and decreases to 60.71 below $\sqrt{k}$ in second buckling mode, because Equation (22) turns into the condition ${\omega }^2>k$.

A typical stretchable electronics example, mechanical response of copper thin film on PDMS substrate subjected to a compressive load, is used to demonstrate the application of theoretical calculation. The buckling of copper thin film on PDMS substrate is widely used in various flexible electronic devices as the bridges in “island-bridge” [51,52]. The elastic modulus E_s, the width b, the thickness h_s, and the length L of PDMS substrate are 2 MPa, 4 mm, 1 mm, and 10 mm, respectively. The elastic modulus E_f and the thickness h_f of copper thin film are 71,000 MPa and 0.01 mm. Meanwhile, the structure with fixed condition at both ends is subjected to a compressive load $\widehat{P}$ = 2.367 × 10⁻⁴ N. The dimensionless quantities k and P are obtained by Equation (6),

$$P=5000,k=3.38\times {10}^6$$

(31)

The mechanical response of this structure can be obtained by the theoretical analysis. The first order critical buckling force is 3677.2, and the first two orders natural frequencies in first buckling mode are 75.40 and 452.22. Figure 6 illustrates the first order buckling mode and its first two vibration modes. The wrinkle mode is helpful to understand the mechanical behavior of stretchable electronics [53,54].

5. Conclusions

In summary, an analytical model was developed for the dynamic behaviors of post-buckled thin films on a compliant substrate. The effects on buckling modes, buckling critical force, and vibration modes via stiffness of the substrate were investigated systematically. The increase of substrate stiffness can obviously reduce the deflection of post-buckling modes and vibration modes. Meanwhile, the first order buckling mode gradually varies from global buckling to local buckling, and the buckling mode is wrinkled. With the increase of substrate stiffness, the natural frequencies in the first two orders buckling modes drop noticeably. The results serve as guidelines for dynamic design of stretchable electronics to avoid resonance in a complicated noise environment, which is of great significance for the accurate measurement and the long-term use of flexible electronic devices [9,10,11,12,13,14] integrated on elastic substrate worked in complex service environments.