FRACTURE ANALYSIS OF ADHESIVE JOINTS USING INTRINSIC COHESIVE ZONE MODELS

M. Alfanoa*, F. Furgiuelea, A. Leonardia, C. Malettaa, G. H. Paulinob Department of Mechanical Engineering, University of Calabria, 87036 Arcavacata di Rende, Italy, e-mail: alfano@unical.it b Department of Civil and Environmental Engineering, University of Illinois at Urbana-Champaign, Newmark Laboratory, 205 N. Mathews Avenue, Urbana, IL 61801-2352, U.S.A.
a

ABSTRACT In this paper Cohesive Zone Model (CZM) concepts are applied in order to study mode I fracture in a pre-cracked bonded Double Cantilever Beam (DCB) specimen. A cohesive surface element is implemented in a Finite Element commercial code using intrinsic cohesive zone models: exponential, bilinear and trapezoidal traction-separation laws. The sensitivity of cohesive zone parameters in predicting the overall mechanical response is examined, then the load displacement curves obtained with the different CZMs are compared and some interesting features concerning the prediction of damage onset in adhesive joints are illustrated. Finally, cohesive parameters are identified comparing numerically simulated load-displacement curves with experimental data retrieved from literature.

1. INTRODUCTION The use of adhesive joints in automotive, aerospace, biomedical, and microelectronics industries is widespread; however, inaccurate joint fabrication or inappropriate curing may cause the presence of bubbles, dust particles or un-bonded areas in the bond line. As a consequence the establishment of reliable integrity assessment methodologies and testing procedures is needed in order to tackle fracture events in adhesive joints. A valuable approach from this point of view is represented by the Linear Elastic Fracture Mechanics (LEFM) [1-3]. So far, both the classic one parameter models of elastic and plastic fracture have been successfully applied but they present some limitations. They require pre-existing crack like defects or notch and, therefore, crack initiation cannot be treated directly. In addition, they neglect a detailed description of what happens in the fracture process zone (FPZ) because, provided the FPZ is small compared to the other specimen dimensions, they lump it all into the crack tip. However, a detailed description of the FPZ is essential, especially to understand fracture mechanisms and to design suitable modifications of the material (e.g. toughening by reinforcement in polymeric structural adhesive [4]). The standard model used to describe the crack tip process zone assumes bonds stretching orthogonal to the crack surfaces until they break at a characteristic stress level. Thus, the singular region introduced from LEFM can be replaced by a lateral region over which non-linear phenomena occur. This model is often mentioned as the Cohesive Zone Model (CZM) and it can be traced back to the works of Dugdale and Barenblatt [5,6]. According to CZM the entire fracture process is lumped into the crack line and is characterized by a cohesive law that relates traction and displacement jumps across cohesive surfaces (T-). Unlike fracture mechanics based strategies, CZM can be used for the analysis of crack initiation and growth that, indeed, are obtained as a natural part of the solution without any a priori or ad hoc assumptions. So far, CZM has been successfully applied to model fracture in metals, concrete, polymers and functionally graded materials (FGMs) [7-11]. Cohesive zone approaches differ in the way by which cohesive surface elements are inserted in the initial geometry. In the intrinsic approach [12] cohesive elements are introduced between volumetric elements from the beginning of the analysis as a network of cohesive surfaces; on the other hand, in

the extrinsic approach, [13] a cohesive element is introduced in the mesh only after the corresponding interface is predicted to have started to fail, i.e. they are inserted adaptively. In this paper CZM concepts have been applied in order to study mode I fracture in a pre-cracked bonded Double Cantilever Beam (DCB) specimen. A cohesive surface element has been implemented in the Finite Element commercial code ABAQUS [14] using intrinsic cohesive zone models: exponential, bilinear and trapezoidal laws. Theoretical and numerical aspects of the CZM are explained. In particular, details regarding the computation of the force vector and the tangent stiffness matrix are presented. The sensitivity of cohesive zone parameters in predicting the overall mechanical response is examined and then cohesive parameters are identified comparing numerically simulated load-displacement curves with experimental data retrieved from literature [2].

2. GENERAL THEORY AND NUMERICAL ASPECTS OF COHESIVE ZONE MODEL (CZM) 2.1 CZM concepts Let consider the domain showed in Fig. 1.a which is divided into two parts, i.e. 1 and 2, by a material discontinuity, c. This last defines the interface between the domains 1 and 2 and represents an internal surface not yet separated. Prescribed tractions, fi, are imposed on the boundary f. The stress field, ij, in , neglecting body forces, is related to the external loading and to the closing tractions, Ti, in the material discontinuity through the following equilibrium equations

