Assessment of cohesive traction-separation relationship according stiffness variation

 

Evaluación de la relación tracción-separación cohesiva según la variación de rigidez

 

 

 

Liliana M. Bustamante-Góez1a, Edwin Chica-Arrieta2, Junes A. Villarraga-Ossa 1b

 

 

1 Grupo de Diseño Mecánico, Universidad de Antioquia, Medellín, Colombia. Orcid: a  0000-0002-1150-1252, b 0000-0002-7893-5362.

Email: a liliana.bustamante@udea.edu.co , b junes.villarraga@udea.edu.co

2Grupo de Energía Alternativa, Universidad de Antioquia, Medellín, Colombia. Orcid: 0000-0002-5043-6414.

Email: edwin.chica@udea.edu.co

Abstract

 

The definition of a traction-separation relationship is a fundamental issue in cohesive zone models because it describes the nonlinear fracture process. Cohesive interactions are generally a function of displacement jump (or separation). If the displacement jump is greater than a characteristic length (n), complete failure occurs. In this study, the softening condition behavior of a cohesive interface between two identical materials was assessed for different stiffness values of solid and cohesive. The cohesive interface was modeled with a traction-separation linear relationship and for the solids continuum elastic constitutive models were used. The softening condition was obtained by analytical and finite element method. The whole system behavior was modeled using ABAQUS 6.14 to obtain stress-displacement relationship. The analytical solution and computational results were compared. The computational results matched the analytical solutions and the simulations allowed to obtain a response in the cases where the analytical solution has singularities “backslash effect”. 

 

Keywords:

Abaqus; cohesive element; cohesive zone model; finite element simulation; traction separation law. 

 

Resumen

 

La definición de la relación tracción-separación es una cuestión fundamental en los modelos de zona cohesiva porque describe procesos de fractura no lineal. Las interacciones cohesivas son generalmente una función del desplazamiento (o separación). Si el cambio en el desplazamiento es mayor que una longitud característica (n), ocurre una falla completa. En este estudio la condición de ablandamiento de la interfaz cohesiva entre dos materiales idénticos fue evaluada para diferentes valores de rigidez del sólido y del cohesivo. La interfaz cohesiva fue modelada con la relación lineal de tracción-separación y para los sólidos se utilizaron modelos constitutivos continuos elásticos. El comportamiento de todo el sistema fue modelado usando ABAQUS 6.14 para obtener la relación esfuerzodeformación. La solución analítica y los resultados computacionales fueron comparados. Los resultados computacionales concordaron con la solución analítica y las simulaciones permitieron obtener una respuesta en los casos donde la solución analítica tiene singularidades "backslash effect". 

 

Palabras clave: modelo de zona cohesiva; ley tracción separación; elemento cohesivo; simulación por elemento finito; Abaqus.

 

Introducción

 

Cohesive zone models have been used to treat fracture nonlinear problems since it provides a more realistic feature of the failure process [1]. The cohesive zone is a surface in a bulk material where displacement discontinuities occur. Thus, continuum is enhanced with discontinuities in the form of displacement jumps.  These displacement jumps require a constitutive description called traction -separations relationship (cohesive laws) to describe cohesive interactions. In general, tractionseparation relationships can be classified into potentialbased models and non-potential-based models. Potentialbased models use the concept of cohesive energy potential, for example, Needleman and Tvergaard [2], [3]. For non-potential-based models, several constitutive relationships of the cohesive zone model with various shapes have been developed, e.g., linear softening [4], trapezoidal shape [5], bilinear softening [6], cubic polynomial, and exponential, as shown in Fig.1. All these models, irrespective of the choice of the elementary functions, are constructed qualitatively as follows: tractions increase, reach a maximum, and then approach zero with increasing separation. This scenario is in harmony with our intuitive understanding of the rupture process. It is analogous to atomic interactions [7]–[9]. 

 

Needleman introduced the cohesive zone models (CZMs) in computational practice. Since then CZMs are used increasingly in finite element simulations of crack tip plasticity and creep; crazing in polymers; adhesively bonded joints; interface cracks in bimaterials; delamination in composites and multilayers; fast crack propagation in polymers, and so on [7], [9].

 

