Simulating Contact Instability in Soft Thin Films through Finite Element Techniques

Abstract

When a thin film of soft elastic material comes in contact with an external surface, contact instability triggered by interaction forces, such as van der Waals, engenders topologically functionalized surfaces. Innumerable technological applications such as adhesives; microelecromechanical systems (MEMS), and nanoelectromechanical systems (NEMS) demand understanding of the physics behind the mechanical contact, relationship between the morphologies, and detachment forces in such films. Indentation tests are important experimental approach toward this; there also exist many simulation procedures to model the mechanical contact. Both atomistic level and analytical continuum simulations are computationally expensive and are restricted by the domain geometries that can be handled by them. Polymeric films also particularly demonstrate a rich variety of nonlinear behavior that cannot be adequately captured by the aforementioned methods. In this chapter we show how finite element techniques can be utilized in crack opening and in contact-instability problems.

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.

Figure 1.

Meshing with tie constraints for cohesive elements in ABAQUS.

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

$F\left(\lambda \right)=\frac{27}{4}{\sigma }_{\text{max}}(1-2\lambda +{\lambda }^2) \text{and} 0\le \lambda \le 1$E1

The nondimensional parameter λ is defined as

$\lambda =\sqrt{{\left(\frac{u_n}{{\delta }_n}\right)}^2+{\left(\frac{u_t}{{\delta }_t}\right)}^2}$E2

For a monotonically increasing value of $\lambda$ the tractions are $T_n=\frac{u_n}{{\delta }_n}F\left(\lambda \right)$ and $T_t=\alpha \frac{u_t}{{\delta }_t}F\left(\lambda \right)$. Here, the normal and tangential tractions are $T_n$ and $T_t$respectively, normal and tangential displacements are $u_n$and $u_{t }$respectively, ${\delta }_n$ and ${\delta }_t$ are the maximum separation distances in normal and traction modes, respectively. Here, we have considered that the two virtual surfaces that form after a crack growth are of unit dimension. For solving this, we have developed our own code using MATLAB and ABAQUS. MATLAB is excellent to call external commands to run the PYTHON scripts for generating the input files for ABAQUS.

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

$C=\frac{\partial \Delta \sigma }{\partial \Delta \epsilon }$E3

where $\Delta \sigma$ is the increment in Cauchy stress and $\Delta \epsilon$ is the increment in strain.

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:

1. Cohesive strength (T_n and T_t, the peak values of the traction-separation curve).

2. Cohesive energy (area under the traction-separation curve).

3. 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).

Figure 2.

Cohesive zone model (CZM) in ABAQUS. Loads and boundary conditions for the case of purely normal traction (T_t = 0).

Figure 3.

Normal traction (T_n) versus separation (u_n) obtained after CZM in ABAQUS.

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:

1. ${\sigma }_{11}$—Direct membrane stress

2. ${\sigma }_{22}$—Direct through-thickness stress

3. ${\sigma }_{12}$—Transverse shear stress

Figure 4.

Cohesive zone model (CZM) in ABAQUS. Loads and boundary conditions shown for the case of purely shear traction (T_n = 0).

Figure 5.

Tangential traction (T_t) versus shear separation (u_t) obtained from CZM in ABAQUS.

