## 1. Introduction

In nature, there are several examples such as lotus leaves, peacock’s feathers, butterfly wings, gecko feet, and moths’ eyes where meso-nano-scaled physical patterns on the surface impart characteristic properties such as hydrophobicity, structural color, adhesiveness, and optical properties, which are absent otherwise, to the surfaces. To engineer such surfaces artificially and to fabricate such features, there are mainly two most commonly used approaches: one a top-down approach and the other a bottom-up approach. The top-down approach consists of the various lithographic techniques such as electron beam lithography, dip-pen lithography, embossing lithography, laser lithography [1], laser photoablation [2], ion implantation [3], micromachining [4], etc. The bottom-up approach, on the other hand, consists of self-assembly of molecules to produce bigger building blocks to a growing structure as in chemical vapor deposition [5]. While the self-assembly procedure involves numerous steps, lithographic methods such as dip-pen lithography, which are used for creating patterns on hard silicon surfaces, are time-consuming, involve great cost and may not at all be suitable to create patterns on soft surfaces such as polymers. Self-organization methods, on the other hand, are for bulk production of desired structures involving lower cost and therefore are more potentially applicable fabrication techniques on soft surfaces. In self-organization technique, materials such as polymeric films, which have highly tunable mechanical properties, can be used to obtain a desired mesoscopic structure due to the influence of internal or external forces. The mechanical and interaction properties of these films can be controlled by the amount of cross-linker added during the synthesis step itself. The instabilities during self-organization can be harnessed to engineer desired periodic arrangements of sub-micron features in technological applications such as in flexible electronics [6], optoelectronics [7], functional coatings [8], sensors [9], microfluidic devices, microelecromechanical systems (MEMS) and nanoelectromechanical systems (NEMS) [10], pressure-sensitive adhesives [11], high-efficiency light-emitting diodes, thin-film transistors, etc. Thus, these thin films are model-mesoscopic systems for understanding physical phenomena such as tunable adhesion, adhesion debonding, pattern formation, wetting dewetting and friction at soft interfaces, etc.

When a thin elastic film (shear modulus of <1 MPa) is brought in contact proximity of an external contactor, it forms miniaturized surface patterns because of a competition between the destabilizing interaction forces between the film and the contactor and the restoring elastic forces present in the bulk of the elastic film. The wavelengths of instabilities scale as 4 h in the case of peeling geometry [12–14] (whereas glass slide rests on a rigidly bonded film in cantilever configuration) and application of normal force at the end of the glass slide leads to well-defined one-dimensional (1D) finger-like patterns at the crack opening and 3 h in the case of adhesive geometry [15–26]. In the adhesive-contact geometry, the free surface of the film in contact proximity with an external contactor spontaneously roughens to form a two-dimensional (2D) labyrinthic pattern when the stiffness of the attractive forces due to van der Waals (VDW) interactions (arising because of the contactor) exceeds the elastic stiffness of the film. This pattern remained unchanged in experiments either when the upper glass slide was silanized with a chemisorbed layer of fluorinated monochlorodecyl silane [15] or when the interactions were switched to electrostatic ones [21, 24, 25]. In recent past, much effort has been devoted to understand how further miniaturization of length scales and hence enhanced functional properties can be made technologically feasible in these already short wavelength exhibiting adhesive soft elastic thin films in contact geometry. From a fabrication point of view, the smallest patterns realized till date at such soft interfaces are ~ 0.1 h and are obtained with elastic bilayers when the surface energies of the films are very low and the top film is considered to be much thinner and of very small shear modulus compared to the bottom film [27, 28]. But pattern formation in elastic bilayers through squeezing modes of instabilities has the disadvantage of delaminating at the film-film interface. To alleviate this problem, one has found that films cast on patterned substrates can give rise to surface patterns, which are at length scales that are about an order of magnitude less than that found in films cast on flat substrates [29–31].

