Buckling and Wrinkling of a Thin Solid Film with Non-uniform Thickness

Masoud Noroozi

Faculty of Mechatronics, Islamic Azad University (Karaj Branch), Karaj, Iran

E-mail: masoud.noroozi@kiau.ac.ir

Received 09 October 2020; Accepted 29 November 2020; Publication 09 January 2021

Abstract

The instability of a strip as a free-standing film and also a deposited film on a substrate is studied in this work. The non-uniform thickness of the film is assumed with a quadratic profile. The problem is categorized under the topic of structural stability and the eigenvalue problem corresponding with the ODE of the system is solved. For the free-standing film, the buckling loads and mode shapes are derived analytically through a closed-form solution. For the substrate-bonded film with a finite length, the uniaxial wrinkling of the film is investigated by using a series solution and a finite difference method and the wrinkling load and wrinkling pattern are characterized. Unlike the wrinkling of thin films with uniform thickness in which the wrinkles propagate along the entire span, it is shown that for the non-uniform film wrinkles are localized near the location with a minimum thickness along the length span; and the wrinkling accumulation is very sensitive to the thickness variations. Therefore, this work is expected to increase the insight into the localization of the wrinkles in thin film-substrate systems in engineering, industry and medical science.

Keywords: Thin solid film, structural stability, Winkler foundation, finite difference method.

1 Introduction

The mechanical instability of the beam/film structures categorized as buckling and wrinkling has attracted great attention of many researchers in sheet metal forming [1, 2], aeronautical systems and light-weight gossamer structures [3], semiconductor technology [4], stretchable electronics [5], woven composite plies [6], skin aging and wound healing [7], material properties metrology [8] and so on [9]. Amongst different patterns for wrinkling of a film on a substrate [4, 10], the uniaxial wrinkling is one of the common patterns in which wrinkles propagate uniaxially along the film length span like a strip [8, 11, 12]. Wrinkles may appear on the free-standing films under tension as studied analytically [13], experimentally and numerically [14]. On the other hand, thin films deposited on the substrate may undergo wrinkling due to compressive loading [15], thermal loading [16] and shrinking of the substrate [4].

The instability problems (i.e. buckling and wrinkling) are considered by many researchers using bifurcation theory [17] through the eigenvalue problem of the system. Birman and Bert [11] considered the wrinkling of the orthotropic films, Coman and Bassom [18] studied the wrinkling in the films under shear loading, Damil and Potier-Ferry [19] investigated the coupling between local and global instabilities of thin films using continuum models as a beam on an elastic foundation and a 3D nonlinear elasticity problem, and Chen and Hutchinson [4] determined the wrinkling parameters through the principle of minimum potential energy of the system. However, in all of these works, the effect of the non-uniformity of the film on the wrinkling problem is neglected. In other words, most of the existing theoretical works on the substrate-bonded films used the uniform assumption for the film thickness with a uniform wrinkling pattern along the entire span, in which the wrinkles propagate with a sinusoidal pattern all over the film with a uniform amplitude. Contrarily, for thin-film structures with finite length, the effect of the boundary conditions at the edges of the film is important on the pattern of the wrinkling and violates the uniformity of the wrinkling along the span [19].

Furthermore, the assumption of the uniform thickness of thin films is unreliable due to tiny thickness of the films, so the importance of the thickness variations and the effects of these changes on the instability parameters (i.e. load and pattern) need to be considered carefully. On the other hand, if the thickness can be varied/engineered in a controlled manner, then its effect on different desired buckling and wrinkling patterns and also other mechanical behaviors of the system seems to be very interesting in developing new ideas on various applications of thin-film structures. In this regard, one may refer to variable thickness superconducting thin films [20], optical systems and mirrors [21–23], sensors and detectors [24, 25], micro electro- mechanical systems [26, 27], etc. Besides, it is shown that the physical properties such as coefficient of thermal expansion [28] are thickness-dependent in thin-film structures, which highlights the importance of the thickness variations of thin films. The study of the non-uniformity in material properties of the film [29] and also the substrate [30] shows that the wrinkling pattern is very sensitive to non- homogeneity and the pattern is completely different from the results of the uniform (i.e. homogeneous) system. The effect of the non-uniform thickness of the film is investigated in this work for finite length film-substrate systems to discover more about the geometrical non-uniformity of the film.

In this study, the buckling and uniaxial wrinkling problem of a variable thickness film with finite length are investigated. The thickness profile is modeled with a quadratic pattern and the eigenvalue problem of the differential equation of the system is solved. In the buckling problem of the free-standing film, an analytical closed-form solution is proposed for the buckling loads and mode shapes. For the substrate-bonded film, the wrinkling load and wrinkling pattern are determined numerically by using a finite difference method and the results are compared with those from a series solution. In contrast with the other studies with uniform thickness assumption, the current work shows that the thickness non-uniformity along the domain is very effective in wrinkling localization, such that wrinkles accumulate around the thinnest location of the film. The tiny thickness of the film increases the significance of this study in thin film technology where the uniformity of the film thickness is not perfectly controllable. The results of the analysis are expected to increase the insight into the physics of the instability problem of thin-film structures and to provide more applications in science and technology of thin film systems.

2 Modeling

