Cohesive/adhesive failure interaction in ductile adhesive joints Part I: A smeared-crack model for cohesive failure

Abstract This paper proposes a new methodology for the finite element (FE) modelling of failure in adhesively bonded joint. Unlike current methods, cohesive and adhesive failures are treated separately. Initial results show the method׳s ability to give accurate prediction of failure of adhesive joints subjected to thickness-induced constraint and complex multi-axial loading using a single set of material parameters. The present paper (part I), focuses on the development of a smeared-crack model for cohesive failure. Model verification and validation are performed comparing the model predictions with experimental data from 3 point bending End Notched Flexure (3ENF) and Double Cantilever Beam (DCB) fracture tests conducted on adhesively bonded composite panels of different adhesive thicknesses.


Introduction
Due to their light weight and high strength and stiffness, composite materials are becoming increasingly popular for the design of many engineering components in areas as diverse as tidal and wind energy, jet-engines, automobiles, etc. Due to the complexity of one-piece component manufacturing, the need to develop good joining techniques is becoming more pressing. As opposed to more widespread joining techniques such as bolted and rivetted joints, adhesive joints do not need holes to be machined, thus reducing the addition of high stress concentrations and the risk of component failure [1]. Furthermore, they are lighter and more economical. Historically, the difficulty to assess their integrity in a fast non-destructive way and to predict their resistance to failure accurately has been a major impediment to a more systematic use of adhesive joints in composites-based engineering components design. However, due to the combined effects of improved adhesive mechanical performance, a better understanding of the failure mechanisms involved and the increased accuracy of the numerical methods available to the designers, this has recently started to change.
The cohesive zone method (CZM) has been the main contributor to the improvement of failure predictions of adhesive joints. The method allows the simulation of damage growth by consideration of energetic principles and allows taking into account phenomena such as mixed-mode loading [2][3][4], rate-dependent effects [5], environmental effects [6,7] and fatigue loading [8]. The mechanical response of the adhesive can either be represented fully with the traction-opening response of the cohesive zone [9][10][11] (this is essentially the strategy used in the modelling of thin adhesive layers) or with a layer of cohesive elements incorporated within a bulk material made of elastoplastic solids [12][13][14] (this suits best the description of failure in thick adhesive layer). In both cases, no formal distinction is made between cohesive (i.e. rupture of the adhesive bulk material) and adhesive failure (i.e. interfacial failure or debonding of the adhesive). In other words, the traction-separation law used in the cohesive elements is set such that it represents the overall behaviour of the bond which results from the interaction between cohesive and adhesive failure.
Failure in brittle adhesives can easily be modelled using a simple bi-linear traction-separation law and one single set of material parameters. On the other hand, failure mechanisms of ductile adhesives involve complex multi-axial plastic deformations. These can be responsible for up to 80% of the load carrying capacity of common structural adhesives [15] and can result in the non-monotonic dependency of bond toughness and strength with joint geometry [16]. In such a case, accurate failure predictions necessitate varying the traction-separation law parameters [14] and shape [17] with the bond dimensions. As a result, full characterisation of a ductile adhesive can be both costly and time consuming. This is particularly the case since adhesive joint behaviour can also depend on hydrostatic pressure [18,19], adding an extra level of complexity. Finally as pointed out in [14], the method assumes the crack path a-priori. As the position of the crack path in the layer affects the plastic strain distribution the contribution due to plastic yielding in the adhesive (which, in return, acts upon the traction-separation law definition) is an additional potential source of inaccuracy of the state-of-the-art modelling technique of adhesively bonded joints.
Here, a new methodology whereby cohesive and adhesive failures are treated separately is proposed. The adhesive deformation (i.e. cohesive failure) is modelled via a smeared-crack approach which allows crack propagation without knowing the crack path in advance and which can easily be used to model phenomena where plasticity and damage coexist [20]. Preliminary analytical calculations (using the method described in [21,22]) and CZM-based finite element analysis (that follows the method described in [14]) made for the predictions of the failure load of the double-lap joint specimens presented in a companion paper [23] suggested that the adhesive plasticity and its dependence to hydrostatic pressure [18,19,24,25] play a key role in the failure mechanisms of the structural adhesive used in the present study. Based on these observations, the adhesive plasticity is taken into account through the use of a pressure-dependent yield criterion. Adhesive failure, on the other hand, is modelled by inserting a layer of cohesive elements at the interface between the adhesive and the adherends. The present paper (part I) focuses on the formulation of a new model for cohesive failure. Particular attention is given to the influence of the bondline thickness and loading mode on the joint apparent toughness. In the companion paper [23], the proposed methodology is applied to the modelling of a double lap-joint specimen with dissimilar adherends and the interaction between cohesive and adhesive failure is studied in more detail.