Apart from direct experiments, these systems are well analyzed both theoretically and numerically. In both cases, atomistic scale models can be considered but are computationally expensive and are restricted by the domain geometries that can be handled by them. For elastic films, linear stability analysis is generally carried out to understand the critical length scales at the incipience of instability by perturbing Navier’s equations with bifurcation modes and obtaining the nontrivial solutions [16, 17]. In nonlinear analysis [18–20], on the other hand, it is possible to find the surface morphology by trying to find the equilibrium displacements that minimize the total energy, which is composed of the interaction, and the elastic energies. Finding the optimized Fourier coefficients, which define the top surface, with the help of conjugate-gradient method, does this. However, for elastic films cast on generally patterned substrates the obtained equations show highly nonlinear behavior, which cannot be solved by the simplistic analytical methods or even by the energy-minimizing schemes developed for the base case. On top of these, polymeric films particularly demonstrate a rich variety of nonlinear behavior ranging from viscoelasticity to plasticity that cannot be adequately captured by the aforementioned methods. Finite element analysis (FEA) methods can be used to alleviate these problems. FEA methods have been used extensively in recent years to study crack opening, crack propagation [32–34], delamination [35], film instability over soft substrates, and bifurcation analysis of elastic contact instability [36–39].

In this chapter we show how to implement cohesive-zone modeling (CZM) for simple shear and normal modes of crack opening in the finite element software and how to simulate contact instability and adhesion debonding in soft thin films cast on complex geometries.

## 2. Cohesive-zone modeling

Most of the failures in engineering structures are the result of cracks that are preexisting (during production) or the cracks that are formed in service. Continuum mechanics calculations fail to handle such cracks and these are handled by fracture mechanics, which primarily studies resistance of a material to fracture. When a crack is formed, stress is usually not transmitted between the two “virtual” surfaces that form the following crack formation. We can consider a cohesive zone, of zero thickness, to bridge these two virtual surfaces and it can be suggested that fracture happens when material strength in this zone decays. The basic premise for cohesive-zone model is that infinite stresses (stress singularity) at crack tip are unrealistic. Dugdale [40] and Barenblatt [41] developed the first models to address this shortcoming and model brittle fracture. Their models divide the crack tip into stress-free region and cohesively stressed regions. More than a decade later, Hillerborg [42] was the first one to use finite element analysis for cohesive-zone modeling of brittle fracture. Tvergaard’s [43–44] modification to Dugdale and Barenblatt’s model ensured that continuum elements remain intact in the cohesive model. It is their first decohesion model where traction is dependent on both normal and tangential separation, whereas in Barenblatt’s model, it is dependent on the distance from the crack tip. They do this by defining the cohesive elements between the continuum elements such that damage initiation is characterized by opening of the cohesive elements (and governed by some material law) and at failure these elements lose their stiffness, resulting in the continuum elements disjointing from each other. This model is also called the cohesive-zone model or damage zone model (DZM). The initiation and evolution of the damage are defined according to material-softening law/cohesive law/bridging law/traction-separation law (TSL). TSL can be defined based on the properties of the cohesive zone. Traction-separation laws are not universal, and these are not physical laws to characterize fracture at bonded interfaces. Rather, they are highly material dependent. We can define the TSLs in many ways but all separation laws have one thing in common: they all suggest that the traction initially increases with separation and beyond a critical separation distance it diminishes to zero (material softening) at which point these cohesive surfaces detach from each other and a real crack tip is defined.

For implementing TSL in CZM, two constitutive relations are considered: one for the continuum or the “bulk behavior” of the blocks and another for the gap between the surfaces that is defined by the traction-displacement relationships. If the cohesive zone is a thin adhesive layer bonding the blocks, we can define the TSL as the traction at this interface with respect to the separation at the interface. In the modeling by commercial finite element software ABAQUS, the cohesive zone is represented by a single layer of elements that connect the two blocks/surfaces and is constrained by one of the blocks as can be seen from **Figure 1**. Each element on the surface of one block is “tied” with its corresponding neighboring element in the other block through a series of “inelastic springs” which govern the interaction between these two surfaces. This tie constraint prevents the propagation of singular modes in deformation that can arise due to lack of membrane stiffness. The debonding/crack-opening process is akin to the failure of these springs.