One of essential aspects in the CMZ is the choice of a traction–separation relation also called tractionseparation law. Because most of these relationships exhibit limitations, especially under mixed-mode conditions, the relationship should be selected with great caution. Many researchers using CZMs consider that the separation work and the cohesive strength are two main parameters characterizing of the separation process. But, if the cohesive element stiffness is less than the stiffness of the surrounding elements, the global response can be affected when the failure process is computed by computational techniques and the back-slash effect is produced during the softening condition. On the other hand, when the analytical solution is considered to solve this situation, it is no possible obtain a response for cases in which the cohesive element stiffness is equal to surrounding elements stiffness. For this reason, we investigated the softening condition behavior of a cohesive interface between two identical materials for different stiffness values of solid and cohesive elements by analytical and computational methods.

 

In this research work, interfaces along two elastic similar solids were studied by using linear form of CZMs and modeled by element finite approach includes in ABAQUS (section 2). The analytical solution was performed considering two blocks bridged by a cohesive zone of zero thickness due to this assembly allows to verify easily the solution (section 3).  In section 4 we concluded, based on the obtained results that simulations allow to obtain a response in the cases where analytical solutions have singularities called “back slash effect”. 

 

The aim of this study is to establish how the numerical modeling using the finite element method can be used to represent the different cases in which the stiffness modulus of the bulk material can vary and obtain adequate representations of the phenomena. These results will be useful in future work modeling composite materials and obtain their mechanical response by using computer simulations.

 

Cohesive zone model (CMZ)

 

The idea for the cohesive model is based on the consideration that infinite stresses at the crack tip are not realistic, the first models to overcome this drawback were porposed by Dugdale [10] and Barenblatt [11]. For practical applications the model became more interesting when numerical methods, mostly the finite element method, were applicable to nonlinear problems, in 1990 Needleman [12] used the model of crack propagation to analyze ductile materials, since then it is a common practice to analyze with this model the growth of cracks in this type of materials.

 

Displacement of the tip position of a crack assumes bonds stretching orthogonal to the crack surfaces until they break According to CZM, the fracture process is lumped into the crack line and is characterized by a cohesive law that relates tractions and displacement jumps across cohesive surfaces, the whole body volume remains elastic while the nonlinearity is embedded in the cohesive law which dictates the interfacial conditions along the crack line (Fig. 1). 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 [13].

 

Cohesive zone model adds a zone of vanishing thickness ahead of the crack tip with the intention of describing more realistically the fracture process without the use of stress singularity. The cohesive zone is idealized as two cohesive surfaces, which are held together by a cohesive traction (see Figure 1). The cohesive crack propagation may consist of four stages: elastic, initiation, softening and complete failure [1], [14].

 

Since the CZMs are a phenomenological model, there is not a rule evidence cohesive law shape most suitable according to failure process. Therefore, the cohesive relationship is assumed independent of a specific material and many authors use the traction separation relationship created by themselves.

 

The magnitude of the parameters in CZMs vary widely, ranging from MPa to GPa for traction, J to kJ for energy, and nanometers to micrometers for separation distance [15].

 

The intrinsic traction-separation relationships used by ABAQUS, is briefly explained below.

 

 

1.1.  Traction-Separation relationship in ABAQUS

 

Since the cohesive model is a phenomenological model there is no evidence which form to take for T (δ). So it has to be assumed independent from the material as a model quality [16]. In the literature it can be found several approaches (Figure 3). ABAQUS considers the traction-separation relationship as the variation of fracture toughness as a function of a mode-mixity ratio. This model was originally developed by Camanho et al. [4]. To describe the mixed-mode condition across the interface, an effective separation (Δ̅) is defined in equation (1).

 

Where Δn and Δtare the normal and tangential separations, respectively.

 

The available traction-separation model in ABAQUS assumes initially linear elastic behavior followed by the initiation and evolution of damage (see Figure 2). The elastic behavior is written in terms of an elastic constitutive matrix that relates the nominal stresses to the nominal strains across the interface, according to constitutive relationship describes below:

 

 