Model formulation
Both continuum damage mechanics (CDM) and cohesive zone modelling make use of history variable that tracks the extent of damage accumulated on the crack plane. d is a damage parameter whose initial value is d 0 ¼ 0 and that remains zero until the damage initiation condition is met. It is a monotonically increasing value, which reaches the failure value 1 when the crack faces fully separate. In CDM, the evolution of d is closely tight to the evolution of the strain ε at the material point under consideration.
Whilst in cohesive zone modelling, the variation of d is dictated by the shape of the traction-separation law. The smeared crack model approach [20,26] used here allows to relate the damage variable d traditionally used in CDM, to the damage variable that would be used to describe an equivalent traction-separation law within a cohesive zone modelling framework. Smeared-crack models are based on the decomposition (see Eq. (1)) of the total strain tensor (ε) into the sum of a part solely due to the material deformation (ε e þ ε p ) and another term,ε c , accounting for the cracking contribution.
In Eq. (1), the part of the total strain tensor due to the material deformation is additively split into an elastic term (ε e ) and a plastic term (ε p ). The additive splitting of the strain used here implies some limitation on the validity of the model concerning large strains. Giving the relatively large strains observed in the adhesive studied (see Fig. 2 for example), this is an evident limitation of the proposed approach. However, due to the complexity of the physical phenomena considered, some simplifications had to be made. The adaptation of the model to a large-deformation framework could possibly be implemented in future. It is also worth noting that this approach is consistent with the small strain assumptions made by a number of other contributions in the literature investigating the failure of ductile adhesive bonds (e.g. [9,14]).
Here, 8 noded brick elements are used. The material degradation is triggered when the plastic strain exceeds a certain threshold. The evolution of a damage parameter (d) is determined from energetic principles via a traction-separation law. d controls the norm of an anisotropic second-order damage tensor, D, that allows the localisation of damage into a plane.

Plastic deformation
The aim of the present contribution is to set a new methodology for failure prediction in adhesive joints and to gain understanding of the way the adhesive plasticity, cohesive and adhesive failure interact with each other. 1 Hence a very pragmatic approach is taken. Prior to degradation, the influence of hydrostatic pressure on the development of plasticity in the adhesive was modelled through a simple linear Drucker-Prager (D.P.) yield criterion (see Eq. (2)).
σ VM and σ H are the Von Mises stress and hydrostatic pressure respectively. ξ, c and η are parameters that can be expressed as a function of the adhesive's yield stresses in pure tension (σ T ) and in pure compression (σ C ).
The evolution equation for the plastic strain is: with n being the direction of the plastic flow and _ γ the plastic multiplier. Inserting Eq. (3) into the Hooke's law and using the additive splitting describe in Eq. (1) (ε c is assumed to be null at the moment) gives rise to: where C e is the elasticity tensor. The evolution equation for the equivalent plastic strain is then given by: The loading and unloading conditions are given by the Kuhn-Tucker equations: In case of plastic loading _ γ 40 and ψ ¼ 0 holds whilst in the case of unloading, _ γ ¼ 0 and ψ o 0. This can be summarized with the consistency condition: In the case when yielding occurs, _ γ can be obtained from the consistency condition: which using (Eqs. (4) and 5) can be rearranged as: 1 The interaction between cohesive and adhesive failure is addressed in the part 2 paper [23].
Rearranging Eq. (9) and inserting the result into Eq. (4) finally gives rise to:

Damage onset and propagation
In cohesive zone modelling, only the three components of the relative displacement acting in the plane where the crack propagates are needed to describe the traction-separation law. Some mixed mode element formulations [27], resort to an effective mixed-mode relative displacement defined as δ m ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi (where δ I and δ II are respectively the normal and shear components of the relative displacement for a pair of interface connecting points in a local orthogonal co-ordinate system) and an effective stress measure defined a σ m ¼ . The degradation of the crack plane properties is performed in the effective space (δ m ,σ m ). In the proposed framework, a continuum formulation is used and the 9 components of any tensors (e.g. the stress σ, or strain ε, etc.) are considered. A similar strategy to mixed mode cohesive zone modelling is used. The deformation state of an element is reduced to the computation of a single scalar value, the "Von Mises" relative displacement, δ VM defined in Eq. (11) where l c is the element length in the direction perpendicular to the plane where the crack propagates and ε VM2 is the Von Mises effective strain given in Eq. (12) where ε D is the deviatoric part of the strain tensor defined as: As shown in Eq. (2), the Drucker-Prager criterion is a modification of the Von Mises criterion, whereby the dependence of the material with the hydrostatic pressure is taken into account by making the Von Mises stress dependent on the hydrostatic stress. Therefore, the use of Von Mises relative displacements and stress is consistent with the choice of yield surface.

Traction-separation law
As illustrated in Fig. 1, the material degradation is performed, in a newly defined effective space formed by σ VM and δ VM , through a traction-separation law which controls the evolution of a damage variable d. σ eff sets for the effective stress tensor. The Von Mises stress σ VM is computed by degrading the Von Mises equivalent effective stress (σ VM eff ) via Eq. (13). Damage initiates when σ VM reaches the material tensile strength σ VM max or equivalently when the equivalent plastic strain, ε p eq , reaches a certain threshold, ε p eq i. At that point, the cohesive displacement takes the value δ VM i and the energy dissipated is G i . After further loading of the material, its strength is assumed to degrade linearly until it fails completely when the energy G dissipated during the process (i.e. the area under the curve σ VM ¼ f ðδ VM Þ in Fig. 1) reaches the critical energy G C , defined as the sum of the stored elastic energy (released when the crack grows) and a dissipation term embracing the plastic dissipation and the energy necessary to create 2 new surfaces [28]. As the constituents of G C are extremely tied to each other, they cannot be separated. Therefore, we worked with G C rather than with each of its constituents individually. The Von Mises cohesive displacement at complete failure (δ VM f ) is determined from Eq. (16).
The parameter m in Eq. (15) simply enforces the linear evolution of the traction separation law once the damage initiation criterion has been met (i.e. for δ VM 4 δ VM i ). It corrects for the fact that the Von Mises effective stress continues to increase after passing this point. By doing so, a small degree of coupling between damage and plasticity (which accounts for the interaction of microcracks and microvoids with the plastic flow) is introduced and damage evolves nonlinearly with respect to δ VM .