During a purely normal crack opening ($u_t=0$), according to our TSL, the normal traction $T_n$ initially increases with the normal separation $u_n$, reaches a maximum value ${\sigma }_{\text{max}}$ and reduces to 0 when $u_n={\delta }_n$ that is, when the final separation occurs (see Figures 2 and 3). Similarly, during a purely shear crack opening ($u_n=0$), the shear traction $T_t$ initially increases with tangential separation $u_t$, reaches a maximum value $\alpha {\sigma }_{\text{max}}$ and then at $u_t={\delta }_t$ the final separation occurs (see 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 $***({\sigma }_{max},{\delta }_n\text{, }{\delta }_t$,$\alpha .$ these are read through props( ) matrix in vumat.f)

*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

${\sigma }_{max}$= props(1) ! max nominal stress at which damage is initiated

${\delta }_n$= props(2) ! max normal displacement at which the material fails and T_n=0

${\delta }_t$ = props(3) ! max tangential displacement at which the material fails and T_t=0

$\alpha$ = props(4) !parameter in $T_t=\alpha \frac{u_t}{{\delta }_t}F\left(\lambda \right)$

mu = props(5)

do ii = 1,nblock

! Reading the old state variables or SDVs.

$u_n$ = stateOld(ii,1) ! normal displacement

$u_n$ = stateOld(ii,2) !tangential displacement

counter = stateOld(ii,3)

$T_t$ = stateOld(ii,4)!tangential traction

$T_n$ = stateOld(ii,5)!normal traction

! Reading the strain Increments.

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

! so $\text{Δ}{u}_t$ will be twice of tensorial shear strain.

$\text{Δ}{u}_n$ = strainInc(ii,1)

$\text{Δ}{u}_t$ = 2.0*strainInc(ii,2)

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

$\begin{array}{lll}\lambda \hfill & =\hfill & \cdots \cdots \cdots \hfill \\ T_{n,\text{new}}\hfill & =\hfill & \frac{u_n}{{\delta }_n}F\left(\lambda \right)\hfill \\ T_{t,\text{new}}\hfill & =\hfill & \alpha \frac{u_t}{{\delta }_t}F\left(\lambda \right)\hfill \end{array}$E102

$\begin{array}{l}{\sigma }_{11}=\left(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{${\delta }_n$}\right.\right)\left(F\left(\lambda \right)+u_n\frac{\partial F\left(\lambda \right)}{\partial u_n}\right); \frac{\partial F\left(\lambda \right)}{\partial u_n}=\frac{\partial F\left(\lambda \right)}{\partial \lambda }\frac{\partial \lambda }{\partial u_n}\hfill \\ {\sigma }_{12}=\left(\raisebox{1ex}{$u_n$}\!\left/ \!\raisebox{-1ex}{${\delta }_n$}\right.\right)\frac{\partial F\left(\lambda \right)}{\partial u_t}; \frac{\partial F\left(\lambda \right)}{\partial u_t}=\frac{\partial F\left(\lambda \right)}{\partial \lambda }\frac{\partial \lambda }{\partial u_t}\hfill \\ \begin{array}{l}{\sigma }_{13}=\alpha \left(\raisebox{1ex}{$u_t$}\!\left/ \!\raisebox{-1ex}{$t$}\right.\right)\frac{\partial F\left(\lambda \right)}{\partial u_n}\\ {\sigma }_{14}=\left(\raisebox{1ex}{$\alpha $}\!\left/ \!\raisebox{-1ex}{${\delta }_t$}\right.\right)\left(F\left(\lambda \right)+u_t\frac{\partial F\left(\lambda \right)}{\partial u_t}\right)\end{array}\hfill \end{array}$E101

!Increment in stress

$\begin{array}{c}\text{Δ}{T}_n=({\sigma }_{11}*\text{Δ}{u}_n)+({\sigma }_{12}*\text{Δ}{u}_t)\\ \text{Δ}{T}_t=({\sigma }_{21}*\text{Δ}{u}_n)+({\sigma }_{22}*\text{Δ}{u}_t)\end{array}$E100

! 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

stateNew(ii,4) =

stateNew(ii,5) =

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 $\mu$<1 MPa, mean thickness 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 $\beta h$ and wavelength ${\lambda }_s=L/n$ where $n$ is the number of substrate troughs or crests. In the asymptotic case when $\beta h\to 0$ or $L/n\to \infty$, we get the classical case of a smooth substrate, though in real systems it is almost impossible to come across a substrate that is atomistically smooth and almost all surfaces have inherent roughness. Again in real situations mechanical/thermal noises and vibrations almost always disturb the film surface at the air interface. To incorporate them for numerical simulations, the film surface is perturbed by sinusoidal disturbances having different frequencies and random amplitudes of the order of 10⁻³ 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 $v_a (~1 nm/s)$ maintained below the critical distance where the instabilities first initiate (adhesion) for around ~15 s. This is necessary to allow the film to have sufficient time to arrange itself under the interplay of imposed force and elasticity of the film. Once the patterns get mature, the debonding phase is started. In the debonding phase, the contactor is pulled off at different velocities (debonding: $v_d$). The contactor positions at different times can be visualized from Figure 6B. The thickness of the film considered for the simulations is 500 nm and length L of the film considered is $8h$. The β 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.

Figure 6.

Schematic diagram of an elastic thin film case on a (A) 3D sawtooth and (B) 2D sawtooth-patterned substrate and in contact proximity to an external contactor. The position of the contactor at different time is also shown.

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:

${\displaystyle \overrightarrow{\nabla }}\cdot {\displaystyle \overline{{\displaystyle \overline{\sigma }}}}+{\displaystyle \overrightarrow{f}}=0$E4

where, ${\displaystyle \overline{{\displaystyle \overline{\sigma }}}}=pI+\mu \left({\displaystyle \overrightarrow{\nabla }}{\displaystyle \overrightarrow{u}}+{\displaystyle \overrightarrow{\nabla }}{{\displaystyle \overrightarrow{u}}}^T\right)$ is the stress developed inside the film, ${\displaystyle \overrightarrow{u}}$ is the displacement vector, p is the pressure across the film, and I is the identity matrix. The second term, ${\displaystyle \overrightarrow{ f,}}$ represents the body force (per unit volume) that the film experiences because of the contactor. The equation as can be seen is a partial differential equation in the unknown ${\displaystyle \overrightarrow{u}}$ and does not yield an analytical solution in the case when the substrate is patterned or when the film material displays geometric or material nonlinearity. A numerical solution can, however, be obtained when the differential equation is discretized by considering a trial/shape function for the displacement and converted into a set of algebraic equations. In the Galerkin finite element formulation, the residual obtained because of the trial function can now be minimized in a weighted fashion according to

${\displaystyle {\displaystyle \underset{\Omega }{\int }}\left({\displaystyle \overrightarrow{\nabla }}.{\displaystyle \overline{{\displaystyle \overline{\sigma }}}}+{\displaystyle \overrightarrow{f}}\right)}\varphi d\Omega =0$E5

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.

Figure 7.

The meshing done inside the film geometry for (A)&(B) a planar substrate, (C)&(D) for a sawtooth substrate, respectively. Figures on the left and right represent meshing in 2D and 3D respectively.

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]

$\text{Ψ}\left(r\right)=-M/r^6$E6

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 $2\pi {\rho }_{c}d{x}_{1}d{x}_2$ where ${\rho }_c$ denotes the number density of atoms/molecules of the contactor. Hence, the total effective potential that a molecule in the film will experience due to a contactor that is semi-infinitely thick is laterally unbounded and is situated at a surface-to-surface separation distance of D which is given by

$w\left(D\right)=-2\pi {\rho }_{c}M{\displaystyle {\int }_{x_2=D}^{x_2=\propto }d{x}_2}{\displaystyle {\int }_{x_1=0}^{x_1=\propto }\frac{x_{1}d{x}_1}{{\left(\sqrt{x_1^2+x_2^2}\right)}^6}=-\frac{\pi M{\rho }_c}{6D^3}}$E7

If the single molecule of the film is now replaced by a solid body of unit volume having number density of atoms/molecules ${\rho }_f$, the body force that the elastic film will experience is given by

$f\left(D\right)=-{\rho }_f\frac{dw}{dD}=-\frac{\pi M{\rho }_f{\rho }_c}{2D^4}=-\frac{A}{2\pi D^4}$E8

where A is known as the Hamaker constant and has the value of 10⁻¹⁹ 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

$f\left(D\right)=-\frac{A}{2\pi D^4}+\frac{9B}{D^{10}}$E9

where $B=A{l}_e{}^6/18\pi$ Jm⁶ (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.

Figure 8.

(A) Schematic diagram of the calculation of force between a contactor surface and a molecule inside the elastic film. (B) It shows the van der Waals body force present between the contacting surface and any nodal point in the mesh inside the film geometry.

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 $\xi$ between a particular node and the surface of the contactor keeps on changing with time as per the following:

During adhesion: $\xi =\left(D_{\text{initial}}-v_a\times t\right)-u_2\left(t\right)$ for $0\le \text{t}\le t_1\text{ where }{t}_1=(d_c+1)/v_a$

During resting:$\xi =\left(D_{\text{initial}}-\left(d_c+1\right)\right)-u_2\left(t\right)$ for $t_1\le \text{t}\le t_2\text{ where} t_2(d_c+1)/v_a+15$

During debonding phase:

$\xi =\left(D_{\text{initial}}-\left(d_c+1\right)\right)+v_d\times \left(t-t_2\right) -u_2\left(t\right)For{t}_2 \le t$E10

In the above set of equations,$D_{\mathrm{initial}}$ denotes the initial gap distance, $u_2\left(t\right)$ the normal position of the particular node at time t, $d_c$ the critical distance at which the film begins to experience substantial force due to the contactor and starts deforming. Its value can be roughly estimated analytically from linear stability analysis on sinusoidal substrates according to [30]

$\frac{hA}{2\pi \mu d_c^4}=\frac{6.22}{\left(1-\beta \right)}$E11

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

$f\left(\xi \right)=-\frac{A}{2\pi {\xi }^4}+\frac{9B}{{\xi }^{10}}$E12

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

$f\left(t\right)=f\left(t-\text{Δ}t\right)+\frac{\partial \xi }{\partial t}\times {\frac{\partial f}{\partial \xi }|}_{\left(t-\text{Δ}t\right)}$E13

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₁ and x₂ 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 $\text{ λ }~2.96\text{ h }$ is at par with the experimental observations. In FEM to obtain the critical wave mode, the geometry is created with the help of a series of cosine waves having random amplitudes and frequencies in the vicinity of a mean value 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.

Figure 9.

(A) The displacement of the top-most corner node of the film for different perturbed film surfaces at a particular time near critical. The dominant wave mode is ~3 for flat substrate. (B) The dominant wave mode and the wavelength for different substrate amplitude averaged over the substrate wavelength. With increased roughness in substrate miniaturization of the surface patterns occurs.

Figure 10.

The trajectory of the corner-most point with time.

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 $\beta$ and 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.

Figure 11.

(A1-E1) Morphologies at different contactor positions. (A2-E2) Stresses in the film at different contactor positions for a film cast on a sawtooth-profiled substrate of amplitude $\beta =0.1$ and $n=4.$

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 $t_1$), the displacements start to have nonzero values and at $t$ just greater than $t_1$ when the contactor is in resting position the film comes in total contact. The node continues to be in contact with the contactor even in the debonding phase when the contactor is pulled up and it gets detached only at snap-off distance ${\text{d}}_{\text{snap}}$. The trajectory of the corner node arising due to different debonding velocity is found to be different in the debonding phase. It can be seen that the snap-off also occurs at different times. Though the snap-off time ${\text{t}}_{\text{snap}}$ for the lower velocity is found to be highest, the snap-off distance ${\text{d}}_{\text{snap}}={\text{v}}_{\text{d}}\times {\text{t}}_{\text{snap}}$ is in reverse order $\left(~30,16,5,10,6.6\text{ nm}\right)$ showing that in case of lower debonding velocity the detachment occurs at the lowest snap-off distance.

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 $\beta =0.1$ and $n=4$) for the different contactor positions. In the inset, a contactor displacement graph is provided for easy understanding of the positions. At position A1 where the contactor is just near critical distance the surface is seen to begin destabilizing where the structures formed have very less asperities. At stage B1, the film surface comes in physical contact with the contactor and the top surface takes the shape of columns bridging the gap with intervening cavities (which form the dispersed phase). At C1, the end of resting phase, the columnar structures mature to form continuous ridges, which upon withdrawal to D1 shrinks in size due to peeling from the edges. Continuous peeling ultimately gives rise to isolated columns (which now forming the dispersed phase) as shown in E1 leading to snap-off. Thus, from adhesion to debonding, the morphologies at the surface go through a phase inversion. The corresponding stress profiles can be seen from 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.

Figure 12.

Contour plots of the work of adhesion (J/m²) for different substrate parameters for a film cast on a rough (sawtooth-profiled) substrate. (A) For debonding velocity $v_d=0.125 nm/s$. (B) For debonding velocity, $v_d=0.25 nm/s$.

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.

Acknowledgments

Funding from the Department of Science & Technology, New Delhi, India, is highly acknowledged.