# Analysis of Experimental Biaxial Surface Wrinkling Pattern Based on Direct 3D Numerical Simulation

Seonho Seok, Hyungdal Park, Jinseok Kim

► To cite this version:

Seonho Seok, Hyungdal Park, Jinseok Kim. Analysis of Experimental Biaxial Surface Wrinkling Pattern Based on Direct 3D Numerical Simulation. Micromachines, 2024, 15 (4), pp.543. 10.3390/mi15040543. [hal-04781276](https://hal.science/hal-04781276v1)

HAL Id: hal-04781276

https://hal.science/hal-04781276v1

Submitted on 13 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# Analysis of Experimental Biaxial Surface Wrinkling Pattern Based on Direct 3D Numerical Simulation

Seonho Seok <fnref n="1" /> ID, HyungDal Park <fnref n="2" /> ID and Jinseok Kim <fnref n="2" />,* ID

**Abstract:** This paper presents a direct 3D numerical simulation of biaxial surface wrinkling of thin metal film on a compliant substrate. The selected compliant substrate is a commercial Scotch tape on which a gold metal thin film has been transferred by using low adhesion between the thin metal film and polyimide substrate. Compared with the previous fabrication of a cylindrical thin-film wrinkling pattern, an undulated wrinkling pattern has been implemented by increasing the width of the thin metal film in order to create biaxial straining in the thin film. To understand the wrinkling behavior due to biaxial loading, a simple direct numerical simulation based on material imperfections defined in the compliant substrate has been conducted. Through modeling and simulation, it was found that the wrinkling mode is determined by the biaxiality ratio (BR), the ratio between transversal strain and longitudinal strain. Depending on the BR, the wrinkling mode belongs to one of the cylindrical, undulated (or herringbone), checkerboard, or labyrinth modes as a function of applied strain. The cylindrical wrinkling is dominant at the input of BR less than 0.5, while the undulated (or herringbone) ones become dominant just after the onset of the wrinkling pattern at BR greater than 0.9. Through the comparison of the wrinkling patterns between simulation and experiment, the applied BR of the fabricated thin film has been successfully estimated.

**Keywords:** biaxial wrinkling; buckling; fem; direct simulation; 3D; transfer

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Seok, S.; Park, H.; Kim, J. Analysis of Experimental Biaxial Surface Wrinkling Pattern Based on Direct 3D Numerical Simulation. *Micromachines* **2024**, *15*, 543. https://doi.org/10.3390/mi15040543

Academic Editor: Peng Gao

Received: 28 February 2024

Revised: 13 April 2024

Accepted: 17 April 2024

Published: 18 April 2024

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Buckling of multi-layered material is a classical subject and has received continuous attention in mechanics due to its various applications, such as sandwich panels. Concerning microfabrication, such buckling has been frequently observed when thin and stiff film is boned to a compliant substrate. This can happen experimentally by depositing a thin metal film onto a compliant substrate, such as an uncured polymer substrate, and letting the system cool down to room temperature. Surface wrinkling driven by mechanical instability is commonly observed in different thin-film structures with a compliant substrate [1–5]. Apart from the metal deposition mentioned earlier, such surface wrinkling patterns have been engineered with different techniques such as surface oxidation of PDMS [6], thin metal film deposition on thermally expanded polymer [7], polymer mold confinement on buckling of thin film [8], etc. Significant progress in theoretical and experimental studies has led to a greater physical understanding of surface instability phenomena, spurring research interest in engineering topologically self-organizing systems across many length scales [9]. Recently, surface wrinkling based on the thin-film buckling phenomenon has been proposed to find the elastic moduli of thin films such as thin polymer films, thin film metals, carbon nanotubes, etc. [10]. This technique allows us to calculate the elastic modulus of a thin film by measuring the wavelength of the wrinkling thin film. Elastomer material such as PDMS is frequently used as a base substrate as it has a low elastic modulus.

<footnote n="1">¹ Center for Nanoscience and Nanotechnology (C2N), Université Paris-Saclay, 91129 Palaiseau, France; seonho.seok@c2n.upsaclay.fr</footnote>

<footnote n="2">² Center for Bionics, Korea Institute of Science and Technology, Seoul 02792, Republic of Korea; hyungdal@kist.re.kr</footnote>

<footnote>* Correspondence: jinseok@kist.re.kr</footnote>
Commercial Scotch tape having adhesive on a supporting membrane such as PVC, PET, and PE was also reported as a buckling substrate with an elastic thin film [11].

For most of the reports on wrinkling studies and material properties extraction from wrinkling patterns, the analytical formulation has been widely employed in analyzing experimental data, for example, wrinkling wavelengths [12–14]. Additionally, numerical studies have been carried out to explore the underlying buckling mechanism as well as to quantitatively control and manipulate the self-assembled patterns for applications in micro/nanofabrication [15–17]. Modeling and simulation based on finite element modeling (FEM) are useful to analyze the surface wrinkling pattern depending on various geometries, different material properties, etc. The FEM models based on 2D or 3D have successfully provided wrinkling wavelengths as a function of material properties associated with surface wrinkling. As for finite element modeling associated with instabilities, it is conventional that linear buckling for wrinkling wavelength based on basic buckling mode and post-buckling for wrinkling amplitude based on buckling mode from pre-performed linear buckling analysis have been used in the studies of the surface wrinkling phenomenon [18–20]. These simulations have been, in general, realized by buckling analysis established in commercial software packages such as ABAQUS and ANSYS. Geometrical imperfection methods have been frequently employed by using mesh, geometry, and boundary condition perturbation techniques for studying surface instability. However, the existing imperfection approaches are relatively complicated, and the implementation can be laborious due to free parameter calibration. Furthermore, the interpretation and verification of results are not straightforward, which makes it less practical for common users. To overcome the drawbacks of the geometrical imperfection methods, a direct FEM simulation method that uses embedded imperfections with perturbed material properties at the film-substrate interface has been recently proposed. The direct simulation can successfully perform the pre/post-buckling simulations in only one analysis step with any common finite element code and analysis platform [21–23]. The direct FEM model can provide the wrinkling wavelength of a surface wrinkling as a function of material properties, and thus the extraction of material properties such as Young's modulus has been achieved [24]. In this paper, a direct 3D numerical simulation of biaxial surface wrinkling of thin metal film on a compliant substrate is presented to comprehend the experimentally created wrinkling pattern. The theoretical background of thin-film wrinkling on the compliant substrate is explained in Section 2, and the fabrication process of the wrinkling pattern is depicted in Section 3. Section 4 shows the results and discussions of the FEM modeling and simulation through comparison with experimental results. Finally, the conclusion will be made in Section 5.

## 2. Brief Theoretical Description of Wrinkling of Thin Film on a Substrate

The force balance approach for surface buckling instability of a thin film on a compliant substrate is here briefly introduced [10]. The wrinkling problem can be described such that a thin elastic film is bound to an elastic foundation, as shown in Figure 1. The source of the wrinkling is the compressive stress of the film.

Considered as a semi-infinite substrate under plane strain deformation, the classical equation for the bending of an elastic film on a compliant elastic substrate is given by Equation (1).

$$ \overline{E}_f I \frac{d^4 z}{dx^4} + F \frac{d^2 z}{dx^2} + k_m z = 0 \qquad (1) $$

where $\overline{E} = E/(1 - v^2)$ is the plane-strain modulus, $E$ is Young's modulus, $v$ is Poisson's ratio, $I = \text{wh}^3/12$ is the moment of inertia (where $w$ is the width of the film and $h$ is its thickness), $F$ is the uniaxially applied force or load, and $k$ is Winkler's modulus of an elastic half-space ($k_m = \overline{E}_s w \pi / \lambda$). The subscripts, $f$ and $s$, denote the film and substrate, respectively.

In Equation (1), the first and third terms are along with the classical Euler-Bernoulli beam-bending equation that equalizes bending forces in the film with the normal force distribution due to the deformation of the substrate above and below its neutral position. The second term means the effects of the force, $F$, on the film. In contrast to the
buckling of a column, the wrinkling of a film on an elastic foundation involves a balance of force between film bending (first term in Equation (1)), which acts to suppress short wavelengths, and substrate deformation (third term in Equation (1)), which acts to suppress large wavelengths [10].

<figure><img src="image_5.png" /><figcaption>Figure 1. Geometrical schematics of the thin-film wrinkling problem.</figcaption></figure>

As the buckling instability of interest here is the first sinusoidal mode, the film deflection can be described by Equation (2).

$$w(x) = A\cos \frac{2\pi x}{\lambda} = A\cos kx \quad (2)$$

Substituting Equation (2) into Equation (1) and solving for the applied force in the thin film on a compliant substrate gives Equation (3).

$$F = 4\overline{E_f}I\left(\frac{\pi}{\lambda}\right)^2 + \frac{\overline{E_s}w}{4}\left(\frac{\pi}{\lambda}\right)^{-1} \quad (3)$$

The film buckling wavelength can be found by minimizing F with respect to λ (or, $\partial F/\partial\lambda = 0$),

$$\lambda = 2\pi h \left( \frac{\overline{E_f}}{3\overline{E_s}} \right)^{\frac{1}{3}} \quad (4)$$

It should be noted here that the wavelength is only a function of the thickness of the film and the elastic properties of the film and substrate. Thus, the wrinkling wavelength can be used to determine the material properties of the thin film if the substrate material properties are known.

Another important parameter is the critical stress or strain needed to induce wrinkling in the system. The critical stress can be found in Equations (4) and (5) by dividing the critical force (Fc) by the cross-sectional area of the thin film.

$$\sigma_c = \frac{F_c}{hw} = \left( \frac{9}{64} \overline{E_f} \overline{E_s}^{-2} \right)^{\frac{1}{3}} \quad (5)$$

Therefore, the critical strain is given below.

$$\epsilon_c = \frac{\sigma_c}{\overline{E_c}} = \frac{1}{4} \left( \frac{3\overline{E_s}}{\overline{E_f}} \right)^{\frac{2}{3}} \quad (6)$$
Under the assumption that the wrinkling wavelength is independent of the applied strain, the wrinkling amplitude can be found in Equation (7).

$$A = h\sqrt{\frac{\varepsilon - \varepsilon_c}{\varepsilon_c}} \qquad (7)$$

where A is wrinkling amplitude, and ε is applied strain.

In the case of biaxial straining, an infinitesimal perturbation can be introduced to the deflection of the film as follows [25]:

$$w(x, y) = A\cos kx + \delta w(x, y) \qquad (8)$$

Note that the first term is the same one in uniaxial straining, as explained above.

By applying the wrinkling behaviors in the x and y axes, the perturbation equation can be expressed as follows [25].

$$\delta w(x, y) = A_1\cos(k_x x + k_y y) + A_2\cos(k_x x - k_y y) + A_3\cos(k_y y) \qquad (9)$$

where k_x and k_y are the associated wavenumbers, and A_1 is the arbitrary amplitude.

Using trigonometric identities, the deflection of the film, including the base solution Acos(kx), can be then written as:

$$\delta w(x, y) = A\cos(k_x x) + b\sin(k_x x)\sin(k_y y) + c\cos(k_x x)\cos(k_y y) + d \cos(k_y y) \qquad (10)$$

where q = k_y/k_x is the rescaled longitudinal wavevector, b = A_2 - A_1, c = A_1 + A_2, and d = A_3. It is noted that the numbers b, c, and d appear as amplitudes of the undulating mode, the varicose mode, and the checkerboard mode, respectively. The amplitudes b, c, and d are taken to be infinitesimal, and thus, linear stability is not considered in this present form.

Biaxial wrinkling patterns depending on parameter q are given in Figure 2. The wrinkling pattern has been defined as λ_x = 50 µm and λ_y = 250 µm. A cylindrical wrinkling pattern is found when q is zero, as expected. When q is 0.25, it is confirmed that the 4-cycle wrinkling pattern is established in the x direction, while the 1-cycle pattern is made in the y direction by definition of the parameter q. The other wrinkling patterns have shown the same tendency, depending on the q.

<figure><img src="image_6.png" /><figcaption>Figure 2. Cont.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 2.** Wrinkling patterns as a function of parameter q.</figcaption></figure>

**3. Fabrication of Biaxial Thin-Film Wrinkling on Scotch Tape**

To make metal wrinkling patterns on Scotch tape, thin metal films have been deposited on a PSPI (Photosensitive Polyimide) layer built on Si substrate by e-beam evaporation (ei-5k, ULVAC, Chigasaki, Japan). As in the previous work, the PSPI layer has been used because it has a low adhesion force to the metal films. Figure 3 shows the fabrication process of the test sample.

<figure><img src="image_8.png" /><figcaption>**Figure 3.** Thin-film metal test pattern fabrication process.</figcaption></figure>

First, the native oxide and organic matter formed on the Si wafer surface were removed through a general piranha cleaning solution, which is a 3:1 mixture of sulfuric acid and 30% hydrogen peroxide. Thereafter, the PSPI was coated to a thickness of 5 µm with a 4500 rpm/1 min condition using a spin coater. The coated PSPI was soft-baked at 110 °C for 6 min, and 200 mJ UV was exposed to PSPI using the aligner equipment (MA-6 Mask Aligner, Karl Suss, Munich, Germany). Afterward, PEB (Post-Exposure Baking) was performed at 110 °C for 5 min, and a pattern was developed using a dedicated developer (401D, HD MicroSystems™, Parlin, NJ, USA) and Rince (AZ Thinner 1500, MERCK, Rahway, NJ, USA). The patterned PSPI was subjected to PDB (Post-Develop Baking) at 180 °C for at least 1 min, and full curing was performed for 30 min at 200 °C and at 300 °C for 1 h in a vacuum drying oven. (Customized Oven, SH Scientific Co., Ltd., Sejong, Republic of Korea).

Au (gold) was deposited on a full-cured PSPI to a thickness of 300 nm using an E-beam evaporator. Subsequently, a commercial thermal adhesive film was applied to safeguard the Au boundary with the Si surface formed at the periphery of the Si wafer. Patterning was performed on the deposited metal thin film using a positive photoresist (GXR-601 46cp, AZ Electronics Materials, Hsinchu Country, Taiwan), and Au was selectively dry etched
using an ICP-RIE (Inductively Coupled Plasma Reactive Ion Etcher, Oxford Instruments, Abingdon, UK). Considering the low adherence of the formed Au pattern to PSPI, the photoresist removal was completed by using O₂ plasma at 300 W for 5 min using Asher equipment (Plasma Finish V15-G, Ebhausen, Germany). Figure 4 shows an Au thin film fabricated on the PSPI polymer with a Si substrate.

<figure><img src="image_9.png" /><figcaption>**Figure 4.** Fabrication result of the thin-film metal test pattern.</figcaption></figure>

The gold thin film was transferred onto Scotch tape by commercial tensile experimental equipment (Shimadzu EZ-S machine, Shimadzu, Kyoto, Japan), and the debonding force was measured at the same time, as shown in Figures 5 and 6, respectively. The maximum applied force is about 0.5 N, as there is no adhesion layer such as Ti or Cr. After metal transfer, metal film wrinkling patterns have been established on the Scotch tape due to a mechanical property mismatch between the two materials.

<figure><img src="image_10.png" /><figcaption>**Figure 5.** Gold thin film from PSPI substrate to Scotch tape.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 6.** Debonding force measurement of the Scotch tape.</figcaption></figure>

Figure 7 shows the undulated wrinkling pattern and the cylindrical one transferred to the Scotch tape. The lines between the undulated wrinkles are created by tape tests due to the angled tape attached to the metal pattern. The wrinkling patterns have been characterized with an optical microscope and SEM (scanning electron microscopy) in order to measure the wrinkling wavelength. The average values for the wavelengths in the direction of *x* and *y* are 57.8 µm (standard deviation of 0.62 µm) and 1115 µm (standard deviation of 78 µm), respectively. It should be noted that the wavelength of the cylindrical wrinkling, shown in Figure 7b, was measured at 55 µm in the previous work. The ratio of wavelengths in the *x* and *y* directions, L<sub>*x*</sub>/L<sub>*y*</sub>, is given by 0.052, which means the wrinkling in the *y* direction occurs just after the critical strain because it is close to *q* = 0.

<figure><img src="image_12.png" /><figcaption>**Figure 7.** Wrinkling patterns of thin metal film are transferred onto Scotch tape, depending on the design dimensions. (a) Undulated wrinkling pattern; (b) Cylindrical wrinkling pattern; (c) dimension of thin metal film for undulated and cylindrical.</figcaption></figure>

**4. Three-Dimensional Direct FEM Modeling, Simulation and Results**

A model consisting of a thin film bonded to a compliant substrate has been built with an embedded imperfection for 3D direct buckling simulation, as shown in Figure 8. It shows the boundary conditions and the applied displacements (or strains). The film surface is basically square or rectangular in shape depending on the loading conditions, such as
equibiaxial or biaxial stressing, and the substrate was sufficiently thick in comparison with the film.

<figure><img src="image_13.png" /><figcaption>Figure 8. Three-dimensional wrinkling model and boundary conditions. (a) Overall view; (b) Top view.</figcaption></figure>

The thin-film thickness is defined as 0.5 µm. Both the film and substrate materials are taken to be isotropic linear-elastic. Material properties needed for thin film and compliant substrates are Young's moduli and Poisson ratios. The compliant substrate has 6.9 MPa and 0.45, while the thin film has 1 GPa and 0.42, respectively. One embedded imperfection, defined by the elastic property of the thin-film material instead of that of the substrate, is placed in the substrate adjacent to the film-substrate interface. To preserve symmetry, the imperfection with the square shape is put at the center of the plane. The dimension of the model has been chosen to be sufficiently large for the formation of different wrinkling patterns, such as cylindrical and undulated, as explained in the previous section. The boundary conditions $U_x = 0$, $U_y = U_z = Free$ at $x = L_x$, $U_y = 0$, and $U_x = U_z = Free$ at $y = L_y$ have been defined as symmetric planes, respectively. The vertex at the center in the bottom plane, $z = 0$, is regarded as fixed. Displacements at $x = 0$ and $y = 0$ are applied for the external loads. As a quarter model is used for the simulation, the results have been extended into a full model.

Biaxial surface wrinkling has been simulated by applying transversal strain ($d_x$) and longitudinal strain ($d_y$) at the same time, and thus, biaxiality ratio (BR), defined as Equation (9), is used as a parameter of simulation with transversal strain. Simulations for various biaxiality ratios were performed in such that $d_x$ was kept constant and $d_y$ was varied within the range of $0 \le d_y \le d_x$ for different cases of BR. As a consequence, the biaxiality varies within the range of $0 \le BR \le 1$, where BR = 0 assimilates with uniaxial compression (applied in the x-direction) and BR = 1 corresponds to perfectly equibiaxial compression.

$$ \text{Biaxiality Ratio (BR)} = \frac{d_y}{d_x} \qquad (11) $$

Verification of the model has been performed with the case of uniaxial compression (BR = 0), as it is considered the initial wrinkling state for the following complex wrinkling patterns. The wavelength variation resulting from the simulation depending on the mesh size has been first investigated, and then wrinkling behaviors have been simulated. Figure 9a shows the cylindrical wrinkling pattern due to the uniaxial straining and its wrinkling characteristics. The wavelength has little dependency on the mesh size, as presented in Figure 9b. It is noteworthy that finer meshes provide a more stable wavelength value. By referring to maximum stress as a function of applied strain ($\epsilon_a$), critical strains initiating the cylindrical wrinkling can be found, as shown in Figure 9c. Figure 9d presents wrinkling parameters in terms of wavelength and amplitude, which can be used for material property extraction of the thin film. Note that the wrinkling wavelength from 3D direct simulation
is greater than that of 2D direct and analytical solutions. As the elasticity of the thin film ($E_f$) increases, the critical strain becomes smaller.

<figure><img src="image_14.png" /><figcaption>**Figure 9.** Cylindrical wrinkling and its wrinkling characteristics. (a) Cylindrical wrinkling; (b) Wrinkling wavelength vs. mesh size; (c) Critical strain vs. $E_f$; (d) Wrinkling parameters vs. $E_f$.</figcaption></figure>

Before the wrinkling behavior study of the model with a full range of biaxial loading, the wrinkling pattern change was investigated in the case of $BR = 0.65$. Figure 10 shows force reactions on the surfaces where longitudinal and transversal displacement loading are applied as a function of the applied strain. The entire wrinkling pattern, from the pre-instability stage to the final stage, has been captured in one simulation run. The reaction force curves introduce sudden changes in slope or load drop whenever the wrinkling mode change occurs. Starting from the flat surface of the thin film before the critical strain, slight thin-film wrinkling has occurred at the critical strain around imperfection sites. The onset of the initial wrinkling is at the applied strain of approximately $1.8 \times 10^{-2}$. The cylindrical wrinkling mode appeared at the applied strain of approximately $2.8 \times 10^{-2}$ and lasted through the strain of approximately $3.2 \times 10^{-2}$ before the onset of the undulated (or herringbone) mode. From a strain of approximately $4.4 \times 10^{-2}$, the wrinkling pattern has been changed to labyrinth mode. Depending on the wrinkling mode, a longitudinal or transversal force reaction can be used to refer to the wrinkling mode change.

To assess the dependency of wrinkle patterns on the biaxiality ratio, extensive simulations have been performed by monitoring the pattern evolution. The resulting wrinkling patterns are presented in two parts, as a more detailed examination is required in the cases of $BR$ between 0.9 and 1.0 due to the nature of thin-film buckling. Figure 11 shows the top view of the wrinkling patterns for biaxiality ratios between 0.3 and 0.8. It is found that cylindrical wrinkling is the dominant mode until a biaxiality ratio of 0.5 in the range of the longitudinal strain of interest. In the case of $BR = 0.65$, 1D cylindrical wrinkling has been changed to a slightly undulated model from $2.8 \times 10^{-2}$ strain input. The 1D cylindrical wrinkling has disappeared from 0.75 of $BR$ as the undulated mode is dominant after the
critical strain, and it becomes labyrinth mode at the highest strain input. The color contours represent the out-of-plane displacement, with the red and yellow (or orange) colors being the highest (peak) and lowest (valley) positions, respectively.

<figure><img src="image_15.png" /><figcaption>**Figure 10.** Force reaction at the loading surfaces when BR is 0.65.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 11.** Wrinkling patterns for biaxiality ratios between 0.3 and 0.8.</figcaption></figure>

The evolution of the wrinkle patterns for BR values greater than 0.9 is shown in Figure 12. It can be seen that 1D cylindrical wrinkling dominant in the BRs lower than 0.9 is suppressed, and undulated or checkerboard-like modes appear just after critical stains. Note that the critical strain value becomes slightly smaller from the BR of 0.92, and initial wrinkling patterns tend to be circular near the imperfection sites with the BR value greater. These primary modes change more complex patterns or more continuous wrinkles, such as various labyrinth structures, as the applied strain increases.
<figure><img src="image_17.png" /><figcaption>**Figure 12.** Wrinkling patterns for biaxiality ratios between 0.9 and 1.</figcaption></figure>

Given various wrinkling patterns depending on the BR, undulated wrinkles close to the experimental result can be extracted, as shown in Figure 13. From the experiment, it is seen that the undulated wrinkle has two different regions: the undulated pattern itself and the arbitrary pattern. The arbitrary pattern is near the boundary line created due to the angle between the tape and the wafer during the tape test. The closest pattern to the experiment in view of the undulated one and the arbitrary one is the wrinkling pattern with BR 0.92, as it clearly shows two different zones due to the boundary conditions. It is also confirmed that the wavelength of the longitudinal direction is much shorter than that of the transversal one, as are the experimental results. The simulation result with a BR of 0.92 gives wavelengths in longitudinal and transversal directions of 10 µm and 50 µm, respectively. Thus, the q value is estimated at 0.2, which is much greater than the experimental one, 0.052.

<figure><img src="image_18.png" /><figcaption>**Figure 13.** Comparison between experiment and undulated simulation patterns of different BRs.</figcaption></figure>

**5. Conclusions and Perspectives**

Thin-film wrinkling due to a mechanical mismatch between elastic thin metal film and a compliant substrate is frequently encountered in microfabrication. A biaxial wrinkling pattern of a thin film has been fabricated with the transfer of thin metal film into Scotch tape by adjusting the width of the thin film from our previous cylindrical wrinkling metal pattern design. To understand the experimental pattern, 3D direct FEM modeling and simulation based on the embedded imperfection approach have been carried out. The technique is easily implemented for structures consisting of thin films above a compliant substrate without treating the surface patterns discretely under specific assumptions. The generation
of wrinkling morphologies and their transformations have been directly obtained from the simulations. The established model has successfully produced the cylindrical, 1D wrinkling pattern with one directional strain input. The biaxial wrinkling has been performed with a parameter defined as BR (biaxiality ratio). It has been found that the 1D wrinkling pattern is dominant when BR approaches 0.9, while the undulated pattern appears at 0.9 of BR. Given that the BR is larger than 0.9, the wrinkling patterns change to cylindrical mode, undulated mode, or checkerboard mode for equibiaxial straining and labyrinth mode. From the undulated patterns depending on the BR, the wrinkling pattern closest to the experiment has been found to be able to estimate the applied strain. It is worth noting that all the results presented here are based on simple isotropic linear elastic material behavior with an initially perfectly flat film layer. The 3D direct model and simulation should be simple and robust for future investigations involving more complex material, geometric, and loading conditions.

**Author Contributions:** Investigation, H.P.; Writing—original draft, S.S.; Funding acquisition, J.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the Bio & Medical Technology Development Program of the National Research Foundation (NRF) funded by the Ministry of Science and ICT (Project No. 2022M3E5E9016506), the Pioneer Research Center Program through the National Research Foundation of Korea funded by the Ministry of Science and ICT (Project No. 2022M3C1A3091627), and the Institutional Program (Project No. 2E33141) of the Korea Institute of Science and Technology (KIST).

**Data Availability Statement:** The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Wang, S.; Song, J.; Kim, D.H.; Huang, Y.; Rogers, J.A. Local versus global buckling of thin films on elastomeric substrates. *Appl. Phys. Lett.* **2008**, *93*, 023126. [CrossRef]

2. Ma, Y.; Xue, Y.; Jang, K.I.; Feng, X.; Rogers, J.A.; Huang, Y. Wrinkling of a stiff thin film bonded to a pre-strained, compliant substrate with finite thickness. *Proc. R. Soc. A Math. Phys. Eng. Sci.* **2016**, *472*, 20160339. [CrossRef] [PubMed]

3. Breid, D.; Crosby, A.J. Surface wrinkling behavior of finite circular plates. *Soft Matter* **2009**, *5*, 425–431. [CrossRef]

4. Barth, J.V.; Costantini, G.; Kern, K. Engineering atomic and molecular nanostructures at surfaces. *Nature* **2005**, *437*, 671–679. [CrossRef] [PubMed]

5. Yoo, P.J.; Suh, K.Y.; Kang, H.; Lee, H.H. Polymer elasticity-driven wrinkling and coarsening in high temperature buckling of metal-capped polymer thin films. *Phys. Rev. Lett.* **2004**, *93*, 034301. [CrossRef] [PubMed]

6. Bowden, N.; Huck, W.T.S.; Paul, K.E.; Whitesides, G.M. The controlled formation of ordered, sinusoidal structures by plasma oxidation of an elastomeric polymer. *Appl. Phys. Lett.* **1999**, *75*, 2557–2559. [CrossRef]

7. Bowden, N.B.; Brittain, S.T.; Evans, A.G.; Hutchinson, J.W.; Whitesides, G.M. Spontaneous Formation of Ordered Structures in Thin Films of Metals Supported on an Elastomeric Polymer. *Nature* **1998**, *393*, 146–149. [CrossRef]

8. Yoo, P.J.; Suh, K.Y.; Park, S.Y.; Lee, H.H. Physical Self-Assembly of Microstructures by Anisotropic Buckling. *Adv. Mater* **2002**, *14*, 1383–1387. [CrossRef]

9. Chen, X.; Hutchinson, J.W. Herringbone buckling patterns of compressed thin films on compliant substrates. *J. Appl. Mech.* **2004**, *71*, 597–603. [CrossRef]

10. Chung, J.Y.; Nolte, A.J.; Stafford, C.M. Surface Wrinkling: A Versatile Platform for Measuring Thin-Film Properties. *Adv. Mater.* **2011**, *23*, 349–368. [CrossRef]

11. Seok, S.; Park, H.D.; Kim, J. Scotch-tape surface wrinkling based thin-film material properties extraction. *J. Micromech. Microeng.* **2022**, *32*, 045002. [CrossRef]

12. Song, J.; Jiang, H.; Choi, W.M.; Khang, D.Y.; Huang, Y.; Rogers, J.A. An analytical study of two-dimensional buckling of thin films on compliant substrates. *J. Appl. Phys.* **2008**, *103*, 014303. [CrossRef]

13. Stafford, C.M.; Harrison, C.; Beers, K.L.; Karim, A.; Amis, E.J.; VanLandingham, M.R.; Kim, H.C.; Volksen, W.; Miller, R.D.; Simonyi, E.E. A buckling-based metrology for measuring the elastic moduli of polymeric thin films. *Nat. Mater.* **2004**, *3*, 545–550. [CrossRef]

14. Stafford, C.M.; Harrison, C.; Karim, A.; Amis, E.J. Measuring the modulus of polymer films by strain-induced buckling instabilities. *Polymer* **2002**, *43*, 1335.

15. Yin, J.; Chen, X.; Sheinman, I. Anisotropic buckling patterns in spheroidal film/substrate systems and their implications in some natural and biological systems. *J. Mech. Phys. Solids* **2009**, *57*, 1470–1484. [CrossRef]
16. Huang, Z.; Hong, W.; Suo, Z. Evolution of wrinkles in hard films on soft substrates. *Phys. Rev. E* **2004**, *70*, 030601.

17. Yin, J.; Chen, X. Elastic buckling of gradient thin films on compliant substrates. *Phil. Mag. Lett.* **2010**, *90*, 423–433. [CrossRef]

18. Huang, Q.; Xu, R.; Liu, Y.; Hu, H.; Giunta, G.; Belouettar, S.; Potier-Ferry, M. A two-dimensional Fourier-series finite element for wrinkling analysis of thin films on compliant substrates. *Thin-Walled Struct.* **2017**, *114*, 144–153. [CrossRef]

19. Mei, H.; Landis, C.M.; Huang, R. Concomitant wrinkling and buckle-delamination of elastic thin films on compliant substrates. *Mech. Mater.* **2011**, *43*, 627–642. [CrossRef]

20. Leifer, J.; Belvin, W. Prediction of wrinkle amplitudes in thin film membranes using finite element modeling. In *Proceedings of the 44th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Norfolk, VA, USA, 7–10 April 2003*.

21. Siavash, N.; Ryu, D.; Shen, Y.-L. Direct numerical simulation of buckling instability of thin films on a compliant substrate. *Adv. Mech. Eng.* **2019**, *11*, 1687814019840470.

22. Siavash, N.; Ryu, D.; Shen, Y.-L. Instabilities of thin films on a compliant substrate: Direct numerical simulations from surface wrinkling to global buckling. *Sci. Rep.* **2020**, *10*, 5728.

23. Nikravesh, S.; Ryu, D.; Shen, Y.-L. Surface Wrinkling versus Global Buckling Instabilities in Thin Film-Substrate Systems under Biaxial Loading: Direct 3D Numerical Simulations. *Adv. Theory Simul.* **2022**, *5*, 2200183. [CrossRef]

24. Seok, S.; Park, H.; Coste, P.; Kim, J. Direct Numerical Simulation of Surface Wrinkling for Extraction of Thin Metal Film Material Propertie. *Micromachines* **2023**, *14*, 747. [CrossRef]

25. Audoly, B.; Boudaoud, A. Buckling of a stiff film bound to a compliant substrate—Part I Formulation, linear stability of cylindrical patterns, secondary bifurcations. *J. Mech. Phys. Solids* **2008**, *56*, 2401–2421. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.