On the 3D extension of failure models for adhesive joints under mixed-mode fracture conditions: LEBIM and CZM

On the 3D extension of failure models for adhesive joints under mixed-mode fracture conditions: LEBIM and CZM L. T´avara, J. Reinoso∗, A. Bl´azquez, V. Mantiˇc Grupo de Elasticidad y Resistencia de Materiales, Escuela T´ecnica Superior de Ingenier´ıa, Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092, Sevilla, Spain Abstract Interface models for triggering damage events have been extensively used in many engineering applications due to their inherent versatility and relatively simple numerical implementations. Usually, the non-linear behaviour of interfaces and thin adhesive joints between solids has been modelled using Cohesive Zone Models (CZMs) although other alternative models have been proposed so far. This is the case of the recently proposed Linear Elastic-Brittle Interface Model (LEBIM). The LEBIM is characterized by a linear constitutive law (traction-relative displacement relationship) in the undamaged state upon failure. Damage occurs when a material point reaches a critical energy release rate leading to an abrupt failure representation. Currently, from the numerical point of view, a great deal of research dealing with interface debonds and adhesive failure is still mainly focused on 2D models. The aim of the present investigation is to extend the LEBIM (originally proposed for 2D problems) for its applications to 3D problems under general fracture conditions and to examine its reliability against experimental data and cohesive-based interfaces pinpointing the main theoretical and numerical differences with respect to alternative methods. Thus, we focus our attention in addressing the main attributes of the proposed LEBIM over nonlinear CZMs: (i) a simple interface law without the need of triggering a gradual stiffness degradation, (ii) a reliable accuracy in comparison with experimental results, endowing a clear identification of the crack front due to the preclusion of the fracture process zone concept, and (iii) a very high efficiency in terms of mesh requirements stemming from the absence of nonlinear cohesive crack representation. This applicability is shown via the implementation of the 3D LEBIM into the general purpose Finite Element Method (FEM) package ABAQUS through the user subroutine UMAT. Finally, special attention is given to the adequate treatment of the mixed mode fracture that may appear in complex 3D crack propagation cases. Keywords: A. Debonding; B. Linear Elastic Interface Model; C. FE-modeling; D. Cohesive Zone Model 1. Introduction Interface formulations constitute plausible numerical tools in order to simulate the development of failure events at adhesive joints and interfaces between solids. These modelling tools are of especial interest in heterogeneous and layered structures, where the existence of significant differences between the mechanical properties of adjacent regions/components can lead to the onset and growth of failure at their interfaces. This is the case of delamination events in layered carbon and glass fiber reinforced polymer (CFRP and GFRP, respectively) structures [1, 2, 3, 4, 5, 6], thin layer-substrate systems [7], dynamic and fatigue loadings in engineering structures and materials [8, 9, 10], among many others. The development of this kind of interface models has attracted an increasing attention in the last decades due to their high versatility for different loadings and geometrical configurations along with their relative ∗Corresponding author Email address: [email protected] (J. Reinoso) Preprint submitted to ——— January 3, 2024 This is an Accepted Manuscript of an article published by Elsevier in Theoretical and Applied Fracture Mechanics, Volume 100, on April 2019, available at: https://doi.org/10.1016/j.tafmec.2019.01.004. Copyright 2019 Elsevier. En idUS Licencia Creative Commons CC BY-NC-ND simplicity for the corresponding numerical implementations into finite element (FE) or boundary element (BE) codes [3, 5, 11, 12, 13]. The basic idea of most of the existing interface formulations exploits the concept of the cohesive zone originally proposed by Barenblatt [14] and Dugdale [15], who introduced the consideration of a cohesive region (process zone) ahead of the actual crack tip. This cohesive zone (CZ) is characterized by the existence of cohesive tractions (σ), which are in general nonlinear functions of the relative displacements (δ) between the two faces of the crack. The relationship σ−δis usually denominated, in the related literature, as traction separation law (TSL), and can be expressed in the following general form: σ=σ(δ, d) (1) where dis an internal variable (history-dependent) which tracks the degradation till complete failure. The area under the TSL curve identifies the fracture energy Gc. Restricting the discussion to pure mode I fracture conditions, with respect to the profile of the TSL, there is a wide range of cohesive laws such as the basic constant law proposed by Dugdale [15], bilinear laws [3, 5, 11, 16], multilinear laws [16, 17, 18], polynomial-based idealizations [19, 20], exponential formulations [1, 8, 12], among many others. Regarding the relationship between the different fracture modes, recently, Dimitri et al. [21] differentiated between two categories of mixed-mode CZMs: (i) uncoupled formulations, and (ii) coupled models. Particularly, uncoupled formulations treat the fracture modes independently from each other, and therefore they could preclude the use of a mixed-mode fracture criterion for triggering the interface or the adhesive failure. Conversely, in coupled formulations, the traction vector components depend on all the components of the separation between both flanks of the crack. From the computational standpoint, most of the existing cohesive-like models adopt two basic definitions in terms of their insertion into the corresponding discretization scheme: (i) the intrinsic cohesive technique through which the interface elements are defined in the pre-process stage, i.e. prior to running the simulations, and is characterized by the consideration of an initial interface stiffness upon damage initiation [6, 22], or (ii) the extrinsic cohesive technique, which exclusively accounts for the softening branch of the TSL. A possible alternative to the previous CZ formulations is the so-called Linear Elastic-Brittle Interface Model (LEBIM), which was originally proposed in [23] and recently extended, e.g., in [24, 25, 26]. This constitutive model is especially suitable for weak or imperfect interfaces and for isotropic thin adhesives according to [27, 28, 29], which are identified as relatively soft thin layers that are confined between two stiffer adjacent bodies. These fundamental modelling hypotheses lead to assume a constant stress distribution over the thickness direction of the interface. From a numerical perspective, the LEBIM does offer several attractive aspects in comparison to standard CZMs such as its robustness, simplicity and computational efficiency, among others. This approach has been examined in different micro-mechanical and macro-mechanical composite applications, leading to very accurate results in comparison with the respective experimental data and computational efficiency with respect to standard intrinsic CZMs using both Boundary Element (BE) and Finite Element (FE) discretization strategies [4, 30, 31]. Inspired by these previous investigations, the current study outlines the 3D extension of the LEBIM to simulate failure of thin adhesive joints with special attention to mixed-mode fracture conditions. In this concern, we aim at pinpointing the principal attributes of the LEBIM model, which correspond to: (i) an interface law that is characterized by an uniform distribution of the linear springs upon failure leading to a simplified numerical implementation, (ii) a reliable accuracy against experimental results with a clear identification of the crack front derived from the lack of fracture process zone concept in the representation of damage process, and (iii) a high computational efficiency in terms of mesh requirements. Moreover, in order to clarify the principal modeling assumptions of the developed formulation, by the first time, a comprehensive comparison with respect to a standard CZM is addressed, whose formulation is omitted here for conciseness reasons, see [32] and the references therein given for a more comprehensive presentation of the CZM. This discussion aims at highlighting the principal differentiating aspects between both models from the conceptual and numerical standpoints. The manuscript is organized as follows. Section 2 introduces the 3D version of the proposed LEBIM, emphasizing its outstanding aspects and the principal differences with respect to a standard CZM obeying 2 a bilinear TSL. The applicability and accuracy of the current formulation to trigger interlaminar damage is examined in Section 3 by several applications, where special attention is devoted to the analysis of 3D fracture events involving mixed mode fracture conditions. Finally the main conclusions of this investigation are summarized in Section 4. 2. Linear Elastic-Brittle Interface Model (LEBIM): 3D constitutive model and finite element formulation In this section, the main aspects of the proposed 3D LEBIM using a finite thickness concept in order to model failure in adhesive joints are introduced. Section 2.1 outlines the constitutive interface law along with the chosen failure criterion to account for decohesion. In order to clarify the main aspects of the proposed model, a direct comparison with a bilinear cohesive interface formulation is presented in Section 2.2. Finally, the numerical treatment of the proposed LEBIM within the context of nonlinear implicit FEM is addressed in Section 2.3. 2.1. Constitutive interface law The 3D interface model under consideration is a direct extension of the 2D constitutive interface law introduced by the authors in [25, 26, 33], so that only the fundamental aspects of the formulation are presented in the sequel. Let consider two adjacent bulk regions B(i)(i= 1,2) with the presence of a weak interface ¯ Bof initial thickness (h > 0), whose points are denoted by the position vector X∈¯ B, see Fig. 1. It is assumed that this initial thickness is significantly smaller than the characteristic thickness of the adjacent bulks, and consequently a constant stress state over the thickness direction can be recalled. This assumption holds for a particular range with regard to the interface thickness hand the stiffness of the adjoining bodies, see [27, 28]. According to [2], from the operative point of view, this initial thickness may be approximately set about 1/100 of the characteristic thickness of the adjoining composite laminates in order to prevent the occurrence of additional eccentricities stemming from the nodal forces at the interface in the subsequent FE discretization. 1 3 2 Figure 1: Interface ¯ Bdefinition between two adjacent bodies B(i)(i= 1,2) whose response is characterized by a continuous distribution of linear elastic springs. Moreover, differing from alternative formulations, which are expressed using the tractions and relative displacements (usually denominated as gaps) between the faces of the interface as fundamental magnitudes, the current model is derived in terms of the stress and strain tensors. In particular, a reduced (degenerated) version of such tensors, which only account for the interlaminar terms (transverse shear and normal stress– strain relationships), is used [34]. 3 The constitutive formulation of the LEBIM can be conceived as a continuous distribution of linear springs (Fig. 1). For undamaged states at any material point x, the governing law that characterizes the interlaminar stress–strain relations at the interface (which is similar to an uncoupled intrinsic CZM formulation) takes the form:    σ(x) = knεn(x), τ1(x) = kt1εt1(x),for G(x)≤ Gc, τ2(x) = kt2εt2(x), (2) where σidentifies the interlaminar normal traction component, whereas τ1and τ2denote the interlaminar shear traction components; εn,εt1,εt2stand for the interlaminar normal and shear strains, respectively, and kn,kt1and kt2correspond to the normal and shear stiffness properties of the springs distribution. Fig. 2 depcits a schematic representation of the LEBIM for the individual fracture modes I and II and the corresponding mixed mode fracture conditions for 2D cases for the sake of simplicity (defining τeq = q⟨σ(x)⟩2+τ2 1). In Eq. (2), G(x) is the energy release rate, whereas Gcis the fracture toughness for general mixed mode fracture conditions. un‐/reloading un‐/reloading un‐/reloading un‐/reloading a b Figure 2: 2D representation for the LEBIM constitutive law: (a) Individual behavior for mode I (left) and mode II (right). (b) 2D mixed-mode conditions. The previous stress and strain components are referred to a local frame defined at the interface, and therefore the construction of the global formulation requires the computation of standard transformation operations [4]. Furthermore, a particular aspect of the current model is the fact that it only identifies undamaged and fully damaged states without intermediate stages. This is a fundamental difference in com- 4 parison with alternative interface formulations relying on the CZ concept, which usually attain a progressive stiffness degradation up to complete failure [5, 12]. Following [25, 26, 33], the energy release rate (ERR) for each individual direction is given by GI(x) = h⟨σ(x)⟩2 + 2kn ;G2(x) = hτ2 1(x) 2kt1 ;G3(x) = hτ2 2(x) 2kt2 , (3) where GI(x), and G2(x) and G3(x) are the ERRs associated with fracture mode I, and tangential directions t1and t2, respectively, see Fig. 1. ⟨δ⟩±=|δ|±+δ 2represents the positive and negative parts of a number δ. Then, for mixed-mode fracture conditions, the corresponding ERR is given by G(x) = GI(x) + G2(x) + G3(x).(4) Notice that G2(x) and G3(x) not necessarily correspond to pure fracture modes GII (x) and GIII (x), although they coincide in some situations. The critical stress and strain values for each individual direction are, respectively, given by the pairs {σc, εnc},{τ1c, εt1c}and {τ2c, εt2c}. Then, the individual critical fracture energies for the fracture mode I and directions t1and t2, respectively, are given by GIc =h⟨σc(x)⟩2 + 2kn ;G2c=hτ2 1c(x) 2kt1 ;G3c=hτ2 2c(x) 2kt2 . (5) According to the constitutive law described in Eq.(2), any material point at the interface suddenly fails when the energy-based failure criterion is reached, i.e. G(x) = Gc(x). Although there exist several proposals for mixed-mode fracture criteria in 2D, see a review in [26], very few fracture criteria for 3D cases have been developed so far [35, 36, 37]. For comparison purposes with the built-in CZM available in ABAQUS using a finite thickness approach, in the present investigation, the interface failure is evaluated by means of the extended version of the Benzeggah-Kenane failure criterion [35] proposed by Reeder et al. [42], and whose expression is given by Gc=GIc+(GIIc −GIc)G2(x) + G3(x) GI(x) + G2(x) + G3(x)η +(GIIIc −GIIc)G3(x) G2(x) + G3(x) G2(x) + G3(x) GI(x) + G2(x) + G3(x)η (6) where ηis the so-called BK-parameter, which is adjusted using experimental data [5]. In addition to the previous aspects, it is worth mentioning that the current LEBIM does no require the definition of any particular stress-based damage initiation criterion, that governs the initiation of the degradation phase, as is typically employed in CZMs. A final remark worthy to comment is that in the present LEBIM formulation, the energy-based criterion is equivalent to a strength-based criterion as was shown previously in [26] and in some other works cited therein. Thus, under pure mode I fracture conditions, the following relation between GIc,σc,E′′ (apparent or constrained Young’s modulus of the layer) and hshould be always fullfiled: GIc =σ2 ch 2E′′ (7) Thus, depending on the imposed parameters several options can be considered:  The first option renders the experimental characterization of the apparent interface strength σc, and consequently for a given hand E′′ of the adhesive layer GIc is determined.  Alternatively, in the case of GIc is measured experimentally, for a given hand E′′ of the adhesive layer, σcis obtained.  Finally, for scenarios in which both σand GIc are characterized in the lab, along with the elastic properties Eand ν, the present LEBIM works for the layer thickness hverifying the condition h=¯ h, with ¯ h= 2GIcE′′/σ2 cdefining a threshold thickness. 5 A way to release this restriction, for h < ¯ h, is the use of the so-called coupled stress-energy criterion (with independent stress and energy criteria) in the spirit of Finite Fracture Mechanics [43]. Nevertheless, previous numerical solutions of several problems (e.g. the debonding between fibre and matrix in composites) showed that the predictions obtained by the current LEBIM and FFM are in a very good agreement [44]. A thorough discussion of these issues can be found in [26] Section 2.3, falling beyond the scope of the present study. The interpenetration between the faces of the interface after damage is prevented through the use of a penalty condition for negative normal strain values. Therefore, once the interface is broken, the following constitutive relations is considered at any interface point x: σ(x) = kn⟨εn(x)⟩−, τ1(x) = τ2(x) = 0,for G(x)≥ Gc(x),(8) The current model considers frictionless contact conditions regarding the tangential components [38]. It should be also noticed that the LEBIM can be considered as a brittle non-smooth limit of a CZM and its use is confined to configurations where the size of the process/softening zone is small in comparison to the crack size. In fact, it was shown in [45, 46], that the effect of the cohesive law profile (i.e. the particular form of the TSL) plays a minor role in practical cases, being the major discrepancies between the different options in terms of the maixum peak load for the corresponding analysis. 2.2. Comparison between 3D formulations of CZM and LEBIM This section addresses the principal differences between the proposed LEBIM and standard CZ formulations based on an intrinsic approach. In order to facilitate the discussion and to identify possible common attributes of both modeling strategies, the subsequent analysis is focused on pure mode I fracture conditions using a CZM with a bilinear TSL [5, 22], and complying with an interface response in terms of the interlaminar stress-strain relationship (i.e. an interface of initial finite thickness). With respect to the bilinear CZM, Fig. 3.a shows the characteristic stress-strain behaviour for 2D cohesive formulation obeying a bilinear TSL in terms of the interlaminar stress–strain relationship. For pure mode I, damage initiation is related to the interface strength, σc, whereas the area under the curve is proportional to the fracture toughness, GIc, and εnc identifies the critical strain for complete decohesion. Prior damage initiation, an initial elastic phase is assumed, which is characterized by the penalty stiffness parameter ¯ kn. This parameter ensures the stiff connection between the interface faces before damage onset. The interfacial strength σcand the initial stiffness ¯ kndefine the normal strain for damage initiation εn0. Standard CZMs are also equipped with an internal damage variable d∈[0,1], which is defined to track the progressive stiffness degradation during the decohesion process at the interface. Particularly, d= 0 represents a pristine condition while d= 1 corresponds to a fully damaged condition. This variable can be related to the dissipated energy up to the current stage Gd(Fig. 4.a) at any intermediate stage. A similar behavior also characterizes the fracture mode II. Finally, this cohesive interface model can incorporate a particular failure criterion for mixed-mode fracture conditions as it is proposed in [39]. Analyzing the previous description and performing a direct comparison with the constitutive formulation corresponding to the LEBIM (Section 2.1), the following principal differences between both modeling approaches can be identified: 1. With reference to the number of material parameters that are required as input data to characterize the interface response, it is noticeable that the bilinear CZM approach requires a higher number of parameters ({σc,¯ kn,GIc}) than the LEBIM ({σc,GIc}) regardless the application. Indeed, for the particular case of the CZM, usually, there is no direct physical interpretation for the initial stiffness whose initial (very high) value could lead to numerical issues during the numerical solution procedure, whereas σcand GIc values can be obtained by standardized tests. For example, in the case of an adhesive layer Double Cantilever Beam (DCB) or Climbing Drum Peel (CDP) tests [25, 47] can be used. 6 loading un‐/reloading loading loading un‐/reloading loading a b Figure 3: 2D representation for the bilinear CZM: (a) Individual behavior for mode I (left) and mode II (right). (b) Mixed-mode conditions. 2. The CZM provides information with regard to the progress of the interface failure along the simulation due to the fact that the corresponding internal damage variable continuously increases from 0 (no damage) to 1 (complete decohesion) in order to preserve its thermodynamic consistency. Therefore, the graphical representation of this variable allows the identification of critical regions at the interface and the prospective interfacial damage evolution path. In contrast to this, the LEBIM is only able to distinguish between the two extreme cases regarding the state of the interface, i.e. either intact or fully damaged states. From the conceptual point of view, whereas the fracture process zone (FPZ) associated with the LEBIM is equal to zero [40], the CZM has an explicit FPZ definition, which is characterized by the characteristic length lCZM. Then, the gradual appearance of a FPZ ahead the crack tip can be accomplished by means of setting the value of the so-called ductility parameter of the interface ¯µ=εnc/εn0[40], see Fig. 4.b. Thus, for given σcand GIC , the case ¯µ= 1 corresponds to the LEBIM, while by increasing this parameter CZ formulation is recovered, see Fig. 4.b for an schematic representation. 3. In line with the previous analysis, ¯µcan be interpreted as a measure of the characteristic length of the FPZ ahead the crack tip. The specific definition of the characteristic length lcz obeys the expression: lcz =αGcE σ2 c ,αbeing a proportionality scalar parameter which also depends on the initial stiffness value of the interface, see Sect. 3.1. This lack of the nonlinear character of the FPZ in the LEBIM leads to: (1) a notable mitigation of mesh requirements ahead the crack tip (typically 7-10 elements are recommended within this region in standard CZMs), since no gradual transition from the pristine 7 loading un-/reloading loading Figure 4: (a) Intermediate degradation stage for bilinear mode I fracture mode. (b) Bilinear cohesive laws for different values of the interface ductility parameter ¯µ. to fully broken states is required, and (2) a clear identification of the predicted crack front based on the status of the corresponding failure map. Interestingly, based on the previous considerations, it is worth noting that an interface response equivalent to the LEBIM could theoretically be achieved by the use of the CZM through calibrating the corresponding input parameters. However, this option usually leads to errors in most of the commercial FE codes that include this modeling capability since this parameter selection is considered as inconsistent with the built-in implementation (Sect. 3.1). 2.3. Finite element formulation In this section, the FE formulation of the proposed 3D version of the LEBIM is derived. In particular, the current implementation is carried out by means of the user-defined capability UMAT of the FE package ABAQUS [41]. This subroutine allows the implementation of the current interface model in a simple and straightforward manner due to the fact that the coding tasks are exclusively confined to the constitutive formulation at integration point level. Thus, the rest of the standard attributes that are usually required for the correct implementation of an interface model, concerning the definition of the local interface frame, to which the interface law is referred to, and the kinematic definition of the interface element, are internally performed by means of built-in ABAQUS procedures. However, it is worth mentioning that the correct use of the LEBIM needs the definition of a local material setting in the preprocessing stage in order to explicitly distinguish between the in-plane and out-of-plane actions. The following developments are focused on the interface contribution (Πintf(u)) to the complete functional of the mechanical system (ΠT(u)), which is composed by the interface and the surrounding bulk material 8 (Πbulk(u)), i.e. ΠT(u)=Πbulk(u)+Πintf(u), where uis the displacement field. The interface contribution to the principle of virtual work of the system can be expressed as δΠintf(u, δu) = Z¯ B σ(u) : δε(δu) dΩ, ∀δu∈V,(9) where uand σrepresent the displacement and stress fields, respectively. δu∈Vis an admissible virtual displacement field, and δεdenotes the corresponding virtual interlaminar strain vector field. Following general derivations, the real (u), virtual (δu) and incremental displacement (∆u) vectors are interpolated via standard shape functions: u≈ nn X A=1 NA(ξ)dA=Nd;δu≈ nn X A=1 NA(ξ)δdA=Nδd; ∆u≈ nn X A=1 NA(ξ)∆dA=N∆d, (10) where NA(ξ) (A= 1, . . . , nn, with nn= 8 for first-order elements) denotes the trilinear shape functions defined at the node Ain the parametric space ξ={ξ1, ξ2, ξ3}with −1≤ξi≤1 (i= 1,2,3); dA,δdAand ∆dAidentify the real, virtual and incremental displacement vectors corresponding to the node A, whereas Nand dare the suitable operators that contain the interpolation functions and the nodal displacements of the interface element. Using a vector arrangement, the virtual (δε) and incremental (∆ε) interlaminar strain components are interpolated at the element level as δε=  δεn δεt1 δεt2  ≈ nn X A=1 BAδdA=Bδd; ∆ε=  ∆εn ∆εt1 ∆εt2  ≈ nn X A=1 BA∆dA=B∆d, (11) where Bidentifies a suitable compatibility strain-displacement operator, whose version at nodal level is denoted by BA. With the previous definitions at hand, the discrete interface contribution to the principle of virtual work, Eq.(9), yields δΠintf(d, δd) = δdTZ¯ B BTσdΩ=δdTR, (12) where Ris the so-called residual vector at the element level. The consistent linearization of Eq.(12) in the direction ∆dreads ∆δΠintf(d, δd,∆d) = δdTZ¯ B BT∂σ ∂ddΩ∆d=δdTZ¯ B0 BTCintfBdΩ∆d=δdT[K] ∆d, (13) where Kis the element stiffness matrix, and Cintf denotes the interface constitutive operator, which takes the simple form (corresponding to the non zero entries of the complete 6×6 matrix that should be defined in the FE implementation for 3D cases): Cintf =  kn0 0 0kt10 0 0 kt2  .(14) As was previously addressed, the integration of the volume integrals at the element level is carried out internally by ABAQUS. In particular, standard Gauss integration rule with 2×2×2 integration points is employed to perform such numerical integration. Using the previous integration scheme, no spurious oscillations of the stress field are observed in the computations. The 3D version of the LEBIM herein outlined has been implemented into the FE package ABAQUS through the user subroutine UMAT. The algorithmic treatment of the developed UMAT routine is thoroughly described in [31], therefore further details are here omitted for the sake of brevity. 9 It is worth mentioning that the definition of these local settings, to which the cohesive law is referred to, represents a crucial point due to the fact that most of the existing investigations employed the global Cartesian setting as reference for the cohesive law. This aspect aims at shedding light to the analysis of complex external loadings that lead to mixed mode fracture conditions, especially combining modes II and III. Table 1 summarizes the characteristics of the developed models. Notice that Initiation criterion is not considered in models using LEBIM because initiation of damage is obtained computing the fracture energies (GI,G2and G3) and comparing them with the fracture toughness, obtained for its corresponding fracture mode mixity. In this paper, the extension of the B-K criterion proposed in [42] was implemented in the LEBIM subroutine. Table 1: Models for the adhesive joint between two cylinders subjected to pure torsion. Designation Damage model Coordinate system Initiation criterion CZ-XYZ-MaxS Cohesive zone (CZM) Cartesian Maximum nominal stress CZ-RQZ-MaxS Cohesive zone (CZM) Cylindrical Maximum nominal stress CZ-XYZ-QuaS Cohesive zone (CZM) Cartesian Quadratic nominal stress CZ-RQZ-QuaS Cohesive zone (CZM) Cylindrical Quadratic nominal stress LEBI-XYZ Linear elastic brittle interface (LEBIM) Cartesian - LEBI-RQZ Linear elastic brittle interface (LEBIM) Cylindrical - Fig. 12 depitcs the stress distributions corresponding to the stress component σθz for an applied angle ϕz= 0.002. Analyzing this graph, it can be observed that the numerical predictions using a bilinear CZM were identical to each other, i.e. no dependency on the local setting of the cohesive law was appreciated. However, these predictions strongly contrasted with those referred to the LEBIM, which were also identical between them. In fact, from the qualitative point of view, the traction level for the LEBI-models (25 MPa) was lower than those corresponding to the CZ-models (38 MPa). This aspect is ascribed to the different initial stiffness values between CZ and LEBI models, which were set for both numerical interface models before damage occurrence in their respective definitions. Figure 12: σθz at the adhesive layer just before the damage onset (ϕz= 0.0020). 16 Figure 13: σθz at the adhesive layer for ϕz= 0.0026. Distributions of σθz for ϕz= 0.0026 are shown in Fig. 13. This value of the rotation angle ϕzcorresponds to the maximum load level reached for CZ-XYZ-MaxS model. Two aspects are worth to be mentioned for the results of this model: (i) zones with shear stresses over τc(40 MPa) were estimated to appear (purple areas in the graphs) when the maximum value of σθz reached 50 MPa along lines forming ±45owith the center of the cylinder; and (ii) the simulation concluded before the final prescribed angle due to convergency issues in achieving equilibrium solutions. These results are a direct consequence of the use of the maximum stress criterion for damage initiation combined with the Cartesian coordinate system for the cohesive law. Thus, the numerical interface model considered that the x- and y-axes corresponded to the s- and t-shear directions of the cohesive law, leading to an inconsistent cohesive law along the previous ±45olines with the following interface properties: Gc= 1500 J/m2,τc= 40√2 MPa and Ett = 42 MPa, as is shown in Fig. 14. Figure 14: Cohesive law for the CZ-XYZ-MaxS model along the ±45o. Figs. 15–16 depict the σθz distributions corresponding to two intermediate damaged stages of the adhesive layer: the first one, ϕz= 0.003, with an incipient crack; and the second one, ϕz= 0.005, including a larger cracked zone (crack extension is depicted in white in these graphs). Results for CZ-XYZ-MaxS model corresponding to a modified model in which Ett =Ess = 125 MPa are shown in these figures only for the sake of completeness. Predictions for the other three CZ-models herein considered were very similar to 17 each other. Based on the data represented in these graphs, it can observed that in the case that the polar coordinate system is used for the reference of the cohesive law, CZ-RQZ-MaxS and CZ-RQZ-QuadS models, rand θaxes can be correctly identified with the s- and t-directions of the cohesive law. Nevertheless, we cannot guarantee at all that the ror θ-directions respectively corresponded to the s- and t-axis or viceversa, so the distinction between mode II and III is not clear so far (recalling that GIIc =GIIIc in these cases and this has a minor influence on the results). Figure 15: σθz at the adhesive layer for ϕz= 0.003. Figure 16: σθz at the adhesive layer for ϕz= 0.005. To clarify this idea, Fig. 17 depicts the stress distributions of the component σθz for three twisting angles: ϕz= 0.002, ϕz= 0.003 and ϕz= 0.005, for the CZ-XYZ-QuadS and CZ-RQZ-QuadS models when GIIc = 1000J/m2and GIIIc = 1500J/m2, i.e. different fracture energies for the fracture modes II and III were defined. 18 Figure 17: σθz at the adhesive layer for CZ-XYZ-QuadS and CZ-RQZ-QuadS models, being GIIc = 1000J/m2and GIIIc = 1500J/m2. Note that, there are not any differences between the distributions at the bottom of Fig. 17 with those of Fig. 12, 15 and 16 corresponding to the CZ-RQZ-QuadS model. On the contrary, comparing distributions at the top of Fig. 17 with those of Figs. 12, 15 and 16 corresponding to the CZ-XYZ-QuadS model, significant differences are noticed for the graph corresponding to ϕz= 0.003 and, particularly, for the one referred to ϕz= 0.005. These discrepancies are attributed to the identification of the built-in B-K criterion used by ABAQUS concerning the shear fracture modes. Thus, we were not able to verify whether the stress components σrz and σxz had been associated to mode II and σθz and σyz to mode III or viceversa. Similar conclusions can obtained for the corresponding results using the LEBI-models, see Fig. 18. These results were also computed using a 3D version of the modified B-K criterion according to [42]. Moreover, performing a direct comparison between the results shown in Fig. 17 and 18, it can be observed that LEBI-results are of the same accuracy with respect to those using the bilinear CZMs of ABAQUS but with a formulation notably simpler. Nevertheless, it is interesting to recall that results would be different if any other failure criterion would have been used, and whose analysis is beyond the scope of the current investigation. Finally, we also compare the global torque-angle evolutions for the previous models. Evolutions of the applied torque versus the twisting angle are shown in Fig. 19. The CZM-model evolutions coincided along the linear stage (till ϕz≈0.002 or M≈0.940 KN m), which is a consequence of the stiffness values used, being the same for all of them. The initial evolution of the LEBI-models differ from the previous ones because their stiffness values cannot be imposed ad-hoc, since it is computed from the strength (σc,τtc and τsc) and toughness (GIc,GIIc and GIIIc) properties. A way to remove this restriction on the layer thickness modelled by LEBIM is the use of a coupled stress and energy criteria as that developed within the framework of Finite Fracture Mechanics (FFM) [44]. In addition to the previous considerations, a special attention should be focused on CZM(RθZ)QuadS and LEBIM(RθZ) results. These models possess the local axes of the cohesive law oriented according to the crack front, then the mixed law of fracture becomes a pure mode III law for this problem. Quantitatively, the peak torque were similar for both interface formulations, being equal to 1 kN m for CZM(RθZ)QuadS and 0.949 KN m for LEBIM(RθZ). However, the qualitative difference is much more significant if we analyze the corresponding rotation angle, which corresponded to ϕz= 0.0024 for the CZM(RθZ)QuadS model and ϕz= 0.0032 for LEBIM(RθZ). Thus, although the selection of any of the two computational models herewith discussed for triggering damage events in adhesives could be irrelevant from the strength point of view (for equivalent properties), the stiffness response can be completely different from each other. Specifically, this 19 Figure 18: σθz at the adhesive layer for LEBI-XYZ and LEBI-RQZ models, being GIIc = 1000 J/m2and GIIIc = 1500 J/m2. is of crucial importance in applications where the stiffness of the adhesive layer plays an essential role in the total stiffness of the complete system. Figure 19: Evolutions of the torque versus rotational angle. 4. Conclusions In this paper, two computational models to study damage events in thin adhesive layers were carefully cosidered: (i) one of the most popular intrinsic CZM complying with a bilinear TSL, and (ii) the new elastic brittle interface model (LEBIM) for 3D applications. Both formulations allow the simulations of the onset and progression of crack-like damage along a specific adhesive layer or interface. 20 Special attention was paid to the 3D extension of the damage adhesive law and the criterion for complete damage identification. To address these points, two problems were presented and investigated: the first application was the standard double cantilever beam (DCB) problem, whereas the second case comprised the analysis of a crack evolution in the adhesive layer between two cylinders subjected to torsional loading. Relying on the results presented, it can be concluded that the LEBIM does offer the following attributes:  Relative simplicity for its numerical implementation into FE- or BE-based codes.  A higher robustness in comparison to intrinsic CZMs since it is more permissive with regard to the mesh size due to the lack of FPZ ahead the crack tip without affecting the accuracy level of the simulations. In this regard, standard CZMs need more than 5 elements throughout the cohesive zone for capturing the nonlinear TSL, whereas this condition is relaxed for the case of the LEBIM.  Nonphysical responses are precluded based on the fundamental assumptions of the LEBI formulation. Therefore, it could follows evolutions with snap-through and snap-back jumps in a quasi-static way because the load-displacement solution for a specific crack size is path-independent.  General intrinsic CZMs can effectively reproduce the stiffness of the adhesive layer because these characteristics have to be explicitly defined. In contrast to this, the stiffness properties for the LEBIM are computed from the strength and fracture toughness data which can produce unreal elastic properties (i.e. elastic modulus and Poisson coefficient) or artificial anisotropy. Thus, it is not surprising that evolutions of stress field using either the CZ formulation or the LEBI approach could be notably different from each other even prior to the appearance of damage. Nevertheless, from the strength and crack size point of view, both models reported similar results. Apart from the previous considerations, the present research has thoroughly analysed the role of the local material reference to which the adhesive model is referred to. In particular, when 3D problems are modeled, the adhesive stiffness, strength and fracture toughness properties should be accounted for in a rotated system parallel to the crack front in order to consider fracture modes I, II and III in a consistent manner. In the particular case of CZM, the strength and fracture toughness properties are independent from each other, thus two mixed criteria are required: (i) one for the tractions (maximum stress criterium and quadratic stress criterium were considered in this paper) for triggering damage onset, and (ii) an additional energy-based criterion for the identification of complete failure of the adhesive. On the contrary, in the case of LEBIM strength and toughness are related, and therefore only a mixed fracture criterion was needed. Moreover, in light of the present results, it was observed that, in the case of CZM, the rotation of the cohesive law affects the evolution of any element in the cohesive zone. This fact contrasts with the LEBIM in which only the element just at the crack front is affected by this rotation. Correspondingly, an additional hypothesis should be adopted for the CZM in order to position the local system (for the evaluation of the TSL and the failure criterion) according to a virtual crack front direction for any element (which can be quite far, in element size units, from the actual crack front) in the cohesive zone. This condition also holds for the LEBIM, however the identification of the crack front is significantly simpler. In this concern, it is worth mentioning that, though results here presented focused on ABAQUS, the authors do not know any commercial software that rotates the cohesive law and the corresponding toughness-based damage criterion to a reference system parallel to the crack front. Finally, in line with the last argument, the authors are conducting additional research activities addressing this point. In particular, an enhanced 3D model with a new proposal where the damage criterion is properly rotated and referred to the actual crack front is currently under development. Due to the simplicity of this task in the case of LEBIM, this model has been chosen for this purpose. The corresponding results are beyond the scope of the present paper and will be discussed in a forthcoming publication. Acknowledgments This research was conducted with the support of the Spanish Ministry of Economy and Competitiveness and European Regional Development Fund (Projects MAT2015-71036-P and MAT2015-71309-P) and Junta 21 de Andaluc´ıa and European Social Fund (Project of Excellence No. P12-TEP-1050 and TEP-7093). This funding is gratefully acknowledged. References References [1] Allix O, Corigliano A, (1999) Geometrical and interfacial non-linearities in the analysis of delamination in composites. Int J Solids Struct 36, 2189–2216. [2] Balzani C, Wagner W (2008) An interface element for the simulation of delamination in unidirectional fiber-reinforced composite laminates. Eng Fract Mech 75, 2597–615. [3] Camanho PP, D´avila CG, de Moura MF (2003). Numerical simulation of mixed-mode progressive delamination in composite materials. J Compos Mater. 37, 1415–1438. [4] Reinoso J, Bl´azquez A, T´avara L, Par´ıs F, Arellano C (2016) Damage tolerance of composite runout panels under tensile loading. Compos Part B: Eng 96, 79–93. [5] Tur´on A, Camanho PP, Costa J, D´avila CG (2006) A damage model for the simulation of delamination in advanced composites under variable-mode loading. Mech Mater, 38, 1072–1089. [6] Vigueras G, Sket F, Samaniego C, Wu L, Noels L, Tjahjanto D, Casoni E, Houzeaux G, Makradi A, Molina-Aldareguia JM, Vazquez M, J´erusalem A (2015) An XFEM/CZM implementation for massively parallel simulations of composites fracture, Compos Struct 125, 542–557. [7] Reinoso, J., Paggi, M., Rolfes, R. (2016) A computational framework for the interplay between delamination and wrinkling in functionally graded thermal barrier coatings, Comp Mater Sci 116, 82–95. [8] de-Andr´es A, P´erez JL, Ortiz M (1999) Elastoplastic finite element analysis of three-dimensional fatigue crack growth in aluminium shafts subjected to axial loading. Int J Solids Struct 36, 2231–2258. [9] Corrado M, Paggi M (2015) Nonlinear fracture dynamics of laminates with finite thickness adhesives. Mech Mater 80, 183–192. [10] Iannucci L (2006) Dynamic delamination modelling using interface elements. Comput Struct 84, 1029–1048. [11] Harper PW, Hallett SR (2008) Cohesive zone length in numerical simulations of composite delamination. Eng Fract Mech, 75(16), 4774–4792. [12] Ortiz M, Pandolfi A (1999) Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis. Int J Num Meth Eng 44, 1267–1282. [13] T´avara L, Mantiˇc V, Salvadori A, Gray LJ, Par´ıs F (2013) Cohesive-zone-model formulation and implementation using the symmetric Galerkin boundary element method for homogeneous solids, Comput Mech 5, 535–551. [14] Barenblatt G (1962) The mathematical theory of equilibrium cracks in brittle fracture. Adv Appl Mech 7, 55-129, 1962. [15] Dugdale D (1960) Yielding of steel sheets containing slits. J Mech Phys Solids 8(2), 100–104. [16] Petersson PE (1981) Crack Growth and Development of Fracture Zones in Plain Concrete and Similar Materials Crack Growth and Development of Fracture Zones in Plain Concrete, Report TVBM-1006, Lund Institute of Technology, Lund. [17] D´avila CG, Rose CA, Camanho PP (2009) A procedure for superposing linear cohesive laws to represent multiple damage mechanisms in the fracture, Int J Fract, 158, 211-223. [18] Vodiˇcka R, Mantiˇc V (2017) An energy based formulation of a quasi-static interface damage model with a multilinear cohesive law, Discret Contin Dyn S, 10, 1539–1561. [19] Xu XP, Needleman A (1993) Void nucleation by inclusion debonding in a crystal matrix. Model Simul Mater Sci Eng 1(2), 111–132. [20] Tvergaard V (1990) Effect of fiber debonding in a whisker-reinforced metal, Mat Sci Eng A, 107, 23–40. [21] Dimitri R, Trullo M, De Lorenzis L, Zavarise G (2015) Coupled cohesive zone models for mixed-mode fracture: a comparative study. Eng Fract Mech 148, 145–179. [22] Tur´on A, D´avila CG, Camanho PP, Costa J (2007) An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models, Eng Fract Mech 74(10), 1665–1682. [23] Prandtl L (1933) Ein Gedankenmodell f¨ur den Zerreißvorgang spr¨oder K¨orper (A thought model for the fracture of brittle solids). Zeitschrift f¨r Angenwandte Mathematik und Mechanik 13(2),129–133. [24] Bruno D, Greco F, Lonetti P (2003) A coupled interface-multilayer approach for mixed mode delamination and contact analysis in laminated composites. Int J Solids Struct 40, 7245–7268. [25] T´avara L, Mantiˇc V, Graciani E, Ca˜nas J, Par´ıs F. (2010) Analysis of a crack in a thin adhesive layer between orthotropic materials. An application to composite interlaminar fracture toughness test. Comp Model Eng Sc 58 (3), 247–270. [26] Mantiˇc V, T´avara L, Bl´azquez A, Graciani E, Par´ıs F (2015) A linear elastic - brittle interface model: Application to the onset and propagation of a fibre-matrix interface crack under biaxial transverse loads, Int J Fract 195, 15–38. [27] Benveniste Y, Miloh T (2001) Imperfect soft and stiff interfaces in two-dimensional elasticity. Mech Mater 33, 309–323. [28] Hashin Z (2002) Thin interphase/imperfect interface in elasticity with application to coated fiber composites. J Mech Phys Solids 50, 2509–2537. [29] Lenci S (2001) Analysis of a crack at a weak interface. Int J Fract, 108, 275–290. [30] T´avara L, Mantiˇc V, Graciani E, Par´ıs F (2016) Modelling interfacial debonds in unidirectional fibre-reinforced composites under biaxial transverse loads, Compos Struct 136, 305–312. [31] T´avara L, Reinoso J, Castillo D, Mantiˇc V (2018) Mixed-mode failure of interfaces studied by the 2D Linear Elastic Interface Model: macro- and micro-mechanical finite element applications to composites, Int. J. Adhesion 94:8, 627–656 . 22 [32] M. Paggi, J. Reinoso (2015) An anisotropic large displacement cohesive zone model for fibrillar and crazing of interfaces. Int. J. Solids Struct. 69-70:106–120. [33] T´avara L, Mantiˇc V, Graciani E, Par´ıs F (2011) BEM analysis of crack onset and propagation along fiber-matrix interface under transverse tension using a linear elastic-brittle interface model. Eng Anal Bound Elem 35, 207–222. [34] Reinoso J, Paggi M, Bl´azquez A (2017) A nonlinear finite thickness cohesive interface element for modeling delamination in fibre-reinforced composite laminates, Compos Part B: Eng 109, 116–128. [35] Benzeggagh ML, Kenane M. (1996) Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixedmode bending apparatus. Comp Sci Tech 56(4), 439–49. [36] Banks-Sills L (2015) Interface fracture mechanics: theory and experiment. International Journal of Fracture 191 (1-2), 131–146. [37] Banks-Sills L (2006) Fracture toughness of the +45◦/-45◦interface of a laminate composite. International Journal of Fracture 141, 195–210. [38] Graciani E, Mantiˇc V, Par´ıs F, Bl´azquez A (2005) Weak formulation of axi-symmetric frictionless contact problems with boundary elements: Application to interface cracks. Comput Struct 83, 836–855. [39] Tur´on A, Camanho PP, Costa J, Renart J (2010) Accurate simulation of delamination growth under mixed-mode loading using cohesive elements: definition of interlaminar strengths and elastic stiffness, Compos Struct 92(8), 1857–1864. [40] Dimitri R, Cornetti P, Mantiˇc V, Trullo M, De Lorenzis L (2017) Mode-I debonding of a double cantilever beam: A comparison between cohesive crack modeling anf finite fracture mechanics. Int J Solids Struct. 124, 57–72. [41] Abaqus-Inc. Abaqus user manual, Version 6.9. Dassault Systems Simulia Corp, Providence, RI, USA; 2010. [42] Reeder J,R, Kyongchan S, Chunchu, P,B, Ambur D,R (2002) Postbuckling and Growth of Delaminations in Composite Plates Subjected to Axial Compression, 43rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Denver, Colorado, vol. 1746, p. 10. [43] Cornetti P, Mantiˇc V, Carpinteri A (2012) Finite Fracture Mechanics at elastic interfaces, Int J Solids Struct. 49, 1022– 1032. [44] Mu˜noz-Reja M., T´avara L., Mantiˇc V., Cornetti P. (2016) Crack onset and propagation at fibre-matrix elastic interfaces under biaxial loading using finite fracture mechanics. Compos Part A: Appl Sci Manuf 82, 267–278. [45] Wang JT (2013) Investigating some technical issues on cohesive zone modeling of fracture. J Eng Mater Technol 135, 011003. [46] Cornetti P, Corrado M, de Lorenzis L, Carpinteri A (2015) An analytical cohesive crack modeling approach to the edge debonding failure of FRP-plated beams. Int J Solids Struct 53,92–106. [47] Ca˜nas J., T´avara L., Bl´azquez A., Estefani A., Santacruz G. (2018) A new in situ peeling test for the characterisation of composite bonded joints. Compos Part A: Appl Sci Manuf 113, 298–310. 23