Multiaxial loading initiation and propagation criteria
As alluded to in the introduction, the dependency of the adhesive plastic behaviour on the hydrostatic pressure implies that ε p eq i and G C vary with σ H . In order to reduce the number of material parameters, two criteria for crack initiation and propagation are proposed. As illustrated in Fig. 2, Eq. (17) assumes that the plastic strain at damage initiation varies linearly with the adhesive yield stress (σ Y ). Despite being quite a simplified assumption this matches existing data collected in [29] fairly well. Furthermore, on the basis that most of the energy is dissipated between the onset of plasticity and the onset of damage and that the stress stays (almost) constant between both events, it is then possible to derive the criterion proposed in Eq. (18). Although fairly simple, Eq. (18) builds up on the work by Li et al. [30] who tested several compressive enhancement criterion for the study of delamination in composite material using interface elements and showed that this kind of criteria were the most adapted to the case.
These two criteria allow the determination of ε p eq i and G c irrespective of the loading mode using the values (ε p eq i 1 and ε p eq i 2 ) taken by ε p eq i in two distinct load configuration and a measure of Fig. 1. Traction-separation law for the new smeared-crack model for cohesive failure in adhesive joints. 2 Due to the nature of failure in polymers (i.e. the polymer chains extend to some distance from the fracture plane), the plastic strain term and the term accounting for the separation of two surfaces in Eq. (1) are closely tied to each other and difficult to separate. It was therefore chosen here to work with the total strain tensor rather than with each of its constituents separately.
in one loading configuration taken as reference.
where a is the slope of the curve ε p were measured under pure tension and pure shear respectively (see Fig. 2) and the reference configuration was set so that σ ref

Stress tensor update procedure
The degradation of the stress tensor can either be isotropic or anisotropic. Isotropic damage assumes that the stiffness degradation, which occurs along with damage development and crack propagation, is the same in every direction and is independent of the loading direction. This assumption is accurate early in the deformation process. However when voids start to coalesce and cracks (that are by nature localised onto one single plane) start to form and propagate, other degradation schemes need to be considered. This has motivated the development of anisotropic damage theories. In the present study, the Desmorat and Cantournet [31] pseudo-scalar damage approach is used. As the creation of new surfaces in ductile media is intimately tied to yield, it is assumed that the principal directions of the anisotropic damage rate ( _ D) coincide with those of the plastic strain rate tensor ( _ ε p ). A scalar damage variable, d, is first calculated using the method described in Section 2.1 and it is then postulated that _ D follows the equation: The notation •j j applied to a tensor A indicates the operation consisting of: (i) diagonalising the tensor A through the change of base matrix P, A diag ¼ P À 1 AP (ii) taking the absolute value of the diagonals terms defining A diag j (iii) turning back the tensor in its initial base A ¼ P A diag P À 1 . The notation max ð Þ applied to a tensor A leads to a scalar that is the maximum of the terms of A diag . The stress tensor is then defined as: The unit tensor, 1, is composed of the column of the three unit vectors. The notation • ð Þ 1 2 applied to a tensor A indicates the operation consisting in: (i) diagonalising the tensor A through the change of base matrix P, A diag ¼ P À 1 AP (ii) taking the square root of the diagonals terms defining A diag 1 2 (iii) turning back the tensor in its initial base A ¼ PA diag 1 2 P À 1 .

Model implementation
The proposed model has been implemented in the form of 8noded three-dimensional elements with a single integration point in the commercial explicit finite element (FE) software LS-Dyna, via the addition of a user-defined material subroutine for the under-integrated hexahedral continuum element. Explicit finite element codes are particularly well suited for the description of impact problems and are therefore also well adapted for the description of highly non-linear damage phenomena that can occur quite suddenly. In order to overcome the traditionally slow running times inherent to the use of an explicit solver, mass was scaled by a factor of 10 5 . This resulted in reduced computational times whilst still avoiding significant dynamic effects.
The hardening curves σ C ¼ f ðε p eq Þ and σ T ¼ f ðε p eq Þ as well as the adhesive critical energy, G ref C , are given as input to the model. At every time increment, t, and every integration point, the values for the strain tensor (ε) and for the cohesive displacement, δ VM , are updated. The algorithm for Drucker-Prager piecewise linear plasticity presented in [32] is then used to update the effective stress and plastic strain tensors (σ eff and ε p ) as well as the effective Von Mises stress (σ VM eff ) and the equivalent plastic strain (ε p eq ). The procedures described in Sections 2.2.1 and 2.2.2 are then used to update the value of the scalar damage variable using Eq. (14). The components of the anisotropic damage tensor (D) are updated using (Eqs. (19) and 20): The computation of the principal values and directions of the plastic strain rate tensor allow determining of ΔD. The current values for the stress tensor (σ) are updated using Eq. (20). In order to avoid having to prescribe crack face contacts, once fully damaged (d ¼ 1), the elements are not deleted. They are prescribed a somewhat artificial residual stiffness of 0.03% that allows the analysis to remain stable.