### 2.1. Cohesive zone implementation in ABAQUS with an example problem

To show the implementation of cohesive zone model with traction-separation law in ABAQUS, we have provided an example problem [43] below. In this, we study the propagation of a crack in a material block. A debonding model is taken from Tvergaard’s work [43] and is described in brief below. This model does not consist of a debonding potential, but is designed to address the effect of fiber debonding from a metal-matrix composite on the tensile properties (as this requires one to consider both normal and tangential separations in the debonding). For this, first, a function *F*(*λ*) was chosen such that

The nondimensional parameter *λ* is defined as

For a monotonically increasing value of

ABAQUS has vast libraries of constitutive material models but this library is neither a complete one nor a completely flexible one. However, fortunately, user-defined constitutive material models can be implemented in ABAQUS through either UMAT or VUMAT. UMAT is an implicit time integration method (used with ABAQUS/Standard) and it requires the material stiffness matrix as input for forming the consistent Jacobian (*C*), which is also called algorithmic tangential stiffness or ATS, of the nonlinear equations. For linearly elastic materials, it is defined as

where

VUMAT, on the other hand, uses explicit time integration and hence does not require Jacobian matrix. It is implemented in ABAQUS/Explicit and is easier to write. In the case of cohesive zone modeling, we can include element failure and progressive damage failure in the VUMAT subroutine. The VUMAT subroutine consists of the strain decomposition of the total strain into elastic (computed from a linear elastic or hyperelastic constitutive model) and plastic components and flow rule to define the material’s plastic behavior. A sample of the input file (filename.inp) and VUMAT subroutine (vumat.f) are annotated and provided at the end of this section. These are run in ABAQUS through the following command:

abaqus job=*filename.inp* user=*vumat.f*

There are a minimum of five parameters that have to be prescribed for the CZM. They are as follows:

Cohesive strength (

*T*_{n}and*T*_{t}, the peak values of the traction-separation curve).Cohesive energy (area under the traction-separation curve).

Characteristic length (

*u*_{n}and*u*_{t}, usually the separation value corresponding to the cohesive strength).

Cohesive strength (related to yield stress) and cohesive energy signify the maximum resistance to fracture and dissipation following material separation, respectively. For the cohesive zone modeling with traction-separation law, we need to specify the initial elastic stiffness (*K*).

We cannot use the exact value of the material’s Young’s modulus (*E*) for this as this stiffness here refers to a “penalty stiffness,” that is, we are free to choose a value that is necessary to prevent interpenetration of the surfaces and for the convergence of our model. One rule of thumb is to take a value of *K* that is equal to *E/h* (*h* is the cohesive layer thickness). Hence, *K* can be taken to be a very large value as the cohesive elements are considered to be a layer of zero thickness. Also, for simplification, identical values of *K* are taken for the crack-opening mode (mode I) and crack-shearing modes (mode II). The chosen damage initiation and damage evolution criteria affect how these two modes of crack opening interact with each other.

The variables we output (stress, strain, etc.) depend on the kind of problems being modeled by the cohesive elements. This is indicated in the .inp file itself in the definition of the section properties by specifying the response type required. The stress components are as follows:

During a purely normal crack opening (**Figures 2** and **3**). Similarly, during a purely shear crack opening (**Figures 4** and **5**). The negative values in **Figure 5** indicate compressive stresses.

### 2.2. Writing ABAQUS .inp file and VUMAT subroutine

#### 2.2.1. Snippet of input file

Snippet of input file (filename.inp) containing the geometry, meshing, material parameters, and output variables

*Heading

** Job name: cohesive-zone Model name: Model-1

** Generated by: Abaqus/CAE 6.9-1

*Preprint, echo=NO, model=NO, history=NO, contact=NO

** PARTS

*Part, name=Part-1

*End Part

**

** ASSEMBLY