ij x j n ij j n ij j ijn j

=0 = fi = ri = Ti

in on on on

= 1 U 2 f u c
(1)

where nj is the outward normal, while ri are the reaction forces on the boundary u. Cohesive zone model is, as previously said, one of the most commonly used tools to model material discontinuities. In the simplest and most usual formulation of CZM, the whole body volume remains elastic while the nonlinearity is embedded in the cohesive law which dictates the boundary conditions along the crack line, c (Fig. 1). The pick stress on the cohesive law is the cohesive strength, cr, of the material while the area under the curve is the cohesive fracture energy, Gc. As a consequence, fracture process can be summarized as illustrated in Fig. 1.b: at first a linear elastic material response prevails (I), as the load increases the crack initiates (T=cr) (II) and then, governed by the non linear cohesive law (softening curve), it evolves from initiation to complete failure (III) with the appearance of new traction free crack surfaces, c and +c (=f) (IV).

Fig. 1 - Cohesive Zone Model concepts

Fig. 2 - Intrinsic versus extrinsic CZM

Therefore, the continuum should be characterized by two constitutive laws: a linear stress-strain relation for the bulk material and a cohesive surface relation (cohesive law) that allows crack spontaneous initiation and growth. It was shown in [15] that the shape of the traction-separation laws plays an important role in the macroscopic mechanical response of the system, therefore, it is

important to properly select the shape of the softening curve. The development of cohesive zone models in FEM framework requires bulk finite elements, for modelling the stage (1) (Fig.1.b), bordered by cohesive surface elements for the remaining three stages. From this point of view there are basically two cohesive zone approaches that differ in the way by which cohesive surface elements are inserted in the initial geometry. In the intrinsic approach [12] cohesive elements are introduced between volumetric elements from the beginning of the analysis as a network of cohesive surfaces. Cohesive tractions increase from zero to a failure point that is represented by the cohesive strength, at which they reach a maximum before they gradually decrease back to zero following the post peak softening behaviour (Fig. 2). This active network of cohesive surfaces introduces a compliance due to the initial slope of the cohesive law, i.e. the penalty stiffness, K0. This drawback could be addressed increasing K0, however, from a numerical standpoint, it is not recommendable to use elements with an extremely stiff initial response because instabilities in the solution procedure arise [16]. Another noteworthy CZM reported in literature is the extrinsic model proposed in [13], which eliminates the artificial compliance typical of the intrinsic models mentioned above. In particular, in the extrinsic approach, a cohesive element is introduced adaptively in the mesh only after the corresponding interface is predicted to have started to fail; as a consequence, it is possible to adopt a cohesive law with a very high initial stiffness like that reported in Fig. 2. In the present paper, the fracture behaviour of adhesive joints is analyzed using intrinsic cohesive zone model. For this particular problem cohesive surface elements are introduced along a predefined fracture path, i.e. the bond line; therefore their number is reduced and this, in turn, yields a lower compliance. 2.3 Traction separation laws examined in the paper In this paper three of the most representative traction-separation laws proposed in literature have been analyzed and compared in order to investigate the fracture behaviour of adhesive joints: bilinear, trapezoidal and exponential cohesive laws (Fig. 3). In particular, only pure-mode I decohesion problems are considered so that the attention has been focused on traction-displacement relationships relating the normal opening displacement to the dual traction component T. In all the traction separation laws presented below, only positive values of are considered, for negative values, a high penalty stiffness is assumed in order to avoid crack faces interpenetration (Fig. 3). In addition, only monotonic loading processes, without local unloading, will be studied in the paper. This assumption strictly makes the laws thermodynamically not consistent (i.e. irreversibility is not accounted for) however it will not affect the results of the comparative analyses. In what follows, the tractionseparation laws analyzed in the paper are briefly summarized. The cohesive traction according to the exponential model is derived trough cohesive potential energy [12]
exp = cr cr 1 1 + exp(1) cr cr

T=

= cr exp1 cr cr

(2a, 2b)

where cr is the critical normal displacement and cr, as mentioned above, is the cohesive strength of the material. As crack face displacements increase, the traction increases, reaches a maximum, and then decays monotonically. The traction integrated to complete separation yields the fracture energy release rate
G c = cr cr exp(1)