Where T is nominal traction stress vector, C is the elasticity matrix and is the strain vector defined by, denoting by T0 the original thickness of the cohesive element.

 

1.1.1.  Damage initiation

 

Damage initiation refers to the beginning of degradation of the response of a material point. When the separations reach the effective critical separation, the state of separation corresponds to the onset of damage and/or crack initiation. ABAQUS has several failure criteria. A maximum nominal stress criterion is employed to determine the onset of damage, which the damage is assumed to initiate when the maximum nominal stress ratio (as defined in the expression below) reaches a value equal to one. This criterion is represented by ecuation (3).

 

 

 

Where 〈 〉 is the Macauley bracket,  which represent the normal and the two represent the peak values of the nominal stress when the shear tractions, respectively. deformation is either purely normal to the interface or  purely in the first or the second shear direction, 

 

 

Figure 3. Effective traction-separation relationships: (a) linear softening, (b) trapezoidal, (c) smoothed trapezoidal, (d) bilinear softening, (e) cubic polynomial, and (f) exponential. Source: authors.

 

When the effective separation is smaller than an effective tractions are proportional to the normal and tangential critical separation the normal and tangential separations, which are given as: 

 

 

 

where Kp is a penalty stiffness.

 

1.1.2.  Damage evolution

 

Damage evolution law describes the rate at which the material stiffness is degraded once the corresponding initiation criterion is reached. ABAQUS uses a scalar damage variable D to represent the overall damage at the contact point. This variable takes 0 and 1 values (if the cohesive element is broken, D = 1 on the contrary, D =

0).

 

When the effective separation is greater than the effective critical separationand smaller than the effective complete failure separation i.e.̅ < Δ̅ < 𝛿𝑓̅ , the state of separation corresponds to the softening condition [7]. Then, the normal and tangential cohesive tractions are defined by equation (5).

 

 

The definition of damage evolution in ABAQUS is specified by two components. The first component involves specifying either the effective complete failure separation (𝛿𝑓̅ ) or energy fracture (GC). The second component to the definition of damage evolution is the specification of the nature of the evolution of the damage variable, D, between initiation of damage and final failure (softening condition). Three types of damage evolution are available in ABAQUS: linear, exponential or tabular. A linear model is employed to describe the softening condition according to the equation (6).

 

 

1.1.3.  Mixed-mode definition

 

The mixed mode of deformation fields in the cohesive zone quantify the relative proportions of normal and shear deformation. ABAQUS uses three measures of mixed mode, two are based on energies and the otherone is based on tractions. Mixed mode definitions based on energies are described by equations (7) and (8).

 

 

Where Gn, Gs, Gt are the work done by the tractions and their conjugate displacements in the normal and shear directions.

 

Analytical solution

 

In this study, a body compromising two symmetric rigid parts bridged by a cohesive zone of zero thickness was modeled as shown Figure 4. In addition, a linear cohesive zone model is used to simulate the interfacial mechanical response. For pure opening (Δ𝑡 = 0) the variation of normal traction with respect to Δ𝑛 to solids and cohesive layer are shown in Figure 5 (a) and (b), respectively, thus, the constitutive response of system is shown in Fig. 6. Increasing forces P or applying a vertical displacement will lead to the nucleation of a crack at center if the stiffness solid is less than the cohesive layer stiffness (case i) or, at the right and left side if the stiffness solid is equal or greater than the cohesive layer stiffness (case ii – iii) see Fig. 6.

 

 

An elastic linear behavior to the solid is assumed and using the traction-separation law mentioned in section 2 the following equations (9) and (10) were obtained in each stage of the damage evolution.

 

 

 

 

Where K and K’ are the penalty stiffness and stiffness degradation to linear traction-separation law of cohesive layer, E and L are the Young’s modulus and length of the solid and Δ are displacement jumps normal to the cohesive zone.

 

To assess the cohesive traction-separation relationship according stiffness variation, the values for Young’s modulus of solid are arbitrarily selected as 1, 10, 70 GPa. Mode I fracture is selected. In addition, the penalty stiffness is 50 MPa/mm, effective critical separation is 10 µm, and the cohesive strength is 10 MPa.