**

*Assembly, name=Assembly

**

*Instance, name=Part-1-1, part=Part-1 ****(two unit blocks/surfaces, each with 4 nodes)*

*Node

⋮

*Element, type=CPE4R****(4 noded continuum plane-strain elements for the blocks)*

⋮

*Element, type=COH2D4 ****(4 noded cohesive elements for the cohesive zone)*

⋮

*Nset, nset=_PickedSet_{i}, internal, generate

⋮

*Elset, elset=_PickedSet_{i}, internal

⋮

** Section: Section-1

*Solid Section, elset=_PickedSet_{i}, material=steel

,

** Section: Section-2

*Cohesive Section, elset=_PickedSet_{i}, material=coh, response=TRACTION SEPARATION ****(specifying response as TRACTION SEPARATION to implement TSL)*

*End Instance

**

⋮

*Surface, type=ELEMENT, name=_PickedSurf21, internal

__PickedSurf21_S3, S3

*End Assembly

*Amplitude, name=Amp-1 ****(time dependent loading conditions)*

0., 0., 0.5, 0.5, 1., 1.

*t=0 amp_at_t=0 t=0.5 amp_at_t=0.5 t=1 amp_at_t=1*

**

** MATERIALS

**

*Material, name=coh ****(cohesive zone parameters)*

*User material, Constants = 5 ****(user defined subroutine with 5 parameters)*

200.0,0.25,0.25,1.0,0.5 **,**

*Density

5.4e-7

*Depvar

7, ****(number of solution-dependent variables SDVs)*

*Initial conditions, type=solution

Part-1-1._PickedSet4,0.0,0.0,0,0.0,0.0,0.0,0.0 ****(initial values of the 7 SDVs)*

*Material, name=material-1 ****(material properties of the continuum surfaces)*

*Density

7.8e-06,

*Elastic

2.5e+8, 0.25****(Young’s modulus E and Poisson’s ratio of Material-1)*

** ----------------------------------------------------------------*

** STEP: Step-1

**

*Step, name=Step-1

*Dynamic, Explicit

, 1.

*Bulk Viscosity

0.06, 1.2

** Mass Scaling: Semi-Automatic

** Whole Model

*Fixed Mass Scaling, factor=50.

**

** BOUNDARY CONDITIONS

**

** Name: BC-2 Type: Displacement/Rotation

*Boundary

_PickedSet17, 2, 2

** Name: BC-3 Type: Displacement/Rotation

*Boundary

_PickedSet18, 1, 1

** Name: BC-4 Type: Displacement/Rotation

*Boundary, amplitude=Amp-1

_PickedSet19, 2, 2, 1.

**

** LOADS

**

** Name: Load-1 Type: Pressure

*Dsload, amplitude=Amp-3

_PickedSurf21, P,50.** ***(Surface-2 of and permitted to have maximum displacement 0.5 in x-direction i.e. =0.5)**

**

** OUTPUT REQUESTS

**

*Restart, write, number interval=1, time marks=NO

**

** FIELD OUTPUT: F-Output-2

**

*Output, field,op=new,number interval=50

*Node Output

RF, U ****(output variables at the nodes: reaction forces and displacements)*

*Element Output, directions=YES

S, SDV ****(output variables over the elements: stresses and SDVs which are defined in the .f file of VUMAT)*

**

** FIELD OUTPUT: F-Output-1

**

*Output, field, variable=PRESELECT

**

** HISTORY OUTPUT: H-Output-1

**

*Output, history, variable=PRESELECT

*End Step

#### 2.2.2. Snippet of vumat.f subroutine

The headers of this subroutine are common to every VUMAT and can be obtained from ABAQUS documentation. We have listed below only portions of the user-defined calculations for the properties.

*! Reading material properties specified in abaqus input file through props() array*

*! max nominal stress at which damage is initiated*

* ! max normal displacement at which the material fails and T*_{n}*=0*

*! max tangential displacement at which the material fails and T*_{t}*=0*