(3)

The analytical expression of the bilinear law [17] is given as follows

(4)

where K0 is the penalty stiffness of the cohesive zone model that can be adjusted selecting proper values of cr and 1 (Fig. 3). The normal work of separation is then given by
G c = cr f / 2

(5)

The analytical expression of the trapezoidal traction-separation law is given as follows [18]

K0 cr = 1 f cr < 2 = 2 f cr T= 2 < f cr ( f ) /( f 2 ) > f 0

(6)

the normal work of separation is given by

(7)

As it can be seen, for this particular (T-) curve, the governing cohesive parameters are the cohesive fracture energy (Gc), the peak stress (cr), the critical opening displacement (f) and the factors 1 and 2 whose values dictate the shape of T=T().

Fig. 3 Cohesive laws examined

Fig. 4 - Cohesive zone element

2.3 Finite element implementation The development of cohesive zone models in FEM framework requires bulk finite elements, for modelling the stage (1) (Fig.1.b), bordered by cohesive surface elements for the remaining three stages: (2) crack initiation, (3) crack evolution and (4) complete failure. In this work, four node cohesive zone elements (CZE) with two integration points have been implemented within the commercial FE code ABAQUS [14] using the user element (UEL) capability. A schematic representation of a CZE is given in Fig. 4. A CZE is made up of two linear line elements (cohesive surfaces) that connect the faces of adjacent elements during the fracture process. The two surfaces initially lie together in the unstressed deformation state (zero thickness) and, subsequently, separate as the adjacent elements deform. In particular, the surface constitutive relations prescribe the evolution of tractions (T) generated across crack faces as a function of displacements jump () at the nodes of the element. Thus cohesive elements do not represent any physical material, but describe the cohesive forces when material elements are being pulled apart. The insertion of cohesive surface elements bridges linear elastic and fracture behaviour allowing for spontaneous crack propagation. In order to carry out the iterations of the method [14], the contribution of cohesive elements to the tangent stiffness matrix as well as to the force vector is acquired from the numerical implementation of the CZM. The implementation of a general cohesive element is explained in what follows. The line cohesive element has eight degrees of freedom. In particular, the nodal displacements vector in the global coordinate system is given as:
ug = ux

( 1)

uy

( 1)

ux

(2)

uy

(2)

ux

(3)

uy

(3)

ux

(4)

uy

(4)

(8)

in which the order follows typical ABAQUS conventions. Crack faces opening for cohesive elements is defined as the difference between top and bottom nodes thereby leading to the following definition in terms of displacements of paired nodes:

 x (1,4 ) (1,4 ) y ( 2,3 ) = x y ( 2,3 )

u x (1) (1) u y (2) u x u y ( 2 )

(4) ux 1 (4) 0 uy = Lu g = (3) 0 ux (3) uy 0

0 1 0 0

0 0 1 0

0 0 0 1 0 0 0 0 0 1 u g 0 1 0 0 0 1 0 1 0 0

(9)

with [L] being an operator matrix. From the nodal positions, the crack face opening is interpolated to the Gauss integration points by means of standard shape functions
1 ~ x 2 = = y 0 1+ 2 0 (1, 4 ) x (1, 4 ) 0 y = N Lu g 1 + x ( 2 , 3 ) 2 y ( 2 , 3 )

0 1 2

(10)

where is the natural coordinate and N is the matrix of shape functions. Since the constitutive relations are based on tractions and displacements in the local coordinate system a transformation from global to local coordinate is needed for the cohesive element. Let [R] defines the orthogonal transformation matrix from global (x,y) reference frame to element specific local coordinate system. Then the relative displacement vector, for an uncoupled cohesive law, is given as:
= RNLug = Bug

(11)

The relative displacements of the element faces create normal and shear displacements, which in turn generate element stresses depending on the constitutive equations of the material. The relationship between tractions and displacements is given by:
Tt T t = 0 0 Tn n

T T=

(12a, 12b)

with the last term being the Jacobian stiffness matrix while the subscripts t and n denote the tangential and normal directions. The constitutive relationships adopted herein have been presented in the previous section and they are independent of the element formulation. The stiffness matrix for cohesive elements can be obtained by minimizing the total amount of potential energy:
= U+ W = 1 T T T dA u g f 2

