3.1 Effect of rigid elements on wrinkles patterns

The relationship between whether the rigid element inclines with a wrinkle and the size of the rigid element was analysed. For tension ${T_1} = {T_2}$ , an eigenvalue buckling analysis of the thin-film structure shown in Fig. 1 was conducted using ABAQUS finite-element analysis (FEA) software. The thin shell element S4R with four nodes and six degrees of freedom was used for the membrane structure, and the hexahedral element C3D8R with eight nodes and three degrees of freedom was used for the rigid element. All finite-element dimensions were 1 mm. The rigid element and membrane were connected by binding constraints. An eigenvalue buckling analysis can be used to determine the buckling mode of a thin film. The eigenvalue buckling results when the location of the rigid element ${R_c} = 30\;{\rm{mm}}$ are shown in Fig. 4.

Through eigenvalue buckling analysis, two buckling modes can be obtained. When rigid elements are present on the thin film, it can have two different modes of wrinkling. The rigid elements in mode 1 are inclined with the appearance of wrinkles, and the rigid elements in mode 2 are not inclined with the appearance of wrinkles. This result is consistent with the experimental results shown in Figs. 2 and 3. Therefore, the wrinkle mode of the nonuniform thin-film structure of rigid elements is mainly determined by the relationship between the radius of the circumscribed circle of the rigid element and the wrinkle wavelength. In other words, when the wrinkle half-wavelength $d \lt 2{\lambda _1}$ , the wrinkle mode of the film is consistent with mode 1; when the wrinkle half-wavelength $d \geqslant 2{\lambda _2}$ , the wrinkle mode of the thin film is consistent with the mode 2. The wavelength ${\rm{\lambda }}$ of the thin-film structure can be expressed as [Reference Bonin16]

(1) \begin{align}\lambda = {\left[ {\frac{{E{\pi ^2}{t^2}{R_w}^2}}{{12\!\left( {1 - {\nu ^2}} \right)\!{\sigma _\rho }}}} \right]^{\frac{1}{4}}}\end{align}

where E is the elastic modulus of the thin film, v is Poisson’s ratio of the thin film, t is the thickness of the thin film, and ${\sigma _\rho }$ is the radial stress of the nonuniform thin film. At the thin-film boundary, the influence of the disturbance stress of the rigid element on the thin-film stress field can be ignored. At this point, the radial stress of the nonuniform thin film is consistent with the radial stress of the uniform thin film, which can be expressed as

(2) \begin{align}{\sigma _\rho } = \frac{{4T\cos \theta }}{{\!\left( {\pi + 2} \right)\!\rho t}}\end{align}

Table 1. Material parameters of thin film and rigid elements

Figure 2. Wrinkling patterns of rigid element at different locations under uniform tension: (a)–(d) show positions R_c of the rigid elements of 30, 40, 50, and 60mm, respectively.

Figure 3. Wrinkling patterns of rigid element at different locations under nonuniform tension: (a)–(d) show positions R_c of the rigid elements of 30, 40, 50, and 60mm, respectively.

Figure 4. Eigenvalue buckling analysis of thin-film structures with rigid elements.

Here, ${R_w}$ is the radius of the wrinkled region in the thin film. The boundary region of the thin film can be obtained [Reference Bonin16] as

(3) \begin{align}{R_w} = {e^{ - v}}R\end{align}

According to the experimental observations, when ${T_1}/{T_2} = 4$ , wrinkles occur across the diagonal region of the thin film. Therefore, the tension ratio is defined as ${\rm{c}} = {T_1}/{T_2}$ . When the wrinkle pattern is shown in mode 2, the global wrinkle radius ${R_{w1}}$ and local wrinkle radius ${R_{w2}}$ can be expressed as

(4) \begin{align}\left\{ \begin{array}{l@{\quad}l} {{R_{w1}} = {{{R_c}} / {\sin \alpha }}} \\[5pt] {{R_{w2}} = {R_w} - {R_{w1}}} \end{array} \right.\end{align}

Here, ${{\rm{\lambda }}_1}$ is the first half-wavelength close to the central line and can be expressed as

(5) \begin{align}{\lambda _1} = \left( {\frac{{{R_c}}}{{\sin \alpha }} + \frac{d}{2}} \right)\!\tan \!\left( {\frac{\alpha }{n}} \right)\end{align}

The wrinkling number n of the thin film can be expressed as

(6) \begin{align}n = \frac{{\alpha {R_w}}}{\lambda }\end{align}

At this time, one can obtain the relationship between the side length and the wrinkle mode when the position of the rigid element ${R_c} = 30\;{\rm{mm}}$ , as shown in Fig. 5. Figure 3 reveals that, when the position of the rigid element ${R_c} = 30\;{\rm{mm}}$ , the side length of the rigid element $2a = 10\;{\rm{mm}}$ , and ${T_1} = {T_2}$ , the wrinkle mode of the nonuniform thin film is mode 2 – that is, the rigid element does not tilt with the appearance of wrinkles.

Figure 5. Relationship between side length of rigid element and wrinkled mode.

3.2 Wrinkle model when rigid element is inclined with the wrinkle

The experimental results in Fig. 2 show that, when the diameter of the circumscribed circle of the rigid element is $\;d \geqslant 2{\lambda _1}$ , the wrinkle mode of the nonuniform thin film is mode 2. To describe the wrinkle pattern in mode 2, the wrinkle configuration function w on a nonuniform thin film can be expressed as

(7) \begin{align}w = \sum\limits_{i = 1}^4 {\left( {{w_{i1}} + {w_{i2}}} \right)} \end{align}

According to thin-film stability theory, the local wrinkle configuration function ${w_{i1}}$ and the global wrinkle configuration function ${w_{i2}}$ can be expressed as

(8) \begin{align}\begin{array}{l}{w_{i1}} = {A_{i1}}\sin\!\left( {\dfrac{{\pi \rho }}{{{c^{{{\left( { - 1} \right)}^i}}}{R_{w1}}}}} \right)\!\sin\!\left( {\dfrac{{\pi \theta }}{\lambda }} \right)\\[10pt] {w_{i2}} = {A_{i2}}\sin\!\left[ {\dfrac{{\pi \!\left( {\rho - l + \dfrac{{{R_{w2}}}}{2}} \right)}}{{{R_{w2}}}}} \right]\!\sin\!\left[ {\dfrac{{\pi \left( {\theta - \beta } \right)}}{\lambda }} \right]\end{array}\end{align}

The geometric strain ${\varepsilon _{\theta w}}$ in the wrinkle direction can be expressed as

(9) \begin{align}{\varepsilon _{\theta w}} = {\varepsilon _{\theta c}} - {\varepsilon _{\theta \rho }} = - \frac{{{A^2}{\pi ^2}}}{{4{\lambda ^2}}}\end{align}

where ${\varepsilon _{\theta \rho }}$ and ${\varepsilon _{\theta c}}$ are the tensile strain in the radial direction and the shrinkage strain in the circumferential direction, respectively, and can be expressed as

(10) \begin{align}{\varepsilon _{\theta \rho }} = - \nu \frac{{{\sigma _\rho }}}{E}{\rm{ }}\end{align}

(11) \begin{align}{\varepsilon _{\theta c}} = \frac{{\int {\frac{{{\sigma _\rho }}}{E}dy} }}{y}\end{align}

Equations (10) and (11) can be substituted into Equation (9) to obtain the wrinkle amplitude $A\!\left( {\rho, \theta } \right)$ of the thin film:

(12) \begin{align}A\!\left( {\rho, \theta } \right) = \frac{{2\lambda \sqrt {\left( {{\varepsilon _{\theta \rho }} - {\varepsilon _{\theta c}}} \right)} }}{\pi }\end{align}

The strain field of a nonuniform thin film can be calculated using Eshelby’s elastic inclusion theory. According to previous research [Reference Li, Niu, Wu, Zhang, Luo and Zhang19], the stress distribution in thin films with rigid inclusions can be expressed as

(13) \begin{align}{\sigma _{ij}} = \sigma _{ij}^0 + {\sigma '_{\!\!ij}} + \sigma _{ij}^c\end{align}

where $\sigma _{ij}^0$ is the stress of the uniform thin film, $\sigma _{ij}{\prime}$ is the disturbance stress caused by the rigid element, and $\sigma _{ij}^c$ is the concentrated stress at the corners of the rigid element, which are respectively expressed as [Reference Sun, Huang, Zhang and Meng22]

(14) \begin{align}{\sigma '_{\!\!ij}} = {C_{ijkl}}\!\left( {{D_{ijkl}}\varepsilon _{ij}^ * } \right)\end{align}

(15) \begin{align}\sigma _{ij}^c = \frac{{{a^2}}}{{\pi \!\left( {{r^2} - {a^2}} \right)}}\frac{{\left| {\left| x \right| - \left| y \right|} \right|}}{r}{\sigma '_{\!\!ij}}\end{align}

(16) \begin{align}\sigma _{ij}^0 = \left\{ {\begin{array}{l}{\sigma _x^0 = - \dfrac{{4T}}{{2\pi R}} + \sum\limits_{i = 1}^{i = 4} {\dfrac{{4T{{\left( {x\cos i\theta - y\sin i\theta } \right)}^2}\left[ {R - \left( {x\sin i\theta + y\cos i\theta } \right)} \right]}}{{{\rho ^4}{h_m}(\pi + 2)}}} }\\[12pt] {\sigma _y^0 = - \dfrac{{4T}}{{2\pi R}} + \sum\limits_{i = 1}^{i = 4} {\dfrac{{4T{{\left[ {R - \left( {x\sin i\theta + y\cos i\theta } \right)} \right]}^3}}}{{{\rho ^4}{h_m}(\pi + 2)}}{\rm{ }}} }\\[12pt] {\tau _{xy}^0 = - \dfrac{{4T}}{{2\pi R}} + \sum\limits_{i = 1}^{i = 4} {\dfrac{{4T\left( {x\cos i\theta - y\sin i\theta } \right){{\left[ {R - \left( {x\sin i\theta + y\cos i\theta } \right)} \right]}^2}}}{{{\rho ^4}{h_m}(\pi + 2)}}} }\end{array}} \right.\end{align}

Here, ${C_{ijkl}}$ is the elastic matrix of the film, $\varepsilon _{ij}^*$ is the characteristic strain of the film, and ${D_{ijkl}}$ is the Eshelby external tensor, which can be expressed as [Reference Sun, Huang, Zhang and Meng22]

(17) \begin{align} {D_{1111}}\!\left( r \right) &= \frac{1}{{8\!\left( {1 - \nu } \right)}}\left[ {\frac{{3{a^4}}}{{{r^4}}} + \left( {1 - 2\nu } \right)\frac{{2{a^2}}}{{{r^2}}}} \right] + \frac{1}{{\left( {1 - \nu } \right)}}\left[ {\left( {1 + \nu } \right)\frac{{{a^2}}}{{{r^4}}}x_n^2 - \frac{{3{a^4}}}{{{r^6}}}x_n^2 - \frac{{2{a^2}}}{{{r^6}}}x_n^4 + \frac{{3{a^4}}}{{{r^8}}}x_n^4} \right]\nonumber \\[5pt] {D_{1122}}\!\left( r \right) &= \frac{1}{{8\!\left( {1 - \nu } \right)}}\left[ {\frac{{{a^4}}}{{{r^4}}} - \left( {1 - 2\nu } \right)\frac{{2{a^2}}}{{{r^2}}}} \right]\nonumber \\[5pt] & \quad + \frac{1}{{2\!\left( {1 - \nu } \right)}}\left\{ {\left[ {\left( {1 - 2\nu } \right)x_n^2 + y_n^2} \right]\frac{{{a^2}}}{{{r^4}}} - \left( {y_n^2 + \frac{1}{{{a^2}}}x_n^2y_n^2} \right)\frac{{{a^4}}}{{{r^6}}} + \frac{{6{a^4}}}{{{r^8}}}x_n^2y_n^2} \right\}\nonumber\\[5pt] {D_{2211}}\!\left( r \right) &= \frac{1}{{8\!\left( {1 - \nu } \right)}}\left[ {\frac{{{a^4}}}{{{r^4}}} - \left( {1 - 2\nu } \right)\frac{{2{a^2}}}{{{r^2}}}} \right]\nonumber \\[5pt] & \quad + \frac{1}{{2\!\left( {1 - \nu } \right)}}\left\{ {\left[ {\left( {1 - 2\nu } \right)y_n^2 + x_n^2} \right]\frac{{{a^2}}}{{{r^4}}} - \left[ {x_n^2 + \frac{1}{{{a^2}}}x_n^2y_n^2} \right]\frac{{{a^4}}}{{{r^6}}} + \frac{{6{a^4}}}{{{r^8}}}x_n^2y_n^2} \right\}\nonumber\\[5pt] {D_{2222}}\!\left( r \right) &= \frac{1}{{8\!\left( {1 - \nu } \right)}}\left[ {\frac{{3{a^4}}}{{{r^4}}} + \left( {1 - 2\nu } \right)\frac{{2{a^2}}}{{{r^2}}}} \right] + \frac{1}{{\left( {1 - \nu } \right)}}\left[ {\left( {1 + \nu } \right)\frac{{{a^2}}}{{{r^4}}}y_n^2 - \frac{{3{a^4}}}{{{r^6}}}y_n^2 - \frac{{2{a^2}}}{{{r^6}}}y_n^2 + \frac{{3{a^4}}}{{{r^8}}}y_n^2} \right]\nonumber\\[5pt] {D_{1212}}\!\left( r \right) &= \frac{1}{{2\!\left( {1 - \nu } \right)}}\left[ {\frac{{{a^4}}}{{{r^4}}} + \left( {1 - 2\nu } \right)\frac{{2{a^2}}}{{{r^2}}}} \right] + \frac{1}{{1 - \nu }}\left[ {\nu\!\left( {x_n^2 + y_n^2} \right)\frac{{{a^2}}}{{{r^4}}} - \left( {x_n^2 + y_n^2} \right)\frac{{{a^4}}}{{{r^6}}}} \right] \end{align}

When $\rho = {R_{w1}}/2$ and $\theta = 0$ , an exponential function obtained through finite-element analysis is introduced to modify the wrinkle amplitude, and the maximum amplitude ${A_{i1}}$ of the global wrinkle configuration function w is obtained as

(18) \begin{align}{A_{i1}} = {c^{{{\left( { - 1} \right)}^i}}}{e^{ - {{\left( {3t\left| {{x_i}} \right|} \right)}^2}}}A\!\left( {\frac{{{c^{{{\left( { - 1} \right)}^i}}}{R_{w1}}}}{2},0} \right)\end{align}

Because the rigid element on the thin film changes the original stress field of the thin film, local wrinkles caused by the disturbance stress of the rigid element appear on both sides of the thin film, and the amplitude of the wrinkles is the largest at the apex of the rigid element. When $\rho = l$ and ${\rm{\;}}\rho = l$ , the maximum amplitude ${A_{i2}}$ of the local wrinkle configuration function w caused by the rigid element is

(19) \begin{align}{A_{i2}} = {c^{{{\left( { - 1} \right)}^i}}}{e^{ - {{\left[ {3t\left(\left| {{x_i}} \right| - \frac{{d + {\lambda _2}}}{2}\right)} \right]}^2}}}A\!\left( {l,\beta } \right)\end{align}

where ${\lambda _2}$ is the half-wavelength of the first wrinkle outside the rigid element, which can be expressed as

(20) \begin{align}{\lambda _2} = \left( {\frac{{{R_c}}}{{\sin \alpha }} + \frac{d}{2}} \right)\!\tan \!\left( {\frac{\alpha }{n} + \beta } \right) - \frac{d}{2}\end{align}

By substituting Equations (18) and (19) into Equation (8), the wrinkle configuration function of the nonuniform thin-film structure of rigid elements corresponding to mode 1 can be obtained. When the rigid element is large, the wrinkle shape shown in mode 1 mostly appears in the case of uniform tension. Therefore, taking uniform stretching as an example, when the distance between the rigid elements ${R_c} = 30{\rm{mm}}$ and tension ${T_1} = {T_2} = 30{\rm{N}}$ , the wrinkle distribution diagrams obtained through the stretching experiment, finite-element simulation, and theoretical analysis are as shown in Fig. 6. The out-of-plane deformation in the path 1 cross section is also shown in Fig. 6.

Figure 6. Wrinkle shapes and out-of-plane deformation in the path cross section: (a) is the experimental result, (b) is the FEA result, (c) is the analytical result, and (d) is the out-of-plane deformation in the path 1 cross section.

3.3 Wrinkle model when the rigid element does not incline with the wrinkle

When the circumscribed circle diameter of the rigid element $d \lt 2{\lambda _1}$ , the wrinkle pattern of the nonuniform thin film is mode 1, and the rigid element tilts with the appearance of the wrinkle. Therefore, to describe the wrinkle morphology shown in mode 1, the wrinkle configuration function w of the nonuniform film can be expressed as

(21) \begin{align}w = \sum\limits_{i = 1}^4 {\left( {{w_{i1}} + {w_{i2}}} \right)} \end{align}

where ${w_{i1}}$ is the local wrinkle configuration caused by rigid elements, and ${w_{i2}}$ is the global wrinkle configuration. According to thin-film stability theory, wrinkles ${w_{i1}}$ and ${w_{i2}}$ can be expressed as

(22) \begin{align} {w_{i1}} &= {A_{i1}}\sin\!\left( {\frac{{\pi \rho }}{{{c^{{{\left( { - 1} \right)}^i}}}{R_w}}}} \right)\!\cos\! \left( {\frac{{\pi \theta }}{\lambda }} \right)\nonumber \\[5pt] {w_{i2}} &= {A_{i2}}\sin\!\left[ {\frac{{\pi \!\left( {\rho - l} \right)}}{{{R_{w2}}}}} \right]\!\cos\! \left[ {\frac{{\pi \!\left( {\theta - \beta } \right)}}{\lambda }} \right] \end{align}

When $\rho = {c^{{{\left( { - 1} \right)}^i}}}/2$ and $\theta = 0$ , the maximum amplitude A of the local wrinkle configuration function w caused by the rigid element is

(23) \begin{align}{A_{i1}} = {c^{{{\left( { - 1} \right)}^i}}}{e^{ - {{\left[ {3t\left( {\left| {{x_i}} \right| - \frac{{{\lambda _2}}}{2}} \right)} \right]}^2}}}A\!\left( {\frac{{{c^{{{\left( { - 1} \right)}^i}}}{R_{w1}}}}{2},0} \right)\end{align}

When $\rho = l$ and $\theta = \beta $ , the maximum amplitude A of the local wrinkle configuration function w caused by the rigid element is

(24) \begin{align}{A_{i2}} = {e^{ - {{\left[ {3t\left(\left| {{x_i}} \right| - \frac{{2{\lambda _1} + {\lambda _2}}}{2}\right)} \right]}^2}}}A\!\left( {l,\beta } \right)\end{align}

When $i$ is odd, ${x_i} = x$ . When $i$ is even, ${\rm{\;}}{x_i} = y$ . By substituting Equations (23) and (24) into Equation (22), the wrinkle configuration function of the nonuniform thin-film structure of rigid elements corresponding to mode 2 can be obtained. The slope of the rigid element can be obtained by substituting the coordinates of the endpoint of the rigid element into Equation (22). The wrinkle pattern shown in mode 1 is prone to occur in the case of nonuniform stretching or small rigid elements. Therefore, taking nonuniform stretching as an example, when the rigid element distance ${R_c} = 30$ and the tensile forces ${T_1} = 5{\rm{N}}$ and ${T_2} = 20{\rm{N}}$ , the wrinkle distribution diagrams obtained by the experiment, finite-element simulation, and theoretical analysis are as shown in Fig. 7. The out-of-plane deformations in the path 1 and 2 cross sections are shown in Fig. 7.

Figure 7. Wrinkle shapes and out-of-plane deformation in the path cross section: (a) is the experimental result, (b) is the FEA result, (c) is the analytical result, (d) is the out-of-plane deformation in the path 1 cross section, and (e) is the out-of-plane deformation in the path 2 cross section.