In order to characterize the buckling of a free-standing film and the uniaxial wrinkling of a substrate-bonded thin film, the instability problem of the system is studied in this work. The classical beam/strip theory with small deformation is used to model a rectangular thin film with thickness t, width b and length L (Figure 1). The substrate is modeled by using a Winkler foundation, in which the interaction between the film and the substrate is modeled by using a spring system. In other words, the effect of the substrate is applied to the film as an external load corresponding with the stiffness of the substrate and the deflection of the film. The Winkler modulus $\overline{K}$ depends on the substrate stiffness (i.e. Young Modulus E${}_s$) and its characteristic depth H as $\overline{K}=E_s/H$ [8, 11, 12]. For the system under uniaxial in-plane loading ${\overline{N}}_x$ in Figure 1, a uniaxial deformation strip-like pattern is developed in the film [10, 11]. The governing equation of this system is derived as [31]

$$\frac{d^2}{d{x}^2}\left[D\frac{d^2\overline{w}}{d{x}^2}\right]+{\overline{N}}_x\frac{d^2\overline{w}}{d{x}^2}+b\overline{K}\overline{w}=0$$

(1)

in which $\overline{w}$ is the deflection and D is the bending rigidity of the film defined by

$$D=\frac{1}{12}\overline{E_f}b{t}^3$$

(2)

where $\overline{E_f}=E_f/(1-v_f^2)$, E${}_f$ is Young’s modulus and v${}_f$ is the Poisson’s ratio of the film. In the current work, the thickness is assumed to vary along the length span as

$$t=t_0(1+\epsilon f(x))$$

(3)

and t${}_{\mathrm{0}}$ is the minimum thickness of the film, $\epsilon$ is thickness gradient denoting the amplitude of the variation of the film thickness, and f(x) is the profile of the thickness variation along the length span x such that Max $[\mathit{\text{f}}(\mathit{\text{x}})]=1$. By substituting Equation (3) into Equation (2), the bending rigidity is derived as

$$D=D_{0}D(x)=D_0{\left(1+\epsilon f(x)\right)}^3$$

(4)

where $D_0=\frac{1}{12}\overline{E_f}bt_0{}^3$. Imposing Equation (4) into Equation (1) and introducing the dimensionless variable $\xi$ defined by $\xi =\mathit{\text{x/L}}$ and normalized deflection $w=\overline{w}/{\overline{w}}_{max}$, the dimensionless governing equation is derived as

$$\left[D(\xi )\right]\frac{d^{4}w}{d{\xi }^4}+\left[2\frac{d}{d\xi }D(\xi )\right]\frac{d^{3}w}{d{\xi }^3}+\left[\frac{d^2}{d{\xi }^2}D(\xi )+N\right]\frac{d^{2}w}{d{\xi }^2}+Kw=0$$

(5)

in which $D(\xi )={(1+\epsilon f(\xi ))}^3$,

$N$	$={\displaystyle \frac{{\overline{N}}_x}{D_0}}{L}^2$	(6a)
$K$	$={\displaystyle \frac{\overline{K}b{L}^4}{D_0}}$	(6b)

In this work, the thickness of the film is assumed with a quadratic profile as illustrated in Figure 2 in which the film takes its minimum thickness at the middle of the length span as

$$f(\xi )=4{\xi }^2$$

(7)

and $-0.5\le \xi \le 0.5$. For this profile, the effect of the thickness variation of the film is investigated on the buckling and wrinkling and the results are compared with those of a film with uniform thickness.

Figure 1 Deposited film with variable thickness on the substrate.

Figure 2 The thickness profile of the film along the length span.

3 Solution Approaches

In this section the buckling problem of a free-standing film and also the uniaxial wrinkling of a film bonded to the substrate are studied. The buckling problem is solved analytically by direct solution of the ODE problem using well-known functions, where the analytical solution is available. For wrinkling of a film-substrate system, numerical methods are used to analyse the problem in the lack of analytical solution.

3.1 Buckling Analysis of a Free-standing Film

For a free-standing film with variable thickness, the governing equation of the system at the onset of instability is represented by

$$\frac{d^2}{d{\xi }^2}\left[D(\xi )\frac{d^{2}w}{d{\xi }^2}\right]+N\frac{d^{2}w}{d{\xi }^2}=0$$

(8)

where $\mathrm{D}(\xi )={(1+4\epsilon {\xi }^2)}^3$ is the bending rigidity of the film with quadratic thickness profile in Figure 2. For a film with non-uniform thickness (i.e. $\epsilon \ne 0$), by using a change of variable as $\mathrm{u}=1/{\sqrt{1+\epsilon \xi }}^2$ and some mathematical treatments of Equation (8), the algebraic form of Mathieu differential equation [32] is obtained as

$(1-u^2){\displaystyle \frac{d^{2}w}{d{u}^2}}-u{\displaystyle \frac{\mathrm{d}w}{\mathrm{d}u}}+[a+2q(1-2u^2)]w=0$	(9a)

which is converted to the canonical form of Mathieu differential equation by substituting $u=cos(\theta )$ as

${\displaystyle \frac{d^{2}w}{d{\theta }^2}}+[a-2q\cos (2\theta )]w=0$	(9b)

and MathieuC ($a$, $q$, $\theta$) and MathieuS ($a$, $q$, $\theta$) are even and odd Mathieu functions respectively, which are the solutions of the Mathieu differential Equations (9a, b). Furthermore, MathieuC and MathieuS are even and odd functions of $\theta$, respectively, and they are normalized such that at $\theta =0$, MathieuC ($a$, $q$, $\theta$) $=$ 1 and $d/d\theta$ [MathieuS ($a$, $q$, $\theta$)] $=$ 1 [32].

Therefore, the general solution of the governing equation of the film in Equation (8) is derived for $\epsilon \ne 0$ as

$$w(\xi )=m_1\mathrm{C}(\xi )+m_2\mathrm{S}(\xi )+m_3\xi +m_4$$

(10)

where $m_i(i=1\mathrm{\dots }4)$ are unknown constants and C($\xi$) and S($\xi$) is even and odd functions of $\xi$ as

$\mathrm{C}(\xi )$	$=\sqrt{1+\epsilon {\xi }^2}\mathrm{MathieuC}[1+{\displaystyle \frac{N}{2\epsilon }},-{\displaystyle \frac{N}{4\epsilon }},\arctan (\xi \sqrt{\epsilon })]$	(11a)
$\mathrm{S}(\xi )$	$=\sqrt{1+\epsilon {\xi }^2}\mathrm{MathieuS}[1+{\displaystyle \frac{N}{2\epsilon }},-{\displaystyle \frac{N}{4\epsilon }},\arctan (\xi \sqrt{\epsilon })]$	(11b)

In order to find the unknown constants $m_i(i=1\mathrm{\dots }4)$ in Equation (10), boundary conditions of the film are imposed. The clamped-clamped boundary conditions at the edges are applied with zero deflection and zero slope as $w=0$ and $\frac{dw}{d\xi }=0$. By imposing the boundary conditions into Equation (10) and after some mathematical treatments, the characteristic equations and mode shapes of the buckling of the film are obtained from the corresponding eigenvalue problem.

On the other hand, for a free-standing film with uniform thickness (i.e. $\epsilon =0$), the solution of the governing Equation (8) is derived by replacing S($\xi$) and C($\xi$) by $\sin (\xi \sqrt{N})$ and $\cos (\xi \sqrt{N})$ in Equation (10), respectively. The characteristic equation of the buckling of the film with clamped-clamped edges is represented by

$$\sin (0.5\sqrt{N})\left[\tan \left(0.5\sqrt{N}\right)-0.5\sqrt{N}\right]=0$$

(12)

which leads to the symmetric and antisymmetric buckling modes with buckling loads and mode shapes as follows:

• Symmetric mode:

$$N=4{\pi }^2\mathit{\hspace{1em}}\mathrm{and}\mathit{\hspace{1em}}w=1+\cos (\sqrt{N}\xi )$$

(13a)

• Antisymmetric mode:

$$N=8.183{\pi }^2\mathit{\hspace{1em}}\mathrm{and}\mathit{\hspace{1em}}w=\sin (\sqrt{N}\xi )-2\xi \sin (0.5\sqrt{N})$$

(13b)

3.2 Wrinkling Analysis of a Substrate Bonded Film

The governing Equation (5) of the system is a fourth-order linear ordinary differential equation with variable coefficients. These types of equations generally do not have a closed-form analytical solution. Therefore, we use semi-analytical solutions as well as numerical solutions to characterize the wrinkling behavior of the film-substrate system with non-uniform thickness. In this work, a series solution and a finite difference approach are used to represent the corresponding eigenvalue problem of the differential Equation (5) in the algebraic form. And by finding the eigenvalues and eigenvectors, one may find the wrinkling loads and wrinkling patterns.

The series solution of the film deflection w($\xi$) is constructed with unknown coefficients in the form of

$$w(\xi )=\sum _{i=0}^{i\to \mathrm{\infty }}{c}_i{\xi }^i$$

(14)

and by plugging w($\xi$) from Equation (14) into Equation (5) and rearranging terms, one may find the corresponding recurrence relation of the differential equation. In addition, for clamped-clamped edges of the film, the boundary conditions are represented as $w=0$ and $\mathit{\text{dw/d}}\xi =0$ at $\xi =\pm 0.5$. Imposing these conditions, the eigenvalue problem of the system is derived and the wrinkling load and wrinkling pattern are determined.

On the other hand, a finite difference method is used to solve the Equation (5) numerically [33] as introduced in Appendix. By discretizing the domain using nodes and adopting a central finite difference method, the differential equation is replaced with a set of algebraic equations in the following format,

$$[A]\{w\}+N[B]\{w\}=0$$

(15)

in which $\{\mathit{\text{w}}\}$ is the vector of the nodal displacement and [A] and [B] are square matrices. For the clamped-clamped boundary conditions at the edges of the film, a symmetric extension of the domain is adopted as explained by Zhao and Wei [34, 35]. The general eigenvalue problem (15) with eigenvector $\{$w$\mathrm{\}}$ and eigenvalue parameter N has a straight forward solution. The eigenvalues correspond with the wrinkling loads and the eigenvectors represent the wrinkling pattern of the system.

4 Results and Discussion

In this section, the numerical results of the instability for a free-standing film and a substrate-bonded film are presented based on the solution approaches discussed in the previous sections including the closed-form solution, the series solution and the numerical finite difference method. For the free-standing film, the buckling loads and mode shapes are determined. For the deposited film on the substrate, the critical wrinkling load is determined from finite difference method and verified by series solution technique. The effect of the thickness non-uniformity of the film on the instability parameters is investigated.

4.1 Buckling of a Free-standing Film