(13)

where f is the external traction vector; after some manipulations it follows

(14)

so that the cohesive element stiffness matrix is given as

(15)

where w is the element width. The contribution of cohesive elements to the global force vector is defined in a variational setting using the principle of virtual work
INT = EXT T T dA = ug FCOH
T

(16)

with INT and EXT being the internal and external virtual work respectively. After some manipulations the equivalent right hand side nodal force vector for cohesive elements, FCOH, is given as follows:
FCOH = ABT TdA = w B T Tdc
c

(17)

3. SIMULATION OF FRACTURE IN DCB CONFIGURATION The specimen analyzed herein is the DCB studied in [2]. A schematic representation is reported in Fig. 5. The beam has a length (L) of 120 mm, a substrate thickness (h) of 15 mm, and a width (B) of 30 mm. The mechanical pre-notch (a0) extends 40 mm from the left to the right edge of the beam. The substrate material is aluminium (Es=70 GPa and s=0.33) and the adhesive is epoxy (Ea=1.3 GPa and a=0.35). The bond line thickness (ha) is equal to 0.3 mm. External loading is imposed under displacement control.

Fig. 5 - Schematic representation of the specimen The measured fracture toughness reported in [2] is Gc= 550 N/m. The cohesive strength, that is a material parameter, is not easily measurable, therefore estimated values are usually adopted [19]. Likewise, in this paper the cohesive strength will be selected using a trial and error procedure until a match between numerical and experimental P- curves will be obtained. The general purpose commercial code ABAQUS has been adopted in order to model the specimen: zero thickness cohesive surface elements were used for the cohesive zone and plane strain four nodes isoparametric elements (CPS4) for the surrounding bulk material. Fig. 6 shows mesh details for the regions where cohesive elements are inserted.

Fig. 6 - Details of the FE model in the region were cohesive elements are inserted

Fig. 7 - Total dissipated fracture energy

In order to properly represent tractions in the cohesive zone, spatial discretization along the crack propagation path is essential and enough elements (three or more) should be inserted in the FPZ. With coarser mesh, the shape of the interface law may have a non negligible influence on the simulated load-displacement curve; however, this problem can be tackled by means of suitable mesh refinements [20]. In the paper, different cohesive element sizes, have been assessed. To illustrate that the element size chosen for the numerical analyses is objective and not sensitive to artefacts of the numerical solution, a global quantity, i.e. total dissipated fracture energy, is evaluated for each element size and then compared. As it can be seen from Fig. 7, with element size ranging from 0.1 to 1 mm, the fracture energy dissipated is almost the same. Therefore, using an element length of 1 mm

the proper mechanics of energy dissipation and crack propagation are ensured. Subsequently, a sensitivity analysis to cohesive fracture parameters has been performed. The influence of the simulated mechanical responses to cohesive parameters, i.e. material strength (cr) and cohesive fracture energy (Gc) has been assessed. Fig. 8 illustrates the sensitivity of P versus the opening displacement curve to different fracture energies and cohesive strengths for exponential and bilinear cohesive laws.

(a)

(b)

(c)

(d)

Fig. 8 - Sensitivity to cohesive strength and fracture energy. Fig.8 clearly shows that when the cohesive strength, cr, is held constant (Fig. 8.a-b) as the fracture energy increases the area under the curve (global fracture energy) and the peak load increases. On the other hand, when the cohesive fracture energy is held fixed (Fig. 8.c-d) as the critical strength increases, the peak load is increased while the global fracture energy is almost constant. The authors obtained similar results for the trapezoidal law [21], but, for the sake of brevity, these results are not reported herein. The load-displacement curves predicted by the cohesive laws using similar fracture parameters are now compared (Fig. 9). As it can be seen, a good agreement between the cohesive laws is observed in the linear elastic range and in the post-peak regions (damage propagation) of the P- curves. However, slight differences in the maximum load region (damage onset) are observed, in particular the trapezoidal law overestimates the peak load. This can be addressed as follows; the trapezoidal law produces a different interfacial stress profile with respect to the other ones. These in turn can greatly affect the simulated P- curve if the FPZ is large, i.e. for very stiff bulk materials [20]. The parameters of the cohesive zone model are then tuned in order to match the experimental results reported in [2]. These last are calibrated by fitting present numerical results to experimental findings in order to account for the differences between experiments and numerical simulations. The simulated load displacement curves along with the corresponding cohesive parameters that give minimum deviations between simulations and experiments are reported in Fig. 10. As the tangent stiffness for the exponential model starts to decrease from the very beginning, the initial slope of the corresponding load-displacement curve is affected by this inherent artificial compliance. When bilinear and trapezoidal models are used this compliance is reduced. It is worth noting that the cohesive strength adopted for the exponential model is higher than that of the bilinear one; this is a drawback of the exponential model that does not allow to control the penalty stiffness, K0, without affecting the cohesive strength (K0 increases if cr increases). On the contrary, in the linear CZMs, K0 can be