Model verification and validation
In this section, the validity of the smeared-crack model for cohesive failure proposed in Section 2 is assessed. 3ENF and DCB tests were performed on adhesively bonded composites panels of varying adhesive thicknesses. The 3ENF test data were used to inversely determine the critical fracture energy under pure shear (G ref C ) of the adhesive under consideration. Model validation was assessed by comparing the model predictions with the experimental data collected from the DCB tests.

The adhesive
The adhesive used here is a commercial high strength epoxy system notably used for the bonding of large structures such as yacht hulls and wind turbine blades. It is often chosen due to its light weight, its good gap filling properties and the relatively low degradation of its mechanical properties when subjected to environmental factors such as moisture and temperature. As illustrated in Fig. 2, the adhesive is highly ductile in shear and a lot more brittle in tension. The hardening curves (where the Von Mises Yield stress is plotted against the equivalent plastic strain) for pure tension and pure shear (in blue and red respectively in Fig. 2) have been derived from experimental testing on butt joints under tensile [29] and torsion loading [33] respectively. The curve corresponding to the purely compressive stress state has been extrapolated from the tension and shear curves. The manufacturer's data gives the adhesive Young's modulus as 2300 MPa and its Poisson ratio as 0.33.

Experiments
3ENF [34] and DCB [35] fracture tests were performed on specimens made of 3.937 mm-thick adhesively bonded aerospace grade carbon-reinforced composite (IM7/8552) panels. Adhesive layers of 0.25 mm ("very thin"), 0.5 mm, 1 mm and 1.5 mm nominal thicknesses were considered. The specimens were made very rigid, in order to stay as much as possible in the range where LEFM assumption is valid. However the obtained values are purely indicative (this is particularly true for the mode I data) as the nature of the adhesive is such that its deformation is expected to be significant. The obtained data are purely used for model verification and validation purposes. 3 The specimen geometry is presented in Fig. 3. In 3 of the 4 batches, the bondline thickness was controlled by inserting 0.5 mm, 1 mm and 1.5 mm thick PTFE sheets between the composite panels as a pre-crack ahead of the crack tip between specimens to act as spacers. At the end of the  PTFE sheets forming the pre-crack, very sharp crack tips were created by inserting a 20 mm long strip of 0.02 mm thick ETFE release film. The remaining batch was manufactured using only the ETFE release film in the pre-crack region between the composite panels in order to achieve a "very thin" adhesive layer. To ensure good bonding [36], the composite panels' surfaces were carefully cleaned with acetone and abraded using plastic mesh scouring pads taking care not to damage the fibres of the composite prior to the adhesive application. To ensure removal of loose particles and degreasing of the material, a second wash with acetone was then performed. For both series of tests, a universal Instron s testing machine with a 25 kN load cell was used in displacement control. The cross-head displacement speed was set to 0.25 mm/min. 8 3ENF fracture tests were performed for each bondline thickness (Fig. 4). The PTFE sheets were kept between the faces of the initial crack to ensure pure mode II crack propagation and to avoid unexpected effects due to crack face contacts. As illustrated in Fig. 4, the procedure described in [34] was followed. The lower roller span was 100 mm. The mode II fracture toughness (G II C ) was determined using the compliance calibration method. The average value of G II C as a function of the adhesive thickness is plotted in Fig. 4b where the error bars represent the standard deviations. The results obtained here are in the same order of magnitude as other data published in the literature for thick, ductile adhesives [37]. The graph shows typical ductile adhesive joint behaviour characterised by an increase of G II C with the bondline thickness [16] at lower thicknesses.
As illustrated in Fig. 5, for each thickness 5 DCB tests were performed. The initial crack length (a 0 ) was fixed to 60 mm. The average values of G I C (derived from this series of tests) are plotted in Fig. 5b (the error bars represent the standard deviations) as a function of the adhesive thickness. The increasing G I C with increasing bondline thickness suggests the existence of some ductility. Fig. 6, shows the crack surfaces of some of the tested DCB specimens and highlights the competition between adhesive failure and cohesive failure. A qualitative assessment suggests that, the thinner the specimens, the more they are prone to adhesive failure. Hence, whilst all but one of the "very thin" specimens failed by debonding of the adhesive, all of the 1.5 mm-thick specimen showed failure in the adhesive.