The buckling load N and mode shapes of a free-standing film with variable thickness for symmetric and antisymmetric buckling modes are studied from Equations (10, 11). The variation of the normalized buckling load N/N${}^{\mathrm{0}}_B$ for symmetric and antisymmetric buckling modes are plotted in Figure 3, where $N_B^0=4{\pi }^2$ is the critical buckling load of a free-standing film with uniform thickness in Equation (13a). The symmetric mode with smaller buckling load corresponds with the first mode, while the antisymmetric mode is the second buckling mode of the film. Both the symmetric and antisymmetric buckling loads increase with the growth of the thickness gradient $\epsilon$. When the parameter $\epsilon$ approaches zero, the symmetric and antisymmetric buckling loads approach the corresponding buckling load of the uniform thickness film in Equations (13a, b) as expected.

Figure 3 Normalized buckling load $N/N_B^0$ versus thickness gradient $\epsilon$ for the clamped-clamped free-standing film from analytical and finite difference (F.D.) solution for the symmetric (1st mode) and antisymmetric (2nd mode) buckling.

In order to propose an explicit expression in terms of the thickness gradient $\epsilon$ for the critical buckling load (i.e. the first mode) of the free-standing film with non-uniform thickness, a regression analysis is performed on the buckling load data. According to Figure 3, the critical buckling load $N_B^{Symm.}$ changes linearly versus parameter $\epsilon$ as

$$N_B^{Symm.}=N_B^0(1+m_0\epsilon )$$

(16)

where $m_0=1.273\pm 0.008$ is a constant parameter obtained from the regression analysis [36] with R-Squared $=$ 0.999 and standard error less than 1%, which shows a high accuracy for the proposed linear relation in Equation (16).

On the other hand, the symmetric buckling mode shapes are shown in Figure 4 for different values of parameter $\epsilon$. For the film with uniform thickness (i.e. $\epsilon =0$), the results in Figure 4 correspond with the mode shape in Equation (13a). By growing the thickness gradient $\epsilon$, the edges of the film get thicker. Increasing the bending rigidity of the film decreases the deflection by squeezing the mode shapes as shown in Figure 4.

Figure 4 Symmetric buckling mode shapes of a free-standing film versus length span $\xi$.

4.2 Wrinkling of a Substrate Bonded Film

For a film with uniform thickness (i.e. $\epsilon =0$) deposited on a substrate, the governing equation (5) of the system with clamped-clamped edges is solved analytically [37], and the critical compressive load of the instability is represented versus the Euler buckling load (i.e. $N_B^0=4{\pi }^2$ in Equation (13a)) and the substrate effect as $N=N_B^0+2\sqrt{K}$. Obviously, by increasing the substrate parameter K defined in Equation (6b), especially for thin-film structure with tiny thickness of the film, the effect of the substrate dominates on the system for large parameters K compared with $N_B^0$. Therefore, the film undergoes a fine wavy pattern and the wrinkling load and wrinkling wavenumber are given in terms of K as

$N_W^0$	$=2\sqrt{K}$	(17)
${\beta }_W^0$	$=\sqrt[4]{K}$	(18)

Figure 5 The wrinkling pattern of a film with uniform thickness $(\epsilon =0)$ deposited on a substrate similar to [12, 19].

Same relations for the load and wavenumber of the wrinkling of a thin film-substrate system are presented in the literature by other researchers [11, 12, 38]. Also, the series solution method and finite difference analysis introduced in this work lead to the same results for the load and wavenumber of the wrinkling in uniform thickness films. According to Equations (17) and (18), the wrinkling load and wrinkling wavenumber are affected only by the substrate parameter for big values of parameter K in Equation (6b). Unlike the buckling of a free-standing film where the boundary conditions affect the load and mode shapes, in the wrinkling problem, the boundary conditions only change the wrinkling pattern. As shown in Figure 5-a, for a thin film with uniform thickness deposited on a substrate with clamped-clamped edges, the wrinkling amplitude decreases by approaching the edges. Same results were presented by Damil and Potier-Ferry [19] when they considered a beam with uniform thickness on a substrate with clamped boundaries. Also, the wrinkling pattern of a uniform thickness film with infinite length is shown in Figure 5-b in which the effect of the boundaries is completely vanished and the wrinkles propagate uniformly all over the domain [4, 7, 12].

On the other hand, for a film with non-uniform thickness deposited on a substrate, the critical load and wavenumber of the wrinkling are sought in terms of substrate parameter K and thickness gradient $\epsilon$. The normalized wrinkling load $F_n=N/N_W^0$ is plotted in Figure 6 from finite difference numerical results, where N${}^{\mathrm{0}}$${}_W$ is the wrinkling load for the film with uniform thickness in Equation (17), and K${}_n$ (K${}_n$ $=$ 10${}^{-9}$ K) is the scaled substrate parameter. It is observed that the variation of the normalized wrinkling load F${}_n$ is too small for different values of parameters K and $\epsilon$ such that it remains almost unchanged around the critical wrinkling load of the system with uniform thickness. The small deviations from $F_n=1$ is observed for very soft substrates, where due to decreasing the substrate parameter K, the magnitude of the critical load of the system with uniform thickness decreases too, and therefore, the small deviations are amplified. By considering these small deviations, a relation for normalized wrinkling load F${}_n$ is proposed as

${\mathrm{F}}_{\mathrm{n}}={\displaystyle \frac{\mathrm{N}}{{\mathrm{N}}_{\mathrm{W}}^0}}=1+{\mathrm{m}}_1{\mathrm{K}}_{\mathrm{n}}^{-{\mathrm{m}}_2}{\epsilon }^{{\mathrm{m}}_3}$

(19)

where the constant parameters m${}_{\mathrm{1}}$, m${}_{\mathrm{2}}$ and m${}_{\mathrm{3}}$ are obtained from a regression analysis. The results of the regression analysis with 85 data-points for $0\le \epsilon \le 1$ and $0.001<K_n<100$ are shown in Table 1 with a high accuracy as R-squared $=$ 0.99 and a standard error less than $\pm$ %3 for the estimated parameters [36]. Clearly, imposing $\epsilon =0$ in Equation (19) leads to $F_n=1$ corresponding with the wrinkling load of the film with uniform thickness in Equation (17).

Table 1 The parameters of the Equation (19) for the normalized wrinkling load obtained from regression analysis

Model summary		R Squared $=$ 0.99
			95% Confidence Interval
Parameter	Estimate	Std. Error	Lower Bound	Upper Bound
m${}_{\mathrm{1}}$	.024	.001	.022	.025
m${}_{\mathrm{2}}$	.213	.004	.206	.220
m${}_{\mathrm{3}}$	.775	.020	.736	.814

Figure 6 Normalized wrinkling load $F{}_n$ versus substrate parameter $K_n$ and thickness gradient $\epsilon$.

Furthermore, the wrinkling pattern of the film-substrate system changes effectively by the non-uniformity of the film thickness (Figure 7). For a film with uniform thickness, the wrinkles propagate all over the length span; while for a film with non-uniform thickness, the wrinkles accumulate locally on the film. The wrinkling pattern is characterized by introducing two parameters, the footprint of the wrinkling and the wavenumber. The footprint of the wrinkles indicates the effective length of the film influenced by the wrinkles and the wavenumber of the wrinkling shows the number of the wrinkles on the film [29]. More wavenumber on shorter footprint shows a denser wrinkling region. Due to the non-uniformity of the thickness, wrinkles accumulate around the location with minimum thickness on the film, and by increasing the thickness gradient $\epsilon$, the wrinkles are compressed even more and the number of them decreases too as shown in Figure 7.

Figure 7 Wrinkling pattern of a film with variable thickness and the effect of the thickness gradient $\epsilon$. (a) $\epsilon =0$; (b) $\epsilon =0.2$; (c) $\epsilon =0.8$.

In order to consider the wavenumber and the footprint of the wrinkles quantitatively, the corresponding numerical results are interpreted by a regression analysis. The normalized wavenumber ${\beta }_n=\beta /{\beta }_W^0$ of the wrinkling for the film with variable thickness versus substrate parameter K${}_n$ (K${}_n$ $=$ 10${}^{-9}$ K) and thickness gradient $\epsilon$ is plotted in Fig. 8, where $\beta$${}^0$${}_W$ is the wrinkling wavenumber for a film with uniform thickness in Equation (18). By increasing the thickness gradient $\epsilon$, the wavenumber of the wrinkling decreases effectively, clearly shown with a descending surface in Figure 8. On the other hand, by increasing the substrate parameter, similar to the case of the uniform thickness film in Equation (18), the wavenumber of the wrinkling increases too, but with a slower rate. Therefore, the normalized wavenumber $\beta$${}_n$ decreases by increasing K${}_n$ as Figure 8 shows.

Figure 8 Normalized wavenumber ${\beta }_n=\beta /{\beta }_W^0$ versus substrate parameter $K_n$ and thickness gradient $\epsilon$.

From the numerical results of the finite difference method for a film with non-uniform thickness, an explicit expression of the normalized wavenumber is proposed as

$${\beta }_n=\frac{\beta }{{\beta }_W^0}=\mathrm{EXP}(-m_1{\epsilon }^{m_2})K_n^{-(m_3{\epsilon }^{m_4})}$$

(20)

where $m_i(i=1\mathrm{\dots }4)$ are constant parameters. Clearly, imposing $\epsilon =0$ in Equation (20) simplifies it to the case of the uniform thickness pattern with $\beta ={\beta }_W^0$ in Equation (18) as expected. By increasing the thickness gradient $\epsilon$, the number of the wrinkles decreases exponentially (Figure 8). The effect of the parameter $\epsilon$ on the normalized wavenumber is more important than the substrate parameter K; hence, the corresponding term with the power of K${}_n$ in Equation (20) can be ignored for a simpler but rough approximation. The best approximation of the normalized wavenumber in Equation (20) is obtained for m${}_{\mathrm{4}}\mathrm{=}\mathrm{1}\mathrm{/}\mathrm{3}$. By using a regression analysis with 68 data-points for $0\le \epsilon \le 1$ and $0.01\le K_n\le 10$, the constant parameters m${}_{\mathrm{1}}$, m${}_{\mathrm{2}}$ and m${}_{\mathrm{3}}$ in Equation (20) are determined [36]. The results of the regression analysis presented in Table 2 show a high accuracy for the proposed relation in Equation (20) with R-squared $=$ 0.980 and a standard error about $\pm$%3 for estimated parameters.

Table 2 The regression analysis results of the wavenumber parameters in Equation (20)

Model summary		R Square $=$ 0.980
			95% Confidence Interval
Parameter	Estimate	Std. Error	Lower Bound	Upper Bound
m${}_{\mathrm{1}}$	.841	.017	.808	.874
m${}_{\mathrm{2}}$	.470	.015	.440	.500
m${}_{\mathrm{3}}$	.122	.004	.114	.130

Figure 9 Normalized wavenumber compared with the predicted values of the proposed relation in Equation (20).

Also Figures 9 and 10 show the diagram of the normalized wavenumber obtained from finite difference solution versus the predicted values form Equation (20) and the histogram of the residual errors, respectively. According to the figures, one may conclude that the presented relation in Equation (20) estimates the wrinkling wavenumbers very well. On the other hand, the finite difference results are compared with the results of the series solution in Figure 11(a) and (b) for normalized load and normalized wavenumber, respectively. The wrinkling loads from both methods agree very well with each other; however, the wrinkling wavenumber of the series solution is estimated more than that of the finite difference method. This may happen due to the fact that the finite difference method proposes a stiffer system due to discretization of the continuous model and truncation errors.

Figure 10 Histogram of the residual errors of the Equation (20) for predicting the wavenumber.

Figure 11 Comparison of the finite difference and series solution results.

The change of the substrate parameter K and thickness gradient $\epsilon$ not only alters the number of wrinkles on the film, but also changes the effective length along the span which undergoes wrinkling. As shown in Figure 7, wrinkles accumulate at the location with minimum thickness (here, at the middle of the film). Introducing the footprint of the wrinkling as a dimensionless parameter, one may study the localization of the wrinkles along the length span. The footprint of the wrinkling is defined as the ratio of the length of the film subjected by the wrinkles to the entire length of the film [29] and indicates the effective length of the film influenced by the wrinkles. The footprint varies between zero and one. When the footprint equals to one, the wrinkles propagate along the entire length span of the film corresponding to the wrinkling of the uniform thickness film (i.e. $\epsilon =0$).

A regression analysis between the footprint and the normalized wavenumber of the wrinkling shows that these parameters are linearly dependent (Figure 12) as

$$\mathit{\text{Footprint}}={\mathit{\text{m}}}_0{\beta }_n$$

(21)

and m${}_{\mathrm{0}}\mathrm{=}\mathrm{1.012}\mathrm{\pm }\mathrm{0.007}$ with R-squared $=$ 0.989 for 45 data-points obtained from finite difference solution. Therefore, by increasing the number of the wrinkles on the film, the effective length of the film subjected by the wrinkles increases too.

Figure 12 The linear relation between footprint and normalized wavenumber $\beta$n of the wrinkling.

Figure 13 Footprint of the wrinkles versus thickness gradient $\epsilon$.

The high sensitivity of the footprint parameter (as well as wrinkling wavenumber) with respect to the thickness gradient $\epsilon$ indicates that even with small disturbances of the uniformity of the film thickness, wrinkles accumulate densely at the thinnest location of the film as shown in Figure 13. The tiny thickness of thin film intensifies the importance of the variation of the thickness as well as the above-mentioned accumulative effect. Any small non-uniformity of the thickness of the film which may unavoidably exist on the film leads to the accumulation of the wrinkles around a region; hence, the wrinkling behavior of a substrate-bonded film with non-uniform thickness may significantly differ from the uniform thickness film. Therefore, this work is expected to increase the insight into the wrinkling of a non-uniform film bonded substrate system.

5 Conclusion

The instability problem of a thin solid film with non-uniform thickness along the span is considered by solving the eigenvalue problem of the differential equation of the system. The buckling problem of the free-standing film is solved analytically, and symmetric and antisymmetric modes are compared for uniform and non-uniform films. For a deposited film on a Winkler substrate, the uniaxial wrinkling pattern is investigated and the effect of the substrate stiffness and thickness non-uniformity on the load and pattern of the wrinkling is studied. Characterizing the localization of the wrinkling shows that for a film with uniform thickness, the wrinkles propagate along the entire length span, while for a film with variable thickness the wrinkles accumulate around the thinnest location of the film. By increasing the gradient of the film thickness, the wrinkles shrink more around the thinnest location of the film. The results of the analysis are promising in predicting and controlling the wrinkling pattern in science and technology of thin-film structures. By controlling the thickness profile of the film, desired instability modes appear which can be employed in design and developing various systems including superconducting films, optical systems, sensors, MEMS and NEMS devices.

Appendix

The finite difference formulation used in this work is a central difference approach with second-order accuracy [33]. The differential Equation (5) is casted into an algebraic form as

$R_{i-2}^{(i)}{w}_{i-2}+R_{i-1}^{(i)}{w}_{i-1}+R_i^{(i)}{w}_i+R_{i+1}^{(i)}{w}_{i+1}+R_{i+2}^{(i)}{w}_{i+2}=0$

(A.1)

represented in matrix form by Equation (15), where the matrices [A] and [B] are given by

and the corresponding coefficients are introduced as (i = 2, 3, …, n-1)

$R_{i-2}^{(i)}$	$={\displaystyle \frac{{\Gamma }_4^{(i)}}{h^4}}-{\displaystyle \frac{{\Gamma }_3^{(i)}}{2h^3}}$
$R_{i-1}^{(i)}$	$={\displaystyle \frac{{\Gamma }_2^{(i)}}{h^2}}+{\displaystyle \frac{{\Gamma }_3^{(i)}}{h^3}}-4{\displaystyle \frac{{\Gamma }_4^{(i)}}{h^4}}$
$R_i^{(i)}$	$={\Gamma }_1^{(i)}-2{\displaystyle \frac{{\Gamma }_2^{(i)}}{h^2}}+6{\displaystyle \frac{{\Gamma }_4^{(i)}}{h^4}}$
$R_{i+1}^{(i)}$	$={\displaystyle \frac{{\Gamma }_2^{(i)}}{h^2}}-{\displaystyle \frac{{\Gamma }_3^{(i)}}{h^3}}-4{\displaystyle \frac{{\Gamma }_4^{(i)}}{h^4}}$
$R_{i+2}^{(i)}$	$={\displaystyle \frac{{\Gamma }_3^{(i)}}{2h^3}}+{\displaystyle \frac{{\Gamma }_4^{(i)}}{h^4}}$	(A.4)

where;

${\Gamma }_1^{(i)}$	$=K$
${\Gamma }_2^{(i)}$	$=3\epsilon f_i^{\prime \prime }+6{\epsilon }^2\left(f_i{f}_i^{\prime \prime }+f_i^{\prime }{}^2\right)+3{\epsilon }^3\left(f_i{}^2{f}_i^{\prime \prime }+2f_{i}f_i^{\prime }{}^2\right)+N$
${\Gamma }_3^{(i)}$	$=6\epsilon f_i^{\prime }{(1+\epsilon f_i)}^2$
${\mathrm{\Gamma }}_4^{(i)}$	$={(1+\epsilon f_i)}^3$	(A.5)

and f${}_i$, f${}_i$’ and $f_i^{\prime \prime }$ are evaluated at $\xi ={\xi }_i$ for $f(\xi )$, its first and second derivatives at each nodes. The corresponding parameters for ${\mathrm{\Gamma }}^{(i)}$ and R${}^{\mathrm{(}i\mathrm{)}}$ are obtained from relations (A.4,5) and the matrices [A] and [B] are constructed row by row. Finally, the eigenvalues and eigenvectors of Equation (15) are easily computed using commercial software.

Nomenclature

Symbol	Description
D	Bending stiffness of the film
E	Young’s Modulus
F	Dimensionless loading parameters
H, t	Thickness
K	Substrate Winkler modulus
L	Length
N	In-plane force parameter
b	Width of the system
f( )	Profile of the film thickness
x, $\xi$	length span coordinate
w	Deflection
$\beta$	Wavenumber of wrinkles
$\epsilon$	Thickness gradient parameter
G, $\theta$, A, B, C, R, S, a, c, h, i, m, q, u	Mathematical dummy variables

References

[1] C. Won, H.G. Kim, S. Lee, D. Kim, S. Park, and J. Yoon, Wrinkling prediction for GPa-grade steels in sheet metal forming process, Int. J. Adv. Manuf. Technol. 102 (2019), pp. 3849–63.

[2] S. Shiri, H. Naceur, and J..M. Roelandt, Numerical modelling of sheet metal forming and crashworthiness of laminated steel structures using multi-layered solid-shell elements, European Journal of Computational Mechanics 21, 3–6, (2012), pp. 351–364.

[3] X. Deng, Y. Xu, and C. Clarke, Wrinkling modelling of space membranes subject to solar radiation pressure, Compos. B. Eng. 157 (2019), pp. 266–75.

[4] X. Chen, and J.W. Hutchinson, A family of herringbone patterns in thin films, Scripta Mater. 50 (2004), pp. 797–801.

[5] B. Wang, S. Bao, S. Vinnikova, P. Ghanta, and S. Wang, Buckling analysis in stretchable electronics, npj Flex. Electron. 1 (2017), pp. 1–9.

[6] N. Hamila, F. Hélénon, P. Boisse, and S. Chatel, Simulation of mono-ply and multi-ply woven composite reinforcements forming, European Journal of Computational Mechanics, 17, 5–7 (2008), pp. 919–931.

[7] E. Cerda, Mechanics of scars, J. Biomech. 38 (2005), pp. 1598–1603.

[8] E. Cerda, and L. Mahadevan, Geometry and physics of wrinkling, Phys. Rev. Lett. 90 (2003), pp. 074302.

[9] J. Genzer, and J. Groenewold, Soft matter with hard skin: From skin wrinkles to templating and material characterization, Soft Matter. 2 (2006), pp. 310–323.

[10] S. Wang, J. Song, D. Kim, Y. Huang, and J.A. Rogers, Local versus global buckling of thin films on elastomeric substrates, Appl. Phys. Lett. 93 (2008), pp. 023126.

[11] V. Birman, and C.W. Bert, Wrinkling of composite-facing sandwich panels under biaxial loading, J. Sandw. Struct. Mater. 6 (2004), pp. 217–237.

[12] K. Niu, and R. Talreja, Modeling of wrinkling in sandwich panels under compression, J. Eng. Mech. ASCE. 125 (1999), pp. 875–883.

[13] C.D. Coman, Global asymptotic approximations for wrinkling of polar orthotropic annular plates in tension, Int. J. Solids Struct. 47 (2010), pp. 1572–1579.

[14] Y. Lecieux, and R. Bouzidi, Experimental analysis on membrane wrinkling under biaxial load –Comparison with bifurcation analysis, Int. J. Solids Struct. 47 (2010), pp. 2459–2475.

[15] M. Watanabe, Striped–pattern formation of a thin gold film deposited onto a stretched elastic silicone substrate, J. Polym. Sci. Pol. Phys. 43 (2005), pp. 1532–1537.

[16] R. Gilat, A 3-D thermoelastic analysis of the buckling of a layer bonded to a compliant substrate and related problems, Int. J. Solids Struct. 47 (2010), pp. 2533–2542.

[17] F. Bloom, and D. Coffin, Handbook of Thin Plate Buckling and Postbuckling, Chapman & Hall/CRC, 2001.

[18] C.D. Coman, and A.P. Bassom, Wrinkling of pre-stressed annular thin films under azimuthal shearing, Math. Mech. Solids. 13 (2008), pp. 513–531.

[19] N. Damil, and M. Potier-Ferry, A generalized continuum approach to predict local buckling patterns of thin structures, European Journal of Computational Mechanics 17, 5–7 (2008), pp. 945–956.

[20] S.J. Chapman, Q. Du, and M.D. Gunzburger, A model for variable thickness superconducting thin films, Z. angew. Math. Phys. 47 (1996), pp. 410–431.

[21] R. Hong, W. Shao, J. Ji, C. Tao, and D. Zhang, Thermal annealing induced the tunable optical properties of silver thin films with linear variable thickness, Superlattice Microst. 118 (2018), pp. 170–176.

[22] P. Tanwar, A.K. Panwar, S. Singh, and A.K. Srivatava, Microstructural and optical properties investigation of variable thickness of Tin Telluride thin films, Thin Solid Films. 693 (2020), pp. 137708.

[23] R. Poplausks, D. Jevdokimovs, U. Malinovskis, D. Erts, and J. Prikulis, Variable thickness porous anodic Alumina/metal film bilayers for optimization of plasmonic scattering by nanoholes on mirror, ACS Omega. 3 (2018), pp. 5783–5788.

[24] O. Çiçek, S. Kurnaz, A. Bekar, and Ö. Öztürk, Comparative investigation on electronic properties of metal-semiconductor structures with variable ZnO thin film thickness for sensor applications, Compos. Part B-Eng. 174 (2019), pp. 106987.

[25] F.G. Hone, and F.B. Dejene, Synthesis lead sulphide thin films from tartaric acid chemical bath: Study the role of film thickness on the structural, optical and electrical properties, Thin Solid Films. 692 (2019), pp. 137600.

[26] E.M. Belozubov, V.A. Vasil’ev, and N.V. Gromkov, Thin-film nano- and micro-electromechanical systems – the basis of contemporary and future pressure sensors for rocket and aviation engineering, Meas. Tech. 52 (2009), pp. 739.

[27] E. Zakar, M. Dubey, R. Polcawich, B. Piekarski, R. Piekarz, J. Conrad, and R. Widuta, Study of PZT film stress in multilayer structures for MEMS devices, Mater. Res. Soc. Symp. Proc. 605 (1999), pp. 287.

[28] J.H. Kim, K.L. Jang, K. Ahn, T. Yoon, T.I. Lee, and T.S. Kim, Thermal expansion behavior of thin films expanding freely on water surface, Sci. Rep. 9 (2019), pp. 1–7.

[29] M. Noroozi, and L.Y. Jiang, Buckling and wrinkling of a functionally graded material (FGM) thin film, Int. J. Appl. Mech. 4 (2012), pp. 1250012.

[30] M. Noroozi, Wrinkling of a homogeneous thin solid film deposited on a functionally graded substrate, Struct. Eng. Mech. 74 (2020), pp. 96052.

[31] S. Timoshenko, and S.W. -Krieger, Theory of Plates and Shells, 2nd ed., McGraw-Hill Book Company, 1959.

[32] M. Abramowitz, and I.A. Stegun, Handbook of Mathematical Functions, Dover Publications, New York, 1972.

[33] F.B. Hildebrand, Finite Difference Equations and Simulations, Prentice Hall, New Jersey, 1968.

[34] S. Zhao, and G.W. Wei, Matched interface and boundary (MIB) for the implementation of boundary conditions in high-order central finite differences, Int. J. Numer. Meth. Eng. 77 (2009), pp. 1690–1730.

[35] G.W. Wei, Discrete singular convolution for solution of the Fokker–Planck equations, J. Chem. Phys. 110 (1999), pp. 8930–8942.

[36] L.H. Kahane, Regression Basics, 2nd ed., Sage Publications, 2008.

[37] Column Research Committee of Japan (CRCJ), Handbook of Structural Stability, Corona, Tokyo, pp. A–1–27, 1971.

[38] L. Pocivavsek, R. Dellsy, A. Kern, S. Johnson, B. Lin, K.Y.C. Lee, and E. Cerda, Stress and fold localization in thin elastic membranes, Science. 320 (2008), pp. 912-916.

Biography

Masoud Noroozi received his Ph.D. degree in Mechanical Engineering from University of Western Ontario, Canada. He is currently a faculty member at the Faculty of Mechatronics, Azad University in Iran. His research interests include thin solid film area, mechanical instability, analytical and numerical methods in solid mechanics, and also computational mechanics.

European Journal of Computational Mechanics, Vol. 29_2-3, 255–278.
doi: 10.13052/ejcm2642-2085.29234
© 2021 River Publishers

Abstract

1 Introduction

2 Modeling

3 Solution Approaches

3.1 Buckling Analysis of a Free-standing Film

3.2 Wrinkling Analysis of a Substrate Bonded Film

4 Results and Discussion

4.1 Buckling of a Free-standing Film

4.2 Wrinkling of a Substrate Bonded Film

5 Conclusion

Appendix

Nomenclature

References

Biography