Finite element approach

 

To compare the analytical solution with numeral response the finite element method is employed by simulations conducted in commercial software ABAQUS 6.14.1. This software allows solving problems in cases where analytical methods present singularities, for example case ii shown in Figure.

 

The simulated model is shown in Figure 4, on the upper edge of the solid a positive displacement Δ𝑢 of the axis y is applied until cohesive elements fail. The lower border is simulated as fixed. A standard analysis is used for case iii whereas to case i and ii, explicit analysis is employed.  The discretized mesh has 512 linear quadrilateral elements of type CPS4R to model solids. 160 linear quadrilateral elements of type COH2D4 with zero thickness in the direction normal to the interface are used to model interface behavior. The mesh and boundary conditions used in the simulation are shown in Figure 6.

 

Cohesive law is defined in ABAQUS specifying the damage evolution as: effective displacement 𝛿𝑓 − 𝛿𝑐 = 20 𝜇𝑚, a variable D linear softening type, modeindependent behaviour and mixed-mode energy definition. An elastic constitutive behavior is assumed for the bulk.

 

Results and discussion

 

To evaluate the influence of stiffness on the tractionseparation relationship, the model is analyzed with three different stiffness values for the solid. The analytical response was obtained using the equations (9) y (10). The evolution of von Mises stress as the displacement imposed on the upper edge occurs is shown in Figures 8 - 10.  A comparison between the analytical and computational solution is shown in Figures. 11 – 13.

 

In the simulation shown in Figure 7 the system presents instabilities due to the bulk is compliant as compared to the rigidity of the cohesive layer, which makes it necessary to use explicit analysis to control the deformation rate and convergence of system. Figure 8 shows a behavior similar to that of Figure 9, but the analytical response is significantly affected since for larger values of the critical separation the system has a singularity called the back slash effect. This is shown in stress-displacements plots computational through waves.

Since at this point the system becomes unstable.  

Observing  the  results  shown  in  the  figures  11-13  it  is evident  that  the  static  and  explicit  analysisreportthe same  tendency  and  the  differences  between  the  graphs are due to the inertial effects of the system for the cases in  which  the bulk  is  more  compliant  than  the  interface. Figure 12shows well agreement between both solutions. In  this  case  the  bulk  is  stiffer  than  the  cohesive  layer hence,    the    global    response    is    stable    and,    the computationalsolution   is   not   depending   of   time increment used. Therefore, standard analysis can be used in this situation

Summarizing, it is recommended to use lower cohesive interface stiffness values than the matrix to avoid convergence problems and stability of simulations. In Figure 12. Stress - displacement plots analytical vs  cases where the interface is more compliant than the computational case iii with E = 70 GPa, max = 10 MPa, matrix or bulk, an explicit type analysis is required using very small time increments in order not to affect the K = 50 MPa/mm, c = 10 µm and f  = 40 µm. Source: authors. global system response and thus the peak strength and the

 fracture energy remain unchanged

 

Conclusions

 

Using the finite element method for the growth of cracks by CMZ is a reliable and economical tool to predict the behavior of these systems.

 

Since the CZM is a phenomenological model is important to select an adequate traction-separation law and carry out a robust validation process to guarantee the quality of the results. It is recommended to use explicit simulations for this type of systems, since these behave better than the implicit simulations.

 

Lower cohesive interface stiffness values than the matrix help to avoid convergence problems and give stability to the simulations.

 

Finally, it is observed that for case iii, the simulations predict very well the behavior of the system, while for the other two cases the dynamic effects have a great impact on the stability, for which it is recommended to use stiffness values corresponding to the case iii, to obtain results with low hysteresis.

 

Acknowledgements

 

The authors wish to acknowledge the Program Education for Research Colciencias for providing financial support to doctoral studies in Colombia.

 

References

 

[1]                C. T. Sun and Z.-H. Jin, “Chapter 9 - Cohesive Zone Model,” in Fracture Mechanics, C. T. Sun and Z.-H. Jin, Eds. Boston: Academic Press, 2012, pp. 227–246.

 

[2]                V. Tvergaard and J. W. Hutchinson, “The relation between crack growth resistance and fracture process parameters in elastic-plastic solids,” J. Mech. Phys. Solids, vol. 40, no. 6, pp. 1377–1397, 1992.

 

[3]                X.-P. Xu and A. Needleman, “Numerical simulations of fast crack growth in brittle solids,” J. Mech. Phys. Solids, vol. 42, no. 9, pp. 1397–1434, 1994.

 

[4]                P. P. Camanho, C. G. Davila, and M. F. de Moura, “Numerical Simulation of Mixed-Mode Progressive Delamination in Composite Materials,” J. Compos. Mater., vol. 37, no. 16, pp. 1415–1438, 2003.

 

[5]                V. Tvergaard and J. W. Hutchinson, “The influence of plasticity on mixed mode interface toughness,” J. Mech. Phys. Solids, vol. 41, no. 6, pp. 1119–1135, 1993. [6] P. H. Geubelle and J. S. Baylor, “Impact-induced delamination of composites: a 2D simulation,” Compos. Part B Eng., vol. 29, no. 5, pp. 589–602, 1998.

 

[7]  K. Park, H. Choi, and G. H. Paulino, “Assessment of cohesive traction-separation relationships in ABAQUS: A comparative study,” Mech. Res. Commun., vol. 78, pp. 71–78, 2016.

 

[8]  G. H. Paulino, “Cohesive Zone Models: A Critical Review of Traction-Separation Relationships Across Fracture Surfaces,” Appl. Mech. Rev., vol. 64, no. 6, 2011. doi: 10.1115/1.4023110.

 

[9]  V. Tvergaard and J. W. Hutchinson, “The relation between crack growth resistance and fracture process parameters in elastic-plastic solids,” J. Mech. Phys. Solids, vol. 40, no. 6, pp. 1377–1397, 1992. doi: 10.1002/cnm.717.

 

 

[10]                       D. S. Dugdale, “Yielding of steel sheets containing slits,” J. Mech. Phys. Solids, vol. 8, no. 2, pp. 100–104, 1960. doi: 10.1016/0022-5096(60)90013-2.

 

[11]                       G. I. Barenblatt, “The Mathematical Theory of Equilibrium Cracks in Brittle Fracture,”, , vol. 7, H. L. Dryden, T. von Kármán, G. Kuerti, F. H. van den Dungen, and L. Howarth, Eds. Elsevier, 1962, pp. 55– 129. doi: 10.1016/S0065-2156(08)70121-2.

 

[12]                       A. Needleman, “An analysis of decohesion along an imperfect interface,” Int. J. Fract., vol. 42, no. 1, pp. 21– 40, Jan. 1990. doi: 10.1007/BF00018611.

 

[13]                       M. Alfano, F. Furgiuele, A. Leonardi, C. Maletta, and G. H. Paulino, “Cohesive Zone Modeling of Mode I Fracture in Adhesive Bonded Joints,” Key Eng. Mater., vol.    348, pp.    13–16,       2007.        doi: 10.4028/www.scientific.net/KEM.348-349.13.

 

[14]                       K. Park and G. H. Paulino, “Computational implementation of the PPR potential-based cohesive model in ABAQUS: Educational perspective,” Eng. Fract. Mech., vol. 93, pp. 239–262, 2012. doi: 10.1016/j.engfracmech.2012.02.007.

 

[15]                       N. Chandra, H. Li, C. Shet, and H. Ghonem, “Some issues in the application of cohesive zone models for metal–ceramic interfaces,” Int. J. Solids Struct., vol. 39, no. 10, pp. 2827–2855, 2002. doi: 10.1016/S00207683(02)00149-X.

 

[16]                       I. Scheider, “Cohesive model for crack propagation analyses of structures with elastic–plastic material behavior Foundations and implementation,” GKSS research center Geesthach, Teltow, apr. 2001.

 

[17]                       Abaqus 6.14 Analysis User’s Manual, Dassault Sistemes, Simulya, 2014.