adjusted ensuring a stiff connection between cohesive and bulk elements and thus a proper representation of the undamaged state, without affecting the cohesive strength. However, the lowest value of the cohesive strength among those that allow to obtain the best fit with experimental data is that of the trapezoidal law; the differences with respect exponential and bilinear laws are -25 % and -14 %, respectively. Indeed, the trapezoidal law overestimates the load for damage onset and then lower values of the cohesive strength must be selected in order to match the experimental results.

Fig. 9 - Comparison among the P- curves keeping cohesive parameters as constants

Fig. 10 - Experimental versus numerical P- curve

4. SUMMARY AND CONCLUSIONS In this paper CZM concepts have been applied in order to study mode I fracture in a pre-cracked bonded Double Cantilever Beam (DCB) specimen. A cohesive surface element has been implemented in a finite element commercial code using intrinsic cohesive zone models: exponential, bilinear and trapezoidal laws. Basically, a good agreement among numerical and experimental results has been observed and some interesting features were reported. In particular, it has been shown that both the linear cohesive zone models are appropriate for modelling the undamaged state of cracked adhesive joints: they can control the pre-peak slope of the traction-separation law allowing to reduce compliance. However, the trapezoidal law overestimates the load for damage onset and then lower values of the cohesive strength must be selected in order to match the experimental results. This problem was addressed to the different interfacial stress profile produced by the trapezoidal law which, in turn, can greatly affect the simulated P- curve for larger FPZ.

REFERENCES [1] Kinloch A.J.: Adhesion and Adhesives Science and Technology, London: Chap. & Hall, 1986. [2] Pirondi A., Nicoletto G. Proc. of the Italian Group of Fracture, Bari, 2000. (in Italian) [3] Alfano M., Furgiuele F., Maletta C. Key Eng Mat, 2006; 325:149-152. [4] Williams J.G.: Fracture Mechanics of Polymers, NY: Halsted Press, John Wiley & Sons, 1984. [5] Dugdale D.S. J Mech Phys Solids, 1960; 8:100-104. [6] Barenblatt G.I. Adv Appl Mech, 1962; 7:55-129. [7] Li W., Siegmund T. Eng Fract Mech, 2002; 69:2073-2093. [8] Roesler J.R., Paulino G.H., Park K., Gaedicke C. Cem & Concr Comp, 2007; 29:300-312. [9] Song S.H., Paulino G.H., Buttlar W.G. ASCE J Eng Mech, 2006; 132:1215-1223. [10] Song S.H., Paulino G.H., Buttlar W.G. Eng Fract Mech, 2006; 73:2829-2848. [11] Zhang Z., Paulino G.H. Int J Plasticity, 2006; 21:1195-1254. [12] Needleman A. ASME J Appl Mech, 1987; 54:525-531. [13] Camacho G. T., Ortiz M. Int J of Solids and Structures, 1996; 33:2899-2938. [14] ABAQUS Users Manual, v6.5-1. Pawtucket, USA, Hibbit, HKS Inc; 2002. [15] Chandra N., Li H., Shet C., Ghonem H. Int J of Solids and Structures, 2002; 39:2827-2855. [16] Alfano M., Technical Report, University of Illinois at Urbana-Champaign, 2006. [17] Geubelle P.H., Baylor J.S. Compos Part B, 1998; 29:589-602. [18] Tvergaard V., Hutchinson J. W. J Mech Phys Solids, 1992: 40:1377-1397. [19] Alfano G., Crisfield M.A. Int J Numer Meth Eng, 2001: 50:1701-1736. [20] Alfano G. Compos Science and Tech, 2006; 66:723-730. [21] Alfano M., Furgiuele F., Leonardi A., Maletta C., Paulino G. H. Key Eng Mat, 2007; 348:13-16.