* !parameter in *

mu = props(5)

do ii = 1,nblock

** ! Reading the old state variables or SDVs**.

*! normal displacement*

*!tangential displacement*

counter = stateOld(ii,3)

*!tangential traction*

*!normal traction*

** ! Reading the strain Increments**.

** ! In VUMAT shear components are stored as tensor components and not as engineering components**.

** ! so **.

*!Depending on pure normal traction or pure shear traction or mixed mode of debonding, calculate the following:*

*!Increment in stress*

*! stress update*

stressNew(ii,1) = stressOld(ii,1)+

stressNew(ii,2) = stressOld(ii,2)+

*! state variable update*

stateNew(ii,1) = stateOld(ii,1) +

stateNew(ii,2) = stateOld(ii,2) +

stateNew(ii,3) = stateOld(ii,3) + 1

enddo

return

end

## 3. Elastic contact instability through FEM simulations

**Figure 6A** and **B** shows the schematic setup of a soft incompressible initially stress-free elastic film of shear modulus *h*, and lateral length *L*. It is cast over a rough substrate. To demonstrate the effect of roughness here, we have chosen a substrate with a uniform sawtooth profile characterized by its amplitude ^{−3} nm. The mean distance from the free surface of the film and the contactor is termed the gap distance and is represented by the symbol *d* (refer to **Figure 6A**). In simulations, the contactor, as in experiments, is brought near the film surface from a finite distance at an approach velocity **Figure 6B**. The thickness of the film considered for the simulations is 500 nm and length *L* of the film considered is *β* is varied from 0.1 to 0.7 and *n* is varied from 1 to 12 (for *L = 8h*) in the simulations. Though in reality the film is three-dimensional (3D) as shown in **Figure 6A**, for ease of computation the simulations are carried out considering a plane-strain model (as because of normal traction from the contactor the displacements can be considered to be confined in the plane of the paper) and a 2D geometry as shown in **Figure 6B** has been mostly considered. To understand the morphology in greater detail, 3D simulations are also carried out in certain cases with the 3D geometry.

If one considers an infinitesimal control volume inside the elastic film, a force balance will give rise to the following governing thin-film (Navier’s) equation:

where, *p* is the pressure across the film, and **I** is the identity matrix. The second term,

where Ω is the discretized domain inside the film and Φ is the shape function. The desired discretization and subsequent simulations to realize the film deformation and morphology due to influence of the external contactor are carried out with the help of the commercial finite element software ABAQUS

### 3.1. Geometry and mesh creation

To mimic the physical geometry and to simulate several cases of substrates (having varying frequency and amplitudes) and to incorporate different perturbation conditions at the film-contactor interface, the graphical user interface in the software proves to be difficult. Python scripting was used instead to generate the multitude of geometries. For the film’s lateral edges and the substrate boundary, straight lines are used, while for achieving the perturbed top surface cubic splines are fitted to get a continuous line out of the random sinusoidal fluctuations at the film-contactor interface. Meshing forms an integral part of the preprocessing steps, which helps in discretizing the geometry. The accuracy that can be obtained from any finite element analysis is directly related to the finite element mesh that is used. Generally, for the large deformations in 2D, eight node elements are used but those are computationally intensive. In 2D as discussed earlier for the present scenario one can consider a plain-strain condition, therefore, continuum plain-strain 4 node-reduced integration elements (CPE4R) are used in 2D. For 3D geometry, continuum three-dimensional eight-node-reduced integration (C3D8R) elements are used. To circumvent the problem of artificial stiffness due to shear locking that can arise due to four node elements, a large number of elements are used while maintaining the aspect ratio of the elements close to 1, which also helps in avoiding problems associated with skewness. Further, for the large deformation problem, reduced integration elements are shown to be successfully implemented without impacting the results and reduce the computational load. On the top of the film, the force is much higher and decreases drastically with depth in the film due to the highly nonlinear dependence of force on distance. Thus, it becomes necessary to have high density of mesh toward the top of the film. The geometry of the film in case of a smooth substrate is seeded with a bias ratio to have smaller elements toward the top as shown in **Figure 7A** and **B** (for 2D and 3D geometries, respectively). For sawtooth profile, the geometry is divided into two blocks. The top block has similar seeding as in case of smooth substrate and the bottom is meshed as shown in **Figure 7C** and **D** (for 2D and 3D geometries, respectively). Since the force applied depends on the amplitude of substrate roughness as well, the meshed geometry is tested with the grid independence test for every case and the finalized meshes are then chosen.

### 3.2. Boundary conditions

To imitate the rigid boundary conditions of the base in the simulations, the bottom edge of the film is pinned to the substrate and movement toward either direction is restricted with the help of an Encastre boundary condition. To incorporate periodic boundary conditions on the lateral sides of the film, x-symmetric boundary conditions are applied.

### 3.3. Force calculations

As discussed earlier, so far the failure at the interfaces was captured via a cohesive zone model. In the model, the contactor (a rigid body) and the elastic film were considered to be in intimate contact separated only by a very thin layer of cohesive zone. In the cohesive zone, a traction-separation law (inbuilt in ABAQUS) is defined in such a way that the traction initially increases with separation and beyond a critical separation distance it diminishes to zero (material softening) at which point these cohesive surfaces detach from each other. The energy release rate at debonding is given by the area under the curve. The problem with the method is except for the energy all the other parameters are not well characterized and the traction-separation law is not based on any physical law. Also, while the model is good for peeling/debonding experiments and almost always considers the initial stage as one where the rigid contactor and the film is in intimate contact, it fails to capture the essential morphologies at critical separation distance.

Thus, to capture both the adhesion and debonding of an elastic film due to the interaction between the film and the contactor, a body-force approach has been considered, where each node in the film experiences a varying magnitude body force depending on its relative position with the contactor. The most generic interaction between two different materials when it comes in contact proximity is given by van der Waals type of potential. The attractive potential of VDW prevalent between two atoms or molecules separated at a distance *r* can be described as [45]

where *M* is the coefficient of attraction. However, a molecule in the film experiences attraction from all the molecules in the contactor. Thus, if one considers a ring-like element in the contactor as shown in **Figure 8A**, the total number of atoms/molecules present in the ring will be *D* which is given by

If the single molecule of the film is now replaced by a solid body of unit volume having number density of atoms/molecules

where *A* is known as the Hamaker constant and has the value of 10^{−19} J for the present simulations. In the simulations, the above body-force term is modified and a born-repulsion term is also added to avoid extremely high forces at total contact such that

where ^{6} (*l*_{e} is an equilibrium distance and has a value of 0.138 nm for the present case). The nature of the body force can be seen in **Figure 8B**.

The contactor in the simulations is virtual in nature, and its presence is felt only through the body-force term *f*. Since the contactor position keeps on changing with time (refer to **Figure 6B**) so also does the body force. Hence, the gap distance

During adhesion:

During resting:

During debonding phase:

In the above set of equations,*t*,

In VDLOAD subroutine of ABAQUS, the body force for a particular node is thus calculated according to

The morphological evolution of the film due to the presence of the contactor is thus dependant on the dynamic contactor position and thus is a time-dependent problem. The solution is obtained by marching forward in time using an explicit scheme such that the body force calculated at a particular time *t* is given by

### 3.4. Structure of VDLOAD subroutine

The headers of this subroutine are common to every VDLOAD and can be obtained from ABAQUS documentation. We have listed below only portions of the user-defined load subroutine

***** Declare new variables used (if any)**

***** Declare the values of contactor approach velocity, retracting velocity, van der Waals potential/force parameters, initial separation distance, cutoff distance**

do k = 1, nblock

u2= curCoords (k,2) ***** curCoords(k,1) and curCoords(k,2) represent the updated coordinates in x**

_{1}

**and**

*x*_{2}

**direction respectively)**

**** Calculate gap distance according to Eq. 10**

***** Calculate forces according to Eq. 12**

**** Calculate force update according to Eq. 13 and input the magnitude in** the array value (*k*), which is the prescribed **ABAQUS declared array to contain forces**

end do

The FE results displaying the effect of the substrate roughness on the morphologies formed, the work of adhesion exhibited, are discussed next in “Results” section.

### 3.5. Results

Linear stability analysis performed with the help of a single Fourier mode [16, 17] and energy-minimizing nonlinear analysis [18–20] has already been able to identify the nontrivial bifurcation modes for which instability is possible. For critical conditions at which instability first initiates, it has been observed that the critical wave length *j*. When the particular geometries are simulated using VDLOAD, only in the adhesive branch, the displacements at the corner-most node, which has the highest amplitude (because of the cosine nature of the waves), is tracked. The displacement of the node thus obtained at any generic time near critical separation distance is plotted against the average wave modes in **Figure 9A**. It can be seen that the curve has a bell shape indicating that near a particular mean surface wave mode the deformation of the surface is highest. This marks the dominant wave mode, whose value matches with those from Linear stability analysis (LSA) and experiments involving flat substrates.

To study whether the substrate roughness has any influence on the emerging surface patterns, similar studies have been conducted for sawtooth-profiled substrates having different substrate parameters of *n*. The results reveal that though the substrate wavelengths have negligible effect, the substrate amplitudes have a profound influence on the emerging wavelengths. The surface wave mode/wave lengths are found to be an increasing/decreasing function of the substrate amplitude averaged over the different substrate wave modes considered, indicating that the patterned substrates do engender miniaturized length scales as can be seen from **Figure 9B**.

The displacement profile of the corner node at different time is shown in **Figure 10**. It can be seen that at initial time *t* = 0 since the initial gap is much greater than the critical value, the film surface remains planar with zero displacement. Near critical distance (near

In ABAQUS, not only a single node but the whole evolution of the film can be traced throughout adhesion-debonding. In **Figure 11**, the labels A1-E1 indicate the morphologies and A2-E2 the stress profiles that develop across a film cast on a sawtooth-profiled substrate (of amplitude **Figure 11** A2-E2. It can be seen that the stresses lie concentrated at the column edges. Thus, more number or greater width of columns (B2 and C2) are seen to have higher stresses at the edges which are responsible for subsequent peeling and stress release as seen in **Figure 11** D2 and E2.

The sum of the stresses across the total top surface area is a signature of the forces required by the contactor to pull the film up to a certain pull-off distance. The area under the force displacement curve from the beginning of debonding phase to snap-off yields the work of adhesion, which is required to separate the two surfaces of the film and the contactor. **Figure 12** reveals the work of adhesion exhibited by films cast on substrates having various degrees of roughness. It can be seen vividly that while the effect of the wavelength of the substrate is negligible, the higher the substrate amplitude, more is the substrate roughness more is the work of adhesion and better is the performance of the elastic film as an adherent. The work of adhesion is also seen to be a function of the debonding velocity. For a higher debonding velocity (**Figure 12B**), it can be seen that the work of adhesion is also high (compared to **Figure 12A**). The fact is intuitive considering in band-aids when we pull at higher velocity higher adhesiveness is felt compared to that felt at lower velocities [46].

## 4. Conclusions

Opening of crack surfaces through traction-separation laws in cohesive-zone modeling is demonstrated for cohesive blocks. In the case of adhesion-debonding of elastic films, simulations are performed for topologically rough substrates to demonstrate the ability of the method to resolve complex geometries, which these analytical methods are incapable of. At incipience of instability when the film surface is perturbed by random fluctuations at the vicinity of prescribed wavelengths, the displacements made by those contacting surface nodes are found highest near the dominant wavelength. These dominant perturbations grow into columnar structures at contact. A thorough FE study reveals that film cast on rough substrate is ideal for fabrication of features/columns at miniaturized length scales and exhibits high work of adhesion. Thus, these films can act as better adhesives and can produce enhanced functional surfaces.