G ref C calibration
In the model formulation, G C (similarly G ref C ) applies to an idealised infinitely thin bond (i.e. all the dissipation phenomenon such as plasticity, damage, etchappening in the vicinity of the crack plane are excluded from the determination of the fracture energy terms). As a result, G C can only be determined through an inverse method.
Pure shear was chosen as the reference loading configuration (see Section 2.2.2). Extrapolating the curve on Fig. 4b) to an adhesive joint of zero thickness, G ref C was found to be in the order of 4000-5000 N/m. As this extrapolation is questionable (a linear function is fitted through the experimental points for the specimens with the 2 lower adhesive thicknesses), the obtained value was subsequently refined using an inverse method.
FE models (see Fig. 7) for 3ENF fracture tests performed on specimens with different bondline thicknesses (i.e. 0.10 mm, 0.25 mm, 0.75 mm, 1.00 mm, 1.25 mm, 1.50 mm, 1.75 mm and 2.00 mm) were set up. The unidirectional IM7/8552 composites panels were modelled using solid elements and an orthotropicelastic material (see properties in Table 1). The PTFE sheets were represented by solid elements. An elastic material, with a Young's modulus of 580 MPa and a Poisson ratio of 0.46 (widely accepted values) were used.
The adhesive was modelled using the new smeared-crack model for cohesive failure of adhesives and the material parameters given in Section 3.1. No account is taken of possible adhesive debonding. Quasi static analysis under displacement control was performed. The bending rig rollers were modelled as rigid bodies. Frictionless contact surfaces were included to prevent penetration of the initial crack surfaces and the PTFE sheet. Tied contacts were used to connect the adhesive with the IM7/8552 panels whilst having a refined mesh in the joint (0.1 mm*0.1 mm*2.86 mm brick elements were used), which is necessary in order to have enough elements in the damaged zone and allow mesh-independent solutions (see [38]). For each adhesive thicknesses, the FE models of the 3ENF tests were run with a range of G ref C values. The corresponding mode II fracture toughnesses were calculated (see Section 3.2) and Fig. 8 plots their evolution with the bondline thickness. The value of G ref C giving the best fit with experimental data was found to be 3.3.2. Non-monotonic variation of the apparent toughness with the adhesive thickness As illustrated in Fig. 8, the model predicts a non-monotonic evolution of the apparent fracture toughness with the adhesive thickness. Figs. 9-11 provide more insight into the model behaviour and the physical phenomena that can explain the existence of this kind of evolution. Fig. 9, where the plastic strain distribution in the adhesive is plotted for 3 different values of the adhesive thickness (i.e. 1 mm, 1.25 mm and 2 mm), illustrates the presence of constraint effects [10,14,[39][40][41]. For the three adhesive thicknesses considered, the 2D cross-section of the region where the plastic strain is maximum was assumed to be elliptical in shape, the area of which was calculated. It is shown that for the 1 mm-thick adhesive, plastic deformation occurs within an area of approximately 22 mm 2 . For the 1.25 mm-and the 2 mm-thick adhesives this value decreases to 17.3 mm 2 and 15.3 mm 2 respectively. Therefore it can be assumed that the local maximum reached in the apparent toughness when the adhesive is approximately 1 mm-thick can be mostly explained using a traditional [16,42] constraint effect argument: When the joint is thinner than the size of the plastic zone in the bulk adhesive material, the development of the plastic zone ahead of the crack tip is altered by the thickness of the joint. As the adhesive thickness increases, more and more energy can be dissipated through plastic straining resulting in an increased toughness. When the adhesive thickness is about the same as the size as the plastic zone in the bulk material, the interaction of the plastic zone with the adhesive thickness is responsible for a slightly thinner but much longer (in the crack propagation direction) plastic zone in comparison to what it would be in the bulk material. This can be observed in Fig. 8 in the case of the 1 mmthick adhesive. In this case, more energy can be dissipated through plastic straining and the toughness reaches a peak value. If the adhesive thickness is increased further, the effect of the bondline thickness is progressively relaxed and the plastic zone is less constrained. It progressively becomes shorter and deeper (see the the 1.25 mm-and the 2 mm-thick adhesives in Fig. 8) and its size starts to decrease. For a thick enough adhesive, the plastic zone would develop as in the adhesive bulk material.
From the discussion in the previous paragraph it could be expected that when the adhesive becomes thicker than 1.25 mm the apparent toughness (see Fig. 8) should plateau around the toughness of the adhesive bulk material. Instead it is predicted to increase again. Fig. 10, gives some understanding of the model behaviour causing this prediction. The figure shows tractionseparation curves in the "Von Mises" space (see Fig. 1) extracted from elements situated at the crack tip in FE models of specimens with 1 mm-, 1.25 mm-and 2 mm-thick bonds. As shown on the curves, the model predicts that an increase in the adhesive thickness will result in an increase of the energy needed to fully fail the material (which is represented by the area under the traction-separation curve). This is a direct consequence of the using Eq. (2) and (Eqs. (17) and 18). The rise in the adhesive's thickness results in an increase in the mode I (compressive) component at the crack tip. Due to the Drucker-Prager criterion (Eq. (2)), the yield stress rises with this increased level of compression. (Eqs. (17) and 18) are then responsible for the increase of the material toughness as its yield stress increases. A similar but much pronounced evolution of the fracture toughness with the mode mixity was observed for rubbers by Giannis et al. [43].
In summary, the predicted non-monotonic evolution of the fracture toughness with the adhesive thickness (see Fig. 8) is the result of a trade-off between the constraint effects and the change of mode mixity that occurs with increasing adhesive thickness (see Fig. 11). For adhesive thickness (t) lower than 1 mm, both constraint effects and mode mixity lead to an increase in toughness with thickness. For 1 mm o t o 1:3 mm, the relaxation of the constraint effect exerted by the bondline thickness on the plastic zone results in the drop in predicted apparent toughness. Finally, for t 4 1:3mm, the constraint effect is completely relaxed. The increased compressive stresses that occur with an increase of t cause the model to predict an associated increase in bond toughness.
A similar variation in the fracture toughness with adhesive thickness has been observed in other literature data, for example for the mode I toughness of toughened epoxy adhesives [16]. An additional experimental point would however be needed to fully validate this in future work.

Model validation against mode I fracture tests
Similar FE models to those in the previous section were set up for the DCB fracture tests. In order to ensure mesh-independent results, the mesh size in the adhesive was set to 0.0625 mm In Fig. 12, model predictions for the evolution of G I C with the adhesive thickness are superimposed on the experimental data. A reasonably good correlation was obtained for higher thicknesses (less than 10% difference from each data point). Furthermore, even though they were relatively far from the mean test values, model predictions stayed within the experimental error range at lower thickness. It is important to note here that the present model aims at capturing the failure response within the adhesive and that the FE models do not account for possible debonding that may arise at the interfaces between the adhesive and the composite panels. The observed discrepancy between the model prediction and experimental results in the thickness range where extensive adhesive debonding occurred is therefore not unexpected. It can be explained by the interfacial failure mechanism present in most of the low bondline thickness specimens, which was thought to be responsible for reduced adhesive fracture toughness compared to the cohesive failure only solution. As a result, when the model predictions are compared to the test data of samples breaking under clear cohesive failure only (see blue squares in Fig. 12) a much better agreement between model predictions and experimental results is obtained. Fig. 13 shows the model's ability to reproduce the very brittle nature of failure in the DCB specimens tested. Here, only the results for 1 mm thick joint specimens are presented but similar agreement was also obtained for the 1.5 mm thick joint samples.

Conclusions
In this paper an alternative method to cohesive elements for fracture prediction in ductile adhesive materials has been    Adhesive thickness (mm) Fig. 11. A trade-off between constraint effects (that can be quantitatively assessed by the plastic zone size A) and the increased compressive stresses observed in thicker joints explains the non-monotonic evolution of the apparent toughness with the bondline thickness predicted by the model. At t¼ 0, the size of the plastic zone is zero and crack tip is under pure shear loading (i.e. the Von Mises stresses at the crack tip is equal to the adhesive yield stress in pure shear).
developed. The new smeared-crack model holds the same versatility as the cohesive zone method and could be easily adapted to phenomena such as: rate-dependent and environmental effects, microstructure gradient, or fatigue loading. In the case of very ductile adhesives (such as the one studied here) the adhesive plasticity cannot be disregarded as it is one of the main mechanisms through which the material dissipates energy whilst deforming. Taking account of this within a CZM framework traditionally used for the modelling of fracture in adhesive bonds requires the use of a trapezoidal cohesive law [14]. This makes the material characterisation and testing required fairly onerous. This characterisation is further complicated when the variation of the material properties with hydrostatic pressure needs to be considered (as it is the case in the present study).
Moreover, in the case where a layer of cohesive elements is introduced within a bulk of elasto-plastic solid elements, plasticity in the adhesive (hardening) ends up counteracting the damage development in the cohesive elements (softening) [44]. The solution to overcome this problem is to "tweak" the cohesive law so that the maximum traction in the interface never exceeds the stress levels in the elasto-plastic solids. This not only introduces some degree of inaccuracy but also makes it more difficult to consider different levels of hydrostatic pressure as it creates the need for variable reduction of the maximum traction in the cohesive law with the hydrostatic pressure.
In the framework proposed here, however, the only material characterisation required is the stress vs strain curves of the adhesive under two different loading modes (e.g. pure tension and pure shear) and the intrinsic work of fracture resulting from the separation processes in the fracture process zone (G C ) in one loading mode. As shown, G C can be determined from an inverse method on ENF samples with different bondline thicknesses. Due to the extreme complexity of the problem, the experimental determination of G C with mode mixity, constraints and possible Rcurve effects, remains very challenging. Therefore a very pragmatic approach was used: the evolution of G C with bondline thickness in the experiments was estimated using only basic LEFM assumptions, not strictly applicable for the case here with significant nonlinearity; to compare like with like, the same data reduction procedure was used on simulated load vs cross-head displacement curves generated by the FE models in which different values of G C were hypothesised. It was concluded that the appropriate value for G C was the one for which the evolution of G C with the bondline thickness predicted by the model matched the one obtained experimentally. If the model were to be applied to a test with parameters significantly different from those used for this calibration, it may be that this procedure is not sufficiently robust to cope with all eventualities. Using the newly defined approach, it has been possible to predict the effect of constraint on the fracture toughness of adhesively bonded composite ENF and DCB specimens with different adhesive thicknesses using one single set of parameters. approach was used. LEFM data reduction procedure approach was used The present study focused on cohesive failure (i.e. within the adhesive layer) of adhesive joints. In a companion part 2 paper, the competition between adhesive and cohesive failure which is very important for adhesive joint failure (as demonstrated in the DCB test cases) is explored in more detail.
Von Mises stress field at crack tip: