Detailed modelling of delamination buckling of thin films under global tension.

Tensile specimens of metal films on compliant substrates are widely used for determining interfacial properties. These properties are identified by the comparison of experimentally observed delamination buckling and a mathematical model which contains the interface properties as parameters. The current two-dimensional models for delamination buckling are not able to capture the complex stress and deformation states arising in the considered uniaxial tension test in a satisfying way. Therefore, three-dimensional models are developed in a multi-scale approach. It is shown that, for the considered uniaxial tension test, the buckling and associated delamination process are initiated and driven by interfacial shear in addition to compressive stresses in the film. The proposed model is able to reproduce all important experimentally observed phenomena, like cracking stress of the film, film strip curvature and formation of triangular buckles. Combined with experimental data, the developed computational model is found to be effective in determining interface strength properties.

Introduction

The interface properties of thin, brittle films on compliant substrates are of great interest in applications such as hard protective coatings, gas barrier coatings, thin film transistors, flexible electronics and sensors. In these applications, interfacial adhesion dominates the mechanical performance where high adhesion is desired over poor adhesion. While there are several techniques available to measure the interfacial properties of thin films on rigid substrates, including nanoindentation [1,2], four-point bending [3] and spontaneous buckle delamination [4], it can be difficult to apply the same techniques to compliant polymer substrate systems [5]. Currently, the most commonly used techniques to assess interface properties of a thin film on a compliant substrate is to use bending experiments [6] or uniaxial tensile straining to induce film fracture and delamination in the form of buckles [7-9]. The tensile straining technique is based on the fragmentation test and described by shear lag theory [10]. Both adhesion measurement approaches (bending and tensile straining) rely on accurate mathematical models of the experimental process to determine the interface properties.

Many theoretical models have been developed over the years to describe brittle thin film fracture [4,11-15]. These models have been able to explain the failure criteria of thin films on a range of substrates due to thermal or mechanical strain. Another set of modelling efforts has focused on the buckling phenomena of films on rigid, i.e. very stiff, substrates [4,16,17]. The most successful theory is that of Hutchinson and Suo [4] to measure the interfacial adhesion using only the dimensions of the buckle and the elastic properties of the film. In this case, the buckles form spontaneously due to high residual stresses in the film or can be triggered by indentation. With respect to films on elastomer substrates, the wrinkling of films [18-20] is of great interest due to the possibility of providing more “stretchability” to flexible electronics [21]. The combination of models to take into account the fracture and buckling of brittle thin films on compliant polymer substrates under tensile strain and to help describe interfacial properties, as in this work, is relatively new. The existing approaches are based on two-dimensional (2-D) models [7] or shear-lag theory [9], and are thus limited in applicability. Particularly in the case of the uniaxial tensile test, as considered in this paper, 3-D effects play a dominant role that never can be achieved by 2-D models.

An analytical model based on an energy balance is able to relate measured buckle geometries from a tensile straining experiment to the adhesion energy of the interface [7]. However, this model was developed as a 2-D model and does not accurately capture the real interface failure or, consequently, buckling behaviour. Experiments have shown [7,8] that buckling and the associated delamination process are initiated and driven by interfacial shear in addition to compressive stresses in the film. Motivated by the deficiencies of the current 2-D models, the goal of this paper is to demonstrate a computational multi-scale 3-D finite element simulation of the experimental process. This new approach uses a two-stage model, which is able to capture all of the important experimentally observed effects, i.e. correct load transfer between film and substrate before and after cracking of the film, shear-stress-induced local out-of-plane deformation of the film and mixed mode delamination. Careful consideration of these 3-D effects is required for a sufficiently accurate determination of initial fracture stress and strain, initial buckling strain, etc. A macromodel comprising the tensile test specimen is analysed first, delivering boundary conditions for the micromodel of a single localized buckle. The modelling strategy involves cohesive zone elements to model the interface between the film and the substrate. The determination of interface properties, i.e. parameters of the cohesive elements, from the finite element model together with experiments is therefore an inverse problem that can be solved by conducting a parameter study.

Finite element modelling of the experimental process

In the uniaxial tension tests performed by Cordill et al. [7,22], specimens of a 100 nm chromium film, subsequently called “film”, on a 50 μm polyimide substrate, subsequently called “substrate”, with effective dimensions of 10 mm times 20 mm, were strained with a constant strain rate in a low-load tensile straining device (Kammrath & Weiss, Dortmund, Germany) and observed inside a scanning electron microscope. Channel cracks transverse to the straining direction initiate at a global strain of 0.4% and, with increasing strain, cracks continue to form until saturation spacing is obtained. The saturation crack spacing is 2.8 ± 0.9 μm (Fig. 1a). The different elastic properties of the two laminate layers lead to an observable curvature of the specimen. Additionally, due to the mismatch of the Poisson’s ratio between the film (ν = 0.21) and substrate (ν = 0.34), uniaxial straining leads to compressive stresses transverse to the tensile direction in the film. At a global strain level of 10%, where both materials are in the plastic regime, buckles form in the film strips, causing interfacial failure between the film and the substrate. Buckles generally form in two shapes, rectangular (Fig. 1b) and triangular (Fig. 1c). The buckles will also form cracks at their apex (the top of the buckle) due to the brittle nature of the film and tensile bending stresses. This study will only address non-cracked buckles as the buckles initially form without cracking.

Fig. 1

(a) Example of the cracks that form perpendicular to the tensile direction and rectangular and triangular buckles which form parallel to the tensile direction; (b) SEM micrograph of a rectangular buckle, cracked at the top; (c) atomic force microscope height image of a triangular buckle of a 500 nm thick film.

Macromodel

Macromodels are used to determine the cracking stress in the film and to gain boundary conditions for a micromodel. For these purposes, two different macromodels are required. In both models shell elements are used to represent a quarter of the symmetric tension test specimen (see Fig. 2). Symmetry boundary conditions are used on the edges in the symmetry planes. The loading edge is clamped and translated in the x-direction, modelling the global straining of the test specimen. The analysis is displacement controlled and starts from an undeformed, stress-free initial state, not taking into account any residual stresses, which might be present in the specimen due to the deposition process.

Fig. 2

Schematic sketch of the macromodel (b = 10 mm, l = 20 mm, t = 50 μm, h = 100 nm).

The first macromodel is used to determine the cracking stress of the film, which is the stress state at the global strain level at which cracking is experimentally observed, being 0.4%. At this strain level both materials, substrate and film, are assumed to be in the elastic regime. The specimen is modelled by composite shell elements with bilinear interpolation functions and a total of 11 integration points over the thickness; one integration point represents the film material and the rest represents the substrate. The determined maximum principal stress in the film at the global strain level of 0.4% is the determined cracking stress.

After the film has cracked transversely to the global straining direction, it loses its stiffness in this direction. Therefore, in a second macromodel, the film strips separated from each other by the cracks perpendicular to the global straining direction can be modelled by truss elements on the shell node lines in the y-direction. Truss and shell elements use shared nodes located on the shell reference surface. To accurately model the bending stiffness of the composite, an offset of (h + t)/2 is used for the truss elements (see Fig. 2). Comparisons of the results achieved from this model with the results obtained from a model, in which the cracked film is modelled as an orthotropic shell with no stiffness in the x-direction, have confirmed the correctness of the “truss” model. The analysis is started from an undeformed, initially stress-free state. This is possible because non-linearities do not occur before the film cracks. Contrary to the first macromodel, plasticity needs to be taken into account now, both in the substrate and the film. The uniaxial stress–strain curve of the substrate is experimentally determined. The material behaviour of the film is assumed to be linear elastic, ideal plastic, with a yield stress slightly higher than the previously determined cracking stress. Stress–strain curves for both materials are depicted in Fig. 3. The local strains at the centre of the specimen, where the compressive stresses in the film are at a maximum, are used as boundary conditions for the micromodel.

Fig. 3

Stress vs. strain curves for the modelled materials: (a) polyimide substrate (Young’s modulus E = 3.704 GPa, Poisson’s ratio ν = 0.34, yield stress σ = 34.3 MPa) and (b) chromium film (E = 280 GPa, ν = 0.21, σ = 1200 MPa).

Micromodel

The micromodel is a 3-D continuum model of a small part of a film strip with a portion of the underlying substrate (see Fig. 4). Due to the assumed double symmetry, only a quarter of the detail is modelled. The detail is cut out from the centre of the tensile test specimen. The size of the micromodel should be sufficiently large to capture the local phenomena, and as small as possible to make parameter studies feasible. The width, d, of the model corresponds to half of the average crack spacing obtained from experiments. The portion of the thickness of the substrate to be modelled is determined by the influence depth of film wrinkles. This is assumed to be d = 1.8 μm, which is much smaller than the complete thickness of the substrate (t = 50 μm).

Fig. 4

Schematic sketch of the micromodel (d = 2.8 μm, h = 100 nm).

The boundary conditions (depicted in Fig. 4) are generated from the macromodel with the cracked film. This simulation delivers a relation between the applied global strain and the local strain state at an arbitrary position in the specimen. The micromodel is extremely small in comparison with the length scale of the macroscopic strain variation, and is located in the centre of the specimen, where the strain variations are smallest. Therefore, spatially constant displacements are prescribed at the faces x = d/2 and y = 4d. The values of these displacements are computed from the macroscopic strains and applied in an incremental–iterative way. At the position of the micromodel the macroscopic strain state is dominated by the membrane strains; the rotations are very small (see Fig. 8) and are therefore ignored in the computation of the local displacement boundary conditions on the y = 4d and x = d/2 planes. To ensure continuity, film and substrate displacements in the y-direction must be equal. Global shear deformations, in plane as well as out of plane, are very small at the position of the micromodel and can be ignored with good accuracy. The nodes on the plane z = −d are constrained to this plane, because disturbances to the trivial strain field in the vertical direction arising due to buckle formation are assumed to be negligible at this distance. This assumption could be confirmed by checking the reaction forces that arise. All other z-displacements in the micromodel are left free, hence applying the z-displacement obtained from the macromodel would only lead to a rigid body motion of the whole micromodel. Thus, the computed z-displacements in the micromodel can be regarded as relative displacements with respect to the z = −d plane.

Fig. 8

Global z-displacement in the symmetry planes of the specimen at different global strain levels. The specimen is clamped at x/l = ±1/2.

The substrate and film are connected by cohesive zone elements with zero thickness in the undeformed configuration. The used elements are based on the concept by Tvergaard and Hutchinson [23]. The in-plane separations, δ1 and δ2, as well as the out-of-plane separation, δ3, are evaluated by calculating the separation vector spanned by the midpoints of the lower and upper sides of the deformed cohesive element. A non-dimensional effective separation is then defined aswhere and are the in-plane and out-of-plane separations, respectively, at which complete interfacial failure would occur if either tangential or normal traction were applied separately. The McCauley bracket is used to distinguish between tension (δ3 ⩾ 0) and compression (δ3 < 0). For δ3 < 0, the expression 〈δ3〉 is zero. A functional relationship between generalized traction and effective separation, t(λ) – the so-called traction-separation law, depicted in Fig. 5 – is used. The relation between the non-dimensional separation and generalized traction is assumed to be bilinear, with an ultimate generalized traction of . The rising slope of the traction–separation relation is kept relatively stiff by setting λ1 = λ2 = 0.01. The non-dimensional separation at which the interface fails is at λ = 1. Increasing the parameter λ1 one order of magnitude to 0.1 and to 0.5, respectively, yields only marginal changes in the simulation results. This insensitivity was also observed in a related problem [24]. The traction potential, defined asis used to calculate the traction vector by t = ∂Φ/∂δ, which can be expressed in matrix formIn the case of compression, the normal traction is calculated byusing a penetration stiffness multiplier p = 10. This means that the initial stiffness of the interface in compression is ten times larger than in tension, effectively preventing any substantial material penetration. The interface toughness per unit area, expressed in terms of an energy release rate, is the area under the traction–separation curve. For the assumed bilinear traction separation law, the energy release rates for mode I, , and for mode II, , are related to the ultimate generalized traction and the ultimate separations byThe properties of the cohesive zone are mainly governed by the parameters and .

Fig. 5

General traction–separation law for the cohesive zone model.

In the current study, the results achieved with three different values, 40 GPa, 80 GPa,and 120 GPa, for the ultimate generalized traction, , as well as various energy release rates, were compared against each other. For the energy release rate in mode I, , the values 0.5 Nm/m2, 1.0 Nm/m2, 4.5 Nm/m2 and 10.0 Nm/m2 were taken. The energy release rate in mode II was set by specifying the ratio . For this, values for v of 0.1, 0.5, 1, 2 and 10 were taken. An important impact of the values of v on the buckling and delamination process was observed.

In order to induce buckling, i.e. deviation from the trivial equilibrium path, a small imperfection needs to be introduced in the computational model. Simulations with two different types of imperfections, depicted in Fig. 6, were performed. The first imperfection consists of an initially debonded area of the width of one element along the yz-symmetry plane of the micromodel. This imperfection was used in most of the simulations. The second imperfection is purely geometrical, in the form of a small initial buckle similar to the buckling mode, at the yz-symmetry plane across the width of one film strip. In this case there is no initial debonding. The initial buckle height is 1/10 of the film thickness. This imperfection was used to check, if the shape of the initial imperfection has an impact on the buckling process.

Fig. 6

Different imperfections to trigger buckling: (a) debonded area, (b) initial buckle.

Results and discussion

Determination of the fracture stress of the film

The fracture stress of the brittle chromium film is determined to be 1200 MPa. This value agrees well with in situ X-ray diffraction measurements of a strained 100 nm thick chromium film on polyimide, which was found to be 1500(±250) MPa). Considering the small grain size of the film (20 nm columnar grains [22]), this high value is realistic. Fig. 7 shows the von Mises equivalent stress in the film at the experimentally determined cracking strain. The stress distribution in the film is almost constant over a wide area in the centre of the specimen.

Fig. 7

Von Mises equivalent stress in the film at a global strain level of 0.4%. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Deformation of the specimen with cracked film

At the unclamped edges the specimen curves downwards, away from the film. The z-displacements of the nodes on the symmetry planes are plotted in Fig. 8. Here it can be seen that the centre of the specimen remains almost flat during the course of the simulation. Since this area was investigated during the experiments, it was chosen for the micromodel. The curvature in the loading direction is generally decreased due to non-linear geometric effects with increasing straining (see Fig. 8). This curvature leads to a higher local strain level in the micromodel compared to the global strain level. The same holds for the strain level transverse to the primary loading direction.

Film strip curvature

As the substrate is strained underneath the cracked film, shear stresses are activated and concentrated at the interface area close to the channel cracks, as described by shear lag models (see e.g. [10]). The shear stresses in the x-direction, i.e. perpendicular to the crack, at the interface also lead to a bending moment in the film. This bending moment is observed as an upward bend of the film at the channel crack edges, as is illustrated in Fig. 9a. Similar bending is found from the simulation (micromodel), as shown in Fig. 9b. The upward bend at the crack edges has been experimentally observed in films held in tension and imaged with an atomic force microscope (AFM) [25,26]. Using a miniaturized, screw-driven tensile frame, it is possible to strain a film–substrate system under an AFM tip to observe the bending of the film between the channel cracks. In Fig. 9c, this is demonstrated for a 200 nm chromium film on a 50 μm thick polyimide substrate strained 15%. It is worth mentioning that this effect cannot be captured by 2-D models.

Fig. 9

(a) Schematic sketch of the upward bending of the edges of a film strip due to shear stresses; (b) cross section of a film strip, far away from the initial imperfection, for a model with and at a global strain level of 4.4%. The simulated model is mirrored with respect to the symmetry planes. Fringe colours show nodal displacements in millimetres; (c) AFM image of a 200 nm Cr film on polyimide under tension strained to 15% and two measured height profiles. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Buckling and delamination

A typical buckle formation is visualized in Fig. 10. The computations have shown that the buckling phenomenon is induced by shear failure of the interfaces along the channel crack edges. The concentration of the shear stresses close to the introduced imperfections triggers localization of the debonding failure at their position. The compressive stresses in the y-direction in the film cause the formation of a buckle at the location of the failed interface. Therefore, a triangular buckle is formed. The two triangular buckles, symmetrically arranged in the strip, grow and propagate to the centre of the strip, where they join and proceed to form a rectangular buckle. The observed process of buckle formation is fundamentally three dimensional. The principal features of this buckling and post-buckling process was found to be widely independent of the used interface parameters and of the type of initial imperfection. Initial mode I failure and the formation of a rectangular buckle, as considered by existing 2-D models [7], has only been achieved by using an initial imperfection similar to the rectangular buckling mode and a very weak interface with values of 20 MPa for the ultimate generalized traction, , and of 0.5 Nm/m2 and 2.5 Nm/m2 for the energy release rates in mode I and mode II, and , respectively.

Fig. 10

Evolution of buckle formation for a model with and : localization around the initial imperfection (a); formation of a triangular buckle (b); creation of a rectangular buckle (c). The simulated model is mirrored with respect to the two symmetry planes. The fringes show the nodal z-displacement in millimetres. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The obtained buckle shapes compare well with the experimentally observed ones for half buckle width (b = 7h) and buckle height (d = 1.4h) [7]. The evolution of the buckle shape is depicted in Fig. 11. The buckle profile remains geometrically similar until the triangular buckle reaches the centre of the film strip. The now rectangular buckle grows wider by rapid delamination (compare Fig. 11 c and d), which is an indication of an unstable portion in the post-buckling path. This supposition is assured by the experimental observation of cracks along the apex of the buckles (see Fig. 1b). The global strain level at which the triangular buckles join at the centre is sensitive to the initial imperfection employed. Therefore, care has to be taken when determining interface properties by comparison of model and experiment.

Fig. 11

Evolution of post-buckling deformation, i.e. of the buckle shape for , and (view from the channel crack perpendicular towards the middle of the strip). The fringes show the nodal z-displacement in millimetres. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Determination of the interface properties

Fig. 12a shows calculated load–displacement curves of the gap width between the film and the substrate at the centre of the strip (at x = y = 0) in relation to the global strain level for a range of different interface properties and initial imperfections. The global strain level at which the triangular buckle reaches the centre of the strip and rapid delamination occurs is clearly visible in the curves as a sudden increase in the gap width at an almost constant global strain. This global strain level will be termed “critical global strain” or “buckling strain” in the following. Fig. 12b shows the dependence of the critical global strain on the energy release rate in mode I for the two modelled initial imperfections. It can be seen that the critical global strain depends on the interface properties. Conducting a parameter study, one can determine the critical global strain as a function of the three interface properties , and . The imperfection similar to the buckling mode is more critical than the one with the initially debonded area, because it leads to earlier buckling with the same interface properties.

Fig. 12

Gap width between film and substrate in the centre of the strip (a); and global strain levels at which separation occurs at the centre (b) for , and different values of and different initial imperfections.

The global strain level at which buckling occurs in the uniaxial tension test experiment can be determined with good accuracy [22]. To estimate the interface properties, parametric studies have to be performed in order to determine which parameter set fits best, when comparing experimental results with computed ones. As a criterion for this comparison, the coincidence between the experimentally obtained and the computed buckling strains is used. This procedure leads to lower bounds for interface strength parameters, as explained in the following: it can be assumed that at least geometrical (buckle-type) imperfections exist in real specimens that are larger than the extremely small one assumed in the model. Since larger imperfections lead to a greater reduction in the buckling strain, it is obvious that a stronger interface than that one assumed in the model must exist in reality if the experimentally determined buckling strain and the one computed for the model coincide. Thus, the interface strength (in terms of interface parameters), which in the simulation leads to the same buckling strain as found in the experiment, represents a conservative estimate. For the specimen under consideration, used as an example, the interface properties of and match best with the experimentally obtained buckling strain of 10%.

The portion of the width of the film strip that delaminates is governed by the value of . Higher values lead to smaller areas of shear failure along the crack. The shear properties of the interface could also be determined from the experimentally measured average crack spacing, as outlined by Xie and Tong [27]. The determination of would allow a reduction in the number of necessary parameter variations. However, for the investigated material combination, the computed crack width is widely independent of the parameter and matches well with experimental measurements. This can be explained by the high stiffness ratio between the film and substrate materials. The straining of the film due to interfacial shear has a minor effect on the crack width, which is primarily dependent on the plastic deformations in the substrate in the vicinity of the crack. However, by a detailed investigation of the curvature of the film strips, it might be possible to gain further insight into the acting shear effects.

Conclusions

The analyses have shown that the debonding between film and substrate is initiated by failure in mode II. In contrast to the assumption of pure mode I failure in 2-D models, the 3-D models show that a combination of mode I and mode II appears in the separation process. The buckling and delamination process is significantly a three-dimensional one. The formation of triangular buckles from the sides of the film strips occurs over a wide range of interface properties, independently of the initial imperfection. Triangular buckles are preferred over the formation of rectangular buckles. The experimentally observed upwards bending of the edges of the film strips are reproduced by the simulation model. The frequently used 2-D models are not able to capture these important effects and are therefore not sufficient to accurately model the buckling and debonding process.

Values of 4.0 Nm/m2 and 8.0 Nm/m2 for the energy release rates in mode I and mode II, respectively, and of 80.0 MPa for the ultimate traction as obtained for the specimen used as an example in the paper, are in reasonable agreement with experimental results and can be considered as estimates of the interface properties.

The presented modelling strategy is applicable to similar material combinations and could be largely automatized to allow an efficient execution of parameter studies to determine the interface properties of thin films on compliant substrates. The modelling concept opens a new way to a more realistic description of the whole experimental procedure of the uniaxial tension test and thus a more accurate determination of the interface properties.

Although the developed model can be used for estimating interface properties, it should be mentioned that the principal aims of this paper were to study the fundamental features of the strip tensile test and to explain the experimental observations by comprehensive theoretical considerations.