Failure load prediction of adhesively bonded single lap joints by using various FEM techniques

Journal Pre-proof Failure load prediction of adhesively bonded single lap joints by using various FEM techniques M.Z. Sadeghi, A. Gabener, J. Zimmermann, K. Savarana, J. Weiland, U. Reisgen, K.U. Schroeder PII:

S0143-7496(19)30242-8

DOI:

https://doi.org/10.1016/j.ijadhadh.2019.102493

Reference:

JAAD 102493

To appear in:

International Journal of Adhesion and Adhesives

Please cite this article as: Sadeghi MZ, Gabener A, Zimmermann J, Savarana K, Weiland J, Reisgen U, Schroeder KU, Failure load prediction of adhesively bonded single lap joints by using various FEM techniques, International Journal of Adhesion and Adhesives, https://doi.org/10.1016/ j.ijadhadh.2019.102493. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Elsevier Ltd. All rights reserved.

Failure load prediction of adhesively bonded single lap joints by using various FEM techniques M.Z. Sadeghi1,*, A. Gabener2, J. Zimmermann1, K. Savarana1, J. Weiland2, U. Reisgen2, K. U. Schroeder1 1

Institute of Structural Mechanics and Lightweight Design, RWTH Aachen University, Wuellnerstrasse 7, 52062 Aachen, Germany 2

Saint-Gobain Performance Plastics Pampus GmbH, Am Nordkanal 37, D-47877, Willich, Germany

3

Welding and Joining Institute, RWTH Aachen University, Pontstrasse 49, 52068, Aachen, Germany

Abstract: Fracture toughness, i.e., GC, is considered an important fracture characteristic parameter by which the state of integrity of the adhesively bonded joints under different loading modes can be assessed. Various finite element (FE) methods are capable of modelling fractures in bonded joints, which are based on the fracture toughness criterion. In this work, single lap joints (SLJ) – as typical mixed-mode loading coupon tests – with different bond-line thickness (0.20 and 0.90 mm) were tested to verify the ability of the proposed FE methods. Four different FE models based on fracture toughness criterion, i.e., cohesive elements, surface-based cohesive, extended finite element methods (XFEM) and virtual crack closure technique (VCCT) were implemented for predicting the failure behaviour of SLJs. This work aims to provide an overview of the applicability of the above-mentioned FE methods (which are based on fracture toughness criterion) for the prediction of fracture behaviour in adhesive SLJ. It was proved that the FE techniques used in combination with correct selection of parameters for constitutive damage law can lead to accurate prediction of the fracture load (Pm). *

Corresponding author. E-mail address: [email protected] (M.Z.Sadeghi).

1

Keywords: Finite element stress analysis, Cohesive zone model, Fracture toughness, Joint design Single lap joint

1. Introduction:

Over the last 70 years, adhesive bonded joints have been increasingly used in structural applications and are observed as alternatives to traditional methods of mechanically fastened and welded joints. Amongst various industries, aerospace has contributed immensely toward development of adhesives with the increase in the use of hybrid structures in new aircraft. This is due to the many advantages they offer compared with other joining techniques, i.e., flexibility in the assembly process, weight reduction, decrease in the stress concentration in the joint, etc. [1]. Structural adhesives can be bonded in different forms such as single and double lap joints, scarf joints, single L and T joints, etc. Adhesively bonded joints are frequently subjected to either static or fatigue loadings. Hence, attempting a conservative design to compensate for the negative effects of such loadings should be considered by designers. Failure prediction of adhesively bonded joints, which is generally divided into different categories, namely, analytical and numerical methods, has also received attention by researchers over the past years. One of the early works about the analytics of bonded joints can be attributed to Volkersen [2], which only considered the effect of shear stress in the adhesive layer in single lap joints. This model was further modified by Goland and Reisnner [3] by accounting for bending behaviour caused due to load eccentricity. Their model led to the determination of peel stress adjacent to the free edges of the adhesive layer. The literature has documented closed-form solutions offered for the adhesively bonded joints [4]. Although various analytical models were developed over the last decades, they have their inherent assumptions (for example, designed for a special type of 2

bonded joints). Hence, analytical models are insufficient to model complex material nonlinearity, different boundary conditions, and complex geometries. On the contrary, numerical methods (such as finite element [FE]) have shown to be promising alternatives to overcome such difficulties. FE models offered for analysis of the adhesively bonded joints cover a wide variety of different topics [5, 6]: Geometrical effects such as stress singularity field at the free edge of the bonded joints [7, 8], effect of the bond-line thickness [9], loading conditions such as fatigue [10-13], environmental effects [14, 15], etc. FE methods used for the analysis interfacial failure can be generally divided into different categories [6]: continuum mechanics, fracture mechanics and the more recent damage mechanics. First attempts on FE analysis of adhesively bonded joints have been reported in early 1970s [16, 17]. However, they are incapable of addressing the issue of stress singularity at the free edge of the overlapped region of adhesive bonded joints, which lead to convergence problems [18]. Furthermore, continuum models assume the material is continuum without any inherent flaws. Methods based on fracture mechanics also suffer limitations such as restriction on linear elastic fracture mechanics (LEFM) and reliance on the initial crack length such as the virtual crack closure technique (VCCT) [19-21] . Various other approaches, which have been developed, do not seem to be promising, as they induce additional empirical parameters [22, 23]. Among existing FE damage models available for modelling fractures in interfacial failures (such as composite delamination, bonded joints), cohesive zone models (CZM) are receiving increasing attention. CZM is able to predict the onset of damage initiation and propagation without defining pre-crack. This has made CZM preferable compared with other methods such as continuum mechanics and fracture mechanics. CZM is nowadays implemented in most of FE solvers (such as Abaqus, Ansys, Ls-Dyna, etc.) The most representative form of CZM in FE solvers is the cohesive element. However, this kind of element can be implemented in other 3

modelling techniques such as cohesive contact in the interface (surface-based cohesive contact) or as a damage model for XFEM. Constitutive response for cohesive elements follows tractionseparation law (TSL), which is capable of simulating gradual degradation of the materials based

on the simple correlation between the traction (T) and relative displacement (separation
). TSL can have different shapes, e.g., trapezoidal, exponential, etc.; however, the simplest approach is a bilinear softening law [24-26]. That being said, experimental work is still required to determine reliable input parameters for the traction-separation law (TSL) to achieve good fit with experimental results [24, 27]. Therefore, reliable determination of the critical strain energy release rate (CSEER) – as the most important fracture characteristic parameters in TSL – is essential when the CZM is used [24, 28, 29]. Cohesive elements have been successfully utilised by researchers in failure prediction of adhesively bonded joints [24, 28, 29]. In the cohesive element approach, the cohesive connection in the interface layer is modelled by using one raw cohesive element. On the contrary, in the surface-based cohesive method, the interfacial connection is modelled with a zero-thickness layer representing the contact in the interface. This approach is also successfully implemented for prediction of failure behaviour of adhesively bonded joint [30]. A major breakthrough made by Belytschko and his colleagues for introducing the fundamental idea of partition of unity led to the eXtended finite element method (XFEM) [31]. XFEM provides the ability to propagate the crack despite the discontinuity induced by the existing crack, which reduces massively the effort required for remeshing in conventional FE modelling. For modelling onset of damage initiation as well as propagation, XFEM is used with constitutive damage models. In Abaqus, for instance, the constitutive damage model considered for XFEM is based on CZM. Many works have successfully implemented the XFEM method using CZM for failure analysis of adhesively bonded joints [32-36]. In this work, four different FE techniques available for modelling interface failures in Abaqus, i.e., cohesive element, surface-based cohesive, XFEM, and VCCT, are used for the modelling 4

fracture behaviour of single lap joints. The first three methods are directly based on CZM in which constitutive theory in damage initiation and propagation in the interface follows TSL. For the damage initiation, cohesive behaviour is based on ultimate stress (strain) in normal and transverse directions. For the damage evolution, the CZM follows the correlation existing between the CSEER under modes I and II fracture (such as Power Law or Benzeggagh–Kenane [B-K] criteria). In such techniques, damage initiation can be modelled without defining existing flaws in the model. However, the later method VCCT is based on the concept of fracture mechanics and requires an existing crack to be considered in the model before fracture analysis. Nevertheless, the constitutive response for the damage evolution is based on CSEER. Hence, the reason fracture toughness GC is a common parameter, which relates these techniques with each other. In this work, the applicability of these techniques in modelling failure behaviour of single lap joints (SLJs) is discussed and compared. The FE results were compared with experimental tests carried out for SLJs under different adhesive bond-line thickness.

2. Experimental Work:

In the present work, a SLJ was manufactured from a steel adherent (S235JRG with yield stress of 265 MPa) and a two-component epoxy adhesive (Araldite 2015- Huntsman Advanced Materials, Basle, Switzerland). The adhesive used is a ductile adhesive with large plastic flow before its failure [37]. The geometry of the SLJs was designed based on ISO4587 (Figure 1). Two different thickness for the adhesive were considered in the test campaign, i.e., 0.20 and 0.90 mm, to investigate the effect of second bending moment induced by the existing eccentricity due to the force pairs in such joints.

5

Figure 1-Geometry of the SLJ (dimensions in mm, width=25 mm)

The plate surfaces prior to the bonding assembly were cleaned with Acetone and Isopropanol and then grit blasted. Prior to bonding, a primer (SACO-SIL by DELO-Germany) was applied on the bonding surfaces. For each thickness, five specimens were produced. The desired bondline thickness was achieved using steel wires. After bonding, the specimens were left in the oven for 1 h at 800 C [38]. Tensile tests on SLJs were carried out by using an electrical Instron machine (Instron 5567 with maximum load capacity of 30 kN). Small metallic tabs were used for both ends of the specimen with the total thickness equal to the thickness of the bond-line were used to assure that SLJs under tensile testing have symmetric geometry. The tests were carried out under displacement control with the rate of 1.75 mm/min. A digital image correlation system (ARAMIS-GOM Germany) was used during the tests to record the local deformations (Figure 2).

6

Figure 2-SLJ installed on the Instron machine ready for the tensile test

3. Finite Element Models:

Numerical simulation was performed by a commercial FEM package Abaqus to predict the loaddisplacement curves of SLJ. Although better stress assessment can be obtained by considering a 3D model, a 2D model is generally preferred for parametric studies. The 2D FE models were developed using Python script, which massively reduced the effort of parametrical studies in this work (including investigation into the effect of adhesive thickness, different damage criterion, etc.). Because the FE models were developed in 2D models, the effect of element types, i.e., plain stress/strain elements, on the failure behaviour of the SLJs was investigated. Therefore, a four-node bilinear plane stress/strain quadrilateral, reduced integration (CPS4R/ CPE4R) with hourglass control was used both for the adherents and the adhesive, except in the cohesive element approach (explained later) in which the adhesive is modelled by using cohesive elements (COHD24).

7

Investigation into the mesh sensitivity was carried out to determine the most appropriate mesh size. To reduce the computation time, the mesh regions were divided into different mesh regions (Figure 3): Fine mesh region (along the overlap region to capture the stress gradients) and a coarse mesh region (area away from the overlap). However, it must be noted that, for cohesive elements, the overlap region modelled as a single row of elements (discussed later). The boundary conditions used in the FE models are shown in Figure 3. The joint is fully constrained at one end; on the other end, the specimen was constrained to only move in the axial direction (tensile loading on displacement control).

Figure 3: FE meshing and boundary conditions used for cohesive surface, VCCT and XFEM

3.1 Damage initiation and propagation criterion:

In particular, for typical fractures under pure mode I (double cantilever test [DCB]) or mode II (end notch failure [ENF]), damage propagation law is independent from the interactions between the fracture modes. However, considering the complex loading scenario of the SLJ, which is a mixed-mode loading (combination of mode I and II) interaction between the modes for the damage initiation and propagation should be considered. 8

As far as the CZM is concerned, at the onset of crack growth process, the constitutive response follows TSL , which is capable of simulating gradual degradation of the materials based on the simple correlation between the traction (T) and relative displacement (separation
). TSL can have different shapes, and the simplest approach is a bilinear softening law [26, 39]. A bilinear TSL was considered for the analysis, as it has been reported that FE techniques such as XFEM revealed accurate results [36]. For the bilinear TSL, there are two damage criteria, i.e., initiation and propagation. For the onset of damage initiation, the failure index of the SLJ was computed using both maximum nominal traction (MAXS) and quadratic traction (QUAD). QUADS Criterion t   t t    +   +   = 1 t t t

(1)

t   t t     ,    ,    = 1 t t t

(2)







MAXS Criterion

Where, tn is traction in pure normal direction, and ts and tt are tractions in shear directions, whereas tn0, ts0, and tt0 are their corresponding maximum values. As long as Eq. (1) or (2) is not met, the interface is in the elastic regime following its stiffness (Knn and Kss), and material does

not undergo damage. The Macaulay bracket   is added to ensure that the pure compression

state has no influence on damage initiation. Once either of these equations is met, the material begins degradation (i.e., damage evolution criteria). Degradation of the material stiffness (adhesive in this case) is characterised by the damage evolution criteria. Damage evolution 9

models are generally based on stress intensity factors (K) or strain energy release rate (G). However, under a mixed mode, energy criteria are preferred, as they often combine tension and shear [40]. Hence, a linearised power criterion was considered as shown below        =  +  +    

(3)

where GI, GII and GIII are the fracture energy values for mode I, II ,and III respectively [41], while their corresponding subscript values, which are denoted by C, influence the accuracy of experimental data from fracture tests under modes I, II, and III. In the present work, the CZM parameters for single joints are presented in Table 1. Table 1-material properties of Araldite 2015

Property

Value

Young modulus (N/mm2)

1850±210 a

Poison’s ratio

0.33b

Tensile failure strength (N/mm2)

21.63±1.61 a

Shear failure strength (N/mm2)

17.9±1.8 a

GIC (N/mm)

0.43±0.02 a

GIIC, GIIIC (N/mm)

4.7±0.34 a

Power Law exponent

1a

a

reference [39], b Manufacturer’s data

10

3.2 FE Methods:

In the present work, four FE techniques were utilised for modelling failure behaviour of the single lap joints, i.e., cohesive element method, cohesive surface method, VCCT, and XFEM. In the following, a description of each method is presented. 3.2.1

Cohesive element method

Initial contributions to cohesive zone models include the pioneering works by Dugdale [42] and Barenblatt [43] as early as 1960s. Among various methods proposed for simulation of the integrity of the interfaces and bonded joints over the last decades, cohesive elements are receiving increasing attention due to the many advantages they offer (such as predicting the state of damage initiation and evolution and the capability of implementing in the component without pre-crack). The material data needed for prediction the onset of damage initiation and propagation by using CZM for Araldite 2015 are shown in Table 1. In the model developed based on cohesive elements, the adhesive part was modelled by a single row of four-node 2D cohesive elements (COH2D4) that hold the two adherent surfaces together. Cohesive behaviour was assigned to the section of the adhesive layer, and MAXS and QUADS damage initiation criteria were incorporated. In cohesive elements, normal and tangential stiffness components (Enn and Ess) are interpreted as penalty stiffness , which is modified considering the thickness of the adhesive layer in the defined traction-separation section as well as in the elastic and shear modulus [44, 45]. As an adhesive layer of non-unit thickness was considered in the SLJ model, the corresponding stiffness value was updated, as mentioned in the Abaqus documentation [44]. This stiffness correction was done by scaling the nominal value by the non-unit thickness.

11

3.2.2

Surface-based cohesive method

The surface-based cohesive model is another approach used in the FE framework of the present study. In this approach, the interaction between the adhesive layer and the adherents are modelled by a cohesive connection with a zero-thickness layer, and its damage constitutive law is based on the TSL. The constitutive damage response is assigned in the property of contact interaction between the adherent and adhesive (contact to contact surface). Similar damage constitutive law (including the initiation and evolution) with the models using the cohesive element approach was used. Although this approach is similar to the cohesive element approach, a surface-based cohesive in comparison with cohesive elements is easier to be implemented and offers more abilities (such as ability to model to sticky surfaces, which can bond when they come into contact within the analysis) [44]. Adherents and adhesives are modelled as continuum elements by using plane stress/strain elements (CPS4R/CPE4R). The properties used for modelling the adhesive are given in Table 1.

3.2.3

Virtual Crack Closure Technique (VCCT)

VCCT is based on Irwin’s [46] crack closure integral and is utilised when plastic dissipation does not exist. Rybicki and Kanninen [20] proposed a modified contour integral that evaluates the work required for crack closure along an element length. This approach can be easily implemented in FE software, and the concerned energy release rates can be computed by using the VCCT method employed in Abaqus. VCCT has been used in various applications such as delamination of composite and failure of adhesively bonded joints [19, 47]. 12

VCCT is based on linear elastic fracture mechanics (LEFM) and thus relies on the existence of an initial crack or debonded region to identify crack tips [24]. Furthermore, there exists no requirement of specifying the damage initiation in this technique. Therefore, a hypothetical crack was introduced as damage was specified in SLJs by means of initial debonding (by not bonding the interfacial nodes between the adhesive layer and upper and lower adherents) at the interface layer (the layer between the adhesive layer and the adherents). For the other nodes in the interface layer, surface-to-surface contact was attributed. VCCT fracture criteria were considered as the property of the defined contact in the interface following Power Law mixed-mode behaviour. In fact, the nodes in the interface during the loading will debond when the fracture criterion, f =

 !  !

(4)

reaches 1.0 with the following tolerance 1 ≤  ≤ 1 + "#$

(5)

in the above relations, GequivC is the value of GIC and GIIC given in Table 1 and Gequiv is the equivalent strain energy release rate induced in the interface while loading. As for f, value of 0.01 (1%) was chosen for our analysis and a suitable desired convergence behaviour was obtained. VCCT is similar to CZM in some aspects: They are both capable of modelling interfacial shearing and propagation of delamination cracks. Also, they use an elastic damage constitutive law for modelling fracture behaviour in the material once the crack is initiated. However, the main difference between VCCT and CZM is that, in VCCT, the onset of damage initiation and propagation is based on fracture energy criterion. In CZM, however, the fracture criterion is only applicable for damage propagation. 13

3.2.4

XFEM

Modelling stationary and propagating cracks by conventional FE methods requires considerable effort for remeshing, as the mesh needs to conform the discontinuities induced by the existence of the crack. XFEM by special enriched features, together with an additional degree of freedoms (DOF), can follow the discontinuities induced by the existence of cracks. The idea of XFEM was initially introduced by Belytschko and Black which is based on the concept of partition of unity [31]. Therefore, XFEM is the combination of conventional FE and enriched functions. Enrichment functions consist of asymptotic crack tip functions (for arresting singularities at the crack tip) and Heaviside functions (for jump displacement happening in the crack faces). An extended finite element displacement field can be expressed as [44]. 1 & u (X) = * N, (X) 0u, + 0 ,∈. /

H(X)a 45657,

89:;< <=9 9>

+

E

L

* FC (X)bC, K 45 565557 K CFG5 H>:?I J

(6)

In which NI and uI are the FE conventional shape functions and nodal DOF, respectively; in addition, the Heaviside enrichment term, which is for the nodes whose elements are cut by crack, consists of H(X), and aI is nodal enriched DOF (representing jump discontinuity). The crack tip enrichment part contributes with the nodes in the elements having crack tip. Fα(X)

is the asymptotic crack tip function, and bC, is the nodal DOF for the corresponding domain. By

this arrangement, the FE foundation can be still fulfilled (sparsity and symmetry of the matrices are reserved). XFEM domain was defined in the adhesive layer which is the region fracture happens. Within this domain, crack initiation criteria can be specified based on QUADS and MAXS criteria 14

Eqs.(1) and (2) respectively. Damage propagation was specified by Power Law (Eq. (3)) (properties listed in Table 1). Furthermore, the maximum number of attempts (in the analysis step control) were increased before abandoning the increment to 20 from 5 (default value) in the analysis step controls as per Abaqus documentation to improve the convergence behaviour of the analysis [44]. For the model using XFEM approach, adherents and the adhesive were modelled with plane stress/strain quadrilateral elements. Non-linear geometric effects were considered in all the developed FE models in order to consider large displacements happening during the loading.

4. Results and discussion 4.1 Experimental:

The output of the tensile tests on the SLJs from the Instron machine was load-displacement curves. Each group of the tests showed a good consistency in terms of load-displacement curves. The average values achieved for the maximum load (Pm) for the thickness of 0.20 mm and 0.90 mm were 8666.74±198.26 N and 7100.56±121.84 N respectively. The representative load-displacement curves for the effect of bond-line thickness is shown in the Figure 4. It can be inferred that the increase in the bond-line thickness decreases the joint strength. This was reported by many researchers [9, 48, 49] and it was attributed to increase in second bending moment induced by the existing eccentricity due to the force pairs in such joints. Symmetry fracture pattern was observed for the adhesive on both adherents (Figure 5). Local deformations during the tensile tests were recorded by ARAMIS system. The deformation of one of the SLJs with the bond-line thickness of 0.90 mm, some moment prior and after the

15

fracture is shown in Figure 6. Final deformation of SLJ with 0.20 mm bond-line thickness at the end of fracture is shown in Figure 7.

Figure 4: Representative experimental load-displacement curves for SLJ with bond-line thickness 0.20 and 0.90 mm

Figure 5: Fracture surface with adhesive bond-line thickness of 0.20 mm (left) and 0.90 mm (right)

Figure 6: deformation of SLJ with 0.90 mm bond-line thickness prior to the failure (left) and after failure (right)

16

Figure 7: SLJ with 0.20 mm bond-line thickness after failure taken by ARAMIS system

4.2 FEM methods:

A mesh study was performed and it was found that influence of mesh refinement (in particular for the mesh sizes less than 0.10 mm) on the increase of Pm was negligible for the frameworks considered in this study. Hence, a mesh size of 0.10 mm was adapted to sufficiently to predict the Pm with reasonable accuracy to experimental results. Furthermore, since the developed FE models are two dimensional, the effect of plane stress and plane strain element were also investigated on the load-displacement curves. SLJ models using plain strain elements revealed a higher stiffness value (from loaddisplacement curves) in comparison to a plane stress elements. In most cases, plane stress under predicted the Pm which is in accordance with the previous study [50]. For the majority of the results, the FE models were not able to capture the response of the load-displacement in post failure state completely like the experimental results. This might be attributed to the softening behaviour used in the present work for the adhesive in damage propagation law. Due to ductile nature of the Araldite 2015, the bi-linear TSL used in this work may not be suitable to simulate the post damage propagation [39]. Cohesive element model: Load-displacement curves for bond-line thickness of 0.20 mm and 0.90 mm are respectively as shown in Figure 8 and Figure 9. 17

Figure 8: Force-displacement curves for SLJ modelled by cohesive elements - 0.20 mm: plane strain (left) and plane stress (right).

Figure 9: Force-displacement curves for SLJ modelled by cohesive elements -0.90 mm: plane strain (left) and plane stress (right).

In comparison to MAXS criterion, the QUAD criterion results were conservative. Furthermore, for a bond-line thickness (ta = 0.20 mm), both plane stress and plane strain yield similar results of the Pm. However, with increase in bond-line thickness (ta = 0.90 mm), both methods yield different Pm.In such a case; plane strain elements provide a more reasonable prediction of fracture load. Comparison of Pm for plane strain and plane stress elements is shown in Figure 10.

18

Figure 10: Fracture load comparison FE vs. experimental for different bond-line thickness: ta = 0.20 mm (left) & ta = 0.90 mm (right)

As an output variable, overall scalar stiffness degradation (SDEG) can be requested in Abaqus. SDEG shows the incremental change of damage and its value can change from zero (representing no damage) and unity (representing fully damaged element). Elements with the value of SDEG equal to unity are removed representing the crack propagation. At the time step in which Pm achieved, element deletion is observed. First element deletion represents the state of crack initiation and followed by crack propagation (deletion of other elements in the path defined for cohesive elements). The deformed SLJ for the bond-line thickness (ta = 0.90 mm) is as shown in Figure 11.

Figure 11: SDEG parameter plot for SLJ with ta: 0.90 mm

The main limitation of cohesive elements is that, unlike XFEM, they cannot model the real fracture behaviour in SLJs which is called cohesive failure (fracture in the adhesive layer in a thin layer close to the interface). 19

Surface-based cohesive Model: Force-displacement curves for bond-line thickness of 0.20 mm and 0.90 mm are respectively as shown in Figure 12 and Figure 13.

Figure 12: Force-displacement curves for SLJ modelled by surface-based cohesive -0.20 mm: plane Strain (left) & plane Stress (right).

Figure 13: Force-displacement curves for SLJ modelled by cohesive surface-0.90 mm: plane Strain (left) & plane Stress (right).

In comparison to MAXS criterion, the QUAD criterion results were conservative. Similar to cohesive elements technique, surface-based cohesive techniques also yields to similar results of Pm for a bond-line thickness of ta = 0.20 mm for both plane stress and plane strain elements.

20

However, with increase in the bond-line thickness (ta = 0.90 mm), both methods yield to different Pm. In such a case; plane strain elements provide a more reasonable prediction of fracture load. Comparison of Pm for plane strain and plane stress elements is shown in Figure 14.

Figure 14: Fracture load comparison FE vs. experimental for different bond-line thickness: ta = 0.20 mm (left) & ta = 0.90 mm (right)

The final deformation of the SLJ model for the bond-line thickness (ta = 0.90 mm) is as shown in Figure 15. Similar to cohesive elements, this technique is also unable to simulate the real cohesive failure.

Figure 15: Deformation of SLJ for ta = 0.90 mm predicted by surface-based cohesive model

21

VCCT: As explained earlier, different initial crack lengths were considered for the FE models developed for SLJs based on VCCT to check the influence of initial crack length on the predicted fracture behaviour. Force-displacement curves for bond-line thickness of 0.20 mm and 0.90 mm for various initial bond lengths are as shown in Figure 16.

Figure 16: Force-displacement curves for SLJ modelled by VCCT for a bond-line thickness -0.20 mm (left) and 0.90 mm (right)effect of different initial crack length

From above, it is evident that an initial debond length of 4.0 mm gives a suitable approximation to the joint of 0.20 mm adhesive thickness. Similarly, for 0.90 mm adhesive thickness, an initial debond of 0.50 mm gives a similar response to the experimental curve. Interestingly enough, the most appropriate initial debond length was achieved by considering the case where it most aptly represents the experimental fracture load at its respective critical displacement. However, the reason for choosing such initial crack lengths is unclear and is a subject of future work. The comparison table considering the variations in Pm is shown in Figure 17.

22

Figure 17: Fracture load comparison FE vs. experimental for different bond-line thickness: ta = 0.20 mm (left) & ta = 0.90 mm (right)

As it was previously discussed, damage criterion for crack initiation is not needed to be implemented in this technique. The deformed SLJ with the bond-line thickness (ta = 0.90 mm which is similar with the surface-based cohesive modelling technique is shown in Figure 18.

Figure 18: Deformation of SLJ with ta = 0.90 mm predicted by VCCT model

In case of VCCT, specification of initial debond length can be a subject of future study where there is sufficient scope to develop effective hypothetical crack length required to study damage initiation.

XFEM:

Force-displacement curves (Pm -δ) for bond-line thickness of 0.20 mm and 0.90 mm are respectively as shown in Figure 19 and Figure 20.

23

Figure 19: Force-displacement curves for SLJ modelled by XFEM-0.20 mm: plane Strain (left) and plane Stress (right).

Figure 20: Force-displacement curves for SLJ modelled by XFEM-0.90 mm: plane Strain (left) and plane Stress (right).

In comparison to MAXS criterion, the QUAD criterion results were conservative (except for plane stress elements bond-line thickness (ta = 0.20 mm). Similar to cohesive based methods (elements & surface), XFEM technique also yields to similar results of Pm for a bond-line thickness (ta = 0.20 mm) for both plane stress and plane strain elements. However, with increase in bond-line thickness (ta = 0.90 mm), both methods yield different Pm. In such a case; plane strain elements provide a more reasonable prediction of fracture load. Comparison of Pm for plane strain and plane stress elements is shown in Figure 21

24

Figure 21: Fracture load comparison FE vs. experimental for different bond-line thickness: ta = 0.20 mm (left) & ta = 0.90 mm (right)

PHILSM was requested was requested as an output variable. It enables Abaqus to display the location of cracked elements belonging to an enriched region. In fact, along an enriched region, these parameters are non-zero, visualising the cracked area. The status of the signed distance function (PHILSM) when the Pm is reached (bond-line thickness of 0.90 mm) is as shown in Figure 22. Furthermore, the values of PHILSM can be interpreted as the nodal coordinates of the enriched nodes with the coordinate system centered at the crack front. If the crack front is along the element edges, the PHILSM value is set to zero, as this would mean that the particular element is not enriched.

Figure 22: PHILSM parameter plot for SLJ with ta: 0.90 mm

25

In XFEM, the crack starts in the first elements at the both ends close to the adherents due to the stress concentration at the both ends of the edges. This is from where crack initiation originates propagating from both ends towards each other. However, in comparison to cohesive based methods, XFEM methods have some inherent limitations. To begin with, the most critical limitation was the crack propagation towards the adherents (not defined in the XFEM domain) for small bond-line thickness (in our case ta = 0.20 mm). Subsequently the load continues to rise till the crack tip encounters the adherents. Right before the point where the crack tip encounters the adherent, Pm is reached and in reality, the bond has lost its load bearing capacity. However, XFEM method fails to predict this behaviour and the load continues to rise even after the crack tip has encountered the adherents. This has been observed in previous studies as well [51]. In contrast, like cohesive based techniques, XFEM predicts well for the bond-line thickness of 0.90 mm. In this case, the cracks would always propagate towards each other meeting in the middle of the bond-line as expected. Another limitation is the mesh dependency .on crack initiation as stated in the previous study [32] However; this was addressed in the current SLJ by considering a fine mesh size of 0.10mm. Finally, XFEM domain has to be defined This is in accordance with the previous finding [52] in which the XFEM method was identified as a method which is suitable for recognising the locus of damage initiation rather than predicting the interfacial failure of the bonded joints.

5. Conclusions:

In this work, the applicability of four different FE techniques (2D) for modelling failure in the adhesively bonded single lap joints was investigated. Three techniques (i.e., cohesive element, surface-based cohesive, XFEM) are directly based on a cohesive zone model and the damage 26

initiation and evolution follows traction-separation law. For the damage evolution, a Power Law criterion was used to define the correlation of the strain energy release rate of the joint under loading. Additionally, VCCT was also used in the FE modelling in which the damage (as existing crack) should be predefined in the model. For the damage evolution constitutive law, Power Law criteria were used. By using Python script, different parametrical studies were developed to investigate the effect of element type, damage initiation criteria, etc. The following conclusions can be drawn: 1- All the FE techniques showed good capability for the prediction of the failure load Pm for different adhesive bond-line thickness. 2- As expected, with increase in bond-line thickness from 0.20 mm to 0.90 mm, the Pm decreased. This can be attributed to the increase in secondary bending moments induced in the joints. The FE models were able to predict successfully this observation. 3- In the developed 2D models, plain strain elements would always lead to higher stiffness in the load displacement curve compared with plain stress elements. 4- In VCCT, no damage initiation criterion needs to be defined. It showed good capability for predicting the failure behaviour of the joints; however, compared with the other techniques, it is less practical, as the pre-crack should be defined in the model, and its length should be defined depending on the bond-line geometry. 5- In all models presented in this study, for a smaller bond-line thickness (ta = 0.20 mm); plane stress and plane strain elements revealed similar Pm. However, plane strain elements lead to more representative results of Pm with increase in bond-line thickness (ta = 0.90 mm). 6- Quadratic damage initiation criterion would always lead to more conservative fracture load except for the plane stress XFEM technique.

27

7- XFEM technique proved inaccurate for modelling SLJ with small bond-line thickness (ta = 0.20 mm). However, for higher bond-line thickness, it can be employed to like other cohesive zone model methods. 8- Except XFEM method, the crack is always initiated at the bond-line for cohesive based methods and VCCT. Hence these methods inaccurately model the locus of crack initiation. 9- An improved accuracy for all the methods could be reached by considering a modification from bi-linear to trapezoidal TSL to improve the damage propagation characteristics in ductile adhesives. Though such techniques showed high potential in modelling failure of adhesively bonded joints, it should be mentioned that the damage evolution criteria is highly dependent on the reliable correlation between strain energy release rate under different mixed-mode loading. In the present work, GIC , GIIC and the damage evolution (Power Law with an index of 1.0) was used for Araldite 2015 were taken from literature which have been successfully in FE framework based on CZM. However, the applicability of different damage evolution laws combined with the FE frameworks discussed (for example Power Law with higher orders) under full envelope of fracture attained by different mixed-mode ratios still needs to be investigated.

28

References:

[1] Da Silva L.F.M., Öchsner A., Adams R.D., Handbook of adhesion technology-2nd edition. Heidelberg: Springer; 2018. [2] Volkersen O., The rivet-force distribution in tension-stressed riveted joints with constant sheet thicknesses, Luftfahrtforschung 1938; 1: 15-41. [3] Reissner E., Goland M., The stress in cemented joints, J Appl Mech 1944; 66: A17-A27. [4] Da Silva L.F., das Neves P.J., Adams R., Spelt J., Analytical models of adhesively bonded joints—Part I: Literature survey, Int J Adhes Adhes 2009; 29: 319-30 https://doi.org/10.1016/j.ijadhadh.2008.06.005. [5] He X., A review of finite element analysis of adhesively bonded joints, Int J Adhes Adhes 2011; 31: 248-64 http://dx.doi.org/10.1016/j.ijadhadh.2011.01.006. [6] Da Silva L.F., Campilho R.D., Advances in numerical modelling of adhesive joints. Heidelberg: Springer; 2012. [7] Van Tooren M., Gleich D., Beukers A., Experimental verification of a stress singularity model to predict the effect of bondline thickness on joint strength, J Adhes Sci Tech 2004; 18: 395-412 https://doi.org/10.1163/156856104323016315. [8] Goglio L., Rossetto M., Stress intensity factor in bonded joints: influence of the geometry, Int J Adhes Adhes 2010; 30: 313-21 https://doi.org/10.1016/j.ijadhadh.2009.12.008. [9] Akhavan-Safar A., Ayatollahi M.R., da Silva L.F.M., Strength prediction of adhesively bonded single lap joints with different bondline thicknesses: A critical longitudinal strain approach, Int J Solid Struct 2017; 109: 189-98 https://doi.org/10.1016/j.ijsolstr.2017.01.022. [10] Khoramishad H., Crocombe A.D., Katnam K.B., Ashcroft I.A., Predicting fatigue damage in adhesively bonded joints using a cohesive zone model, Int J Fatigue 2010; 32: 1146-58 http://dx.doi.org/10.1016/j.ijfatigue.2009.12.013. [11] Imanaka M., Hamano T., Morimoto A., Ashino R., Kimoto M., Fatigue damage evaluation of adhesively bonded butt joints with a rubber-modified epoxy adhesive, J Adhes Sci Tech 2003; 17: 98194 https://doi.org/10.1163/156856103322112888. [12] Monteiro J., Akhavan-Safar A., Carbas R., Marques E., Goyal R., El-zein M., et al., Mode II modeling of adhesive materials degraded by fatigue loading using cohesive zone elements, Theor Appl Fract Mech 2019; 103: 102253 https://doi.org/10.1016/j.tafmec.2019.102253. [13] Graner Solana A., Crocombe A.D., Ashcroft I.A., Fatigue life and backface strain predictions in adhesively bonded joints, Int J Adhes Adhes 2010; 30: 36-42 https://doi.org/10.1016/j.ijadhadh.2009.08.001. [14] Crocombe A.D., Hua Y.X., Loh W.K., Wahab M.A., Ashcroft I.A., Predicting the residual strength for environmentally degraded adhesive lap joints, Int J Adhes Adhes 2006; 26: 325-36 https://doi.org/10.1016/j.ijadhadh.2005.04.003. [15] da Silva L.F.M., Adams R.D., Adhesive joints at high and low temperatures using similar and dissimilar adherends and dual adhesives, Int J Adhes Adhes 2007; 27: 216-26 https://doi.org/10.1016/j.ijadhadh.2006.04.002. [16] Carver D., Wooley G., Stress concentration factors for bonded lap joints, J Aircraft 1971; 8: 817-20 https://doi.org/10.2514/3.44305. [17] Adams R., Peppiatt N., Stress analysis of adhesive-bonded lap joints, J Strain Anal 1974; 9: 185-96 https://doi.org/10.1243/03093247V093185. [18] Towse A., Potter K., Wisnom M., Adams R., The sensitivity of a Weibull failure criterion to singularity strength and local geometry variations, Int J Adhes Adhes 1999; 19: 71-82 https://doi.org/10.1016/S0143-7496(98)00058-X. [19] Jokinen J., Wallin M., Saarela O., Applicability of VCCT in mode I loading of yielding adhesively bonded joints—a case study, Int J Adhes Adhes 2015; 62: 85-91 https://doi.org/10.1016/j.ijadhadh.2015.07.004. [20] Rybicki E.F., Kanninen M.F., A finite element calculation of stress intensity factors by a modified crack closure integral, Eng Fract Mech 1977; 9: 931-8 https://doi.org/10.1016/0013-7944(77)90013-3. 29

[21] Eder M.A., Bitsche R., Fracture analysis of adhesive joints in wind turbine blades, Wind Energy 2015; 18: 1007-22 https://doi.org/10.1002/we.1744. [22] Abdel Wahab M.M., Hilmy I., Ashcroft I.A., Crocombe A.D., Evaluation of Fatigue Damage in Adhesive Bonding: Part 1: Bulk Adhesive, J Adhes Sci Tech 2010; 24: 325-45 https://doi.org/10.1163/016942409X12508517390798. [23] Abdel Wahab M.M., Hilmy I., Ashcroft I.A., Crocombe A.D., Evaluation of Fatigue Damage in Adhesive Bonding: Part 2: Single Lap Joint, J Adhes Sci Tech 2010; 24: 325-45 https://doi.org/10.1163/016942409X12508517390833. [24] Campilho R.D.S.G., Banea M.D., Neto J.A.B.P., da Silva L.F.M., Modelling of Single-Lap Joints Using Cohesive Zone Models: Effect of the Cohesive Parameters on the Output of the Simulations, Int J Adhes Adhes 2012; 88: 513-33 https://doi.org/10.1080/00218464.2012.660834. [25] Fernandes R., Campilho R., Testing different cohesive law shapes to predict damage growth in bonded joints loaded in pure tension, J Adhes 2017; 93: 57-76 https://doi.org/10.1080/00218464.2016.1169176. [26] Alfano G., On the influence of the shape of the interface law on the application of cohesive-zone models, Compos Sci Technol 2006; 66: 723-30 http://dx.doi.org/10.1016/j.compscitech.2004.12.024. [27] Ferreira L.R.F., Campilho R.D.S.G., Rocha R.J.B., Barbosa D.R., Geometrical and material optimization of tensile loaded tubular adhesive joints using cohesive zone modelling, J Adhes 2019; 95: 425-49 https://doi.org/10.1080/00218464.2018.1551136 [28] Abrate S., Ferrero J.-F., Navarro P., Cohesive zone models and impact damage predictions for composite structures, Meccanica 2015; 50: 2587-620 https://doi.org/10.1007/s11012-015-0221-1. [29] Xu Y., Li X., Wang X., Liang L., Inverse parameter identification of cohesive zone model for simulating mixed-mode crack propagation, Int J Solid Struct 2014; 51: 2400-10 https://doi.org/10.1016/j.ijsolstr.2014.03.008. [30] Sadeghi M.Z., Zimmermann J., Gabener A., Schroeder K.U., The applicability of J-integral approach in the determination of mixed-mode fracture energy in a ductile adhesive, Int J Adhes Adhes 2018; 83: 2-8 https://doi.org/10.1016/j.ijadhadh.2018.02.027. [31] Belytschko T., Black T., Elastic crack growth in finite elements with minimal remeshing, Int J Numer Meth Eng 1999; 45: 601-20 https://doi.org/10.1002/(SICI)1097-0207(19990620)45:5<601::AIDNME598>3.0.CO;2-S. [32] Campilho R.D.S.G., Banea M.D., Chaves F.J.P., Silva L.F.M.d., eXtended Finite Element Method for fracture characterization of adhesive joints in pure mode I, Comp Mater Sci 2011; 50: 1543-9 http://dx.doi.org/10.1016/j.commatsci.2010.12.012. [33] Campilho R.D.S.G., Banea M.D., Pinto A.M.G., da Silva L.F.M., de Jesus A.M.P., Strength prediction of single- and double-lap joints by standard and extended finite element modelling, Int J Adhes Adhes 2011; 31: 363-72 http://dx.doi.org/10.1016/j.ijadhadh.2010.09.008. [34] Mubashar A., Ashcroft I.A., Crocombe A.D., Modelling Damage and Failure in Adhesive Joints Using A Combined XFEM-Cohesive Element Methodology, J Adhes 2013; 90: 682-97 https://doi.org/10.1080/00218464.2013.826580. [35] Bendemra H., Compston P., Crothers P.J., Optimisation study of tapered scarf and stepped-lap joints in composite repair patches, Compos Struct 2015; 130: 1-8 http://dx.doi.org/10.1016/j.compstruct.2015.04.016. [36] Xará J., Campilho R., Strength estimation of hybrid single-L bonded joints by the eXtended Finite Element Method, Compos Struct 2018; 183: 397-406 https://doi.org/10.1016/j.compstruct.2017.04.009. [37] Sadeghi M.Z., Zimmermann J., Gabener A., Schröder K.U., On the bondline thickness effect in mixed-mode fracture characterisation of an epoxy adhesive, J Adhes 2019: 1-23 https://doi.org/10.1080/00218464.2019.1656614. [38] A. G., Methodik zur Auslegung von Klebverbindungen mit variabler Klebschichtdick [PhD], RWTH Aachen University, 2016.

30

[39] Campilho R.D.S.G., Banea M.D., Neto J.A.B.P., da Silva L.F.M., Modelling adhesive joints with cohesive zone models: effect of the cohesive law shape of the adhesive layer, Int J Adhes Adhes 2013; 44: 48-56 http://dx.doi.org/10.1016/j.ijadhadh.2013.02.006. [40] Ridha M., Tan V., Tay T., Traction–separation laws for progressive failure of bonded scarf repair of composite panel, Compos Struct 2011; 93: 1239-45 https://doi.org/10.1016/j.compstruct.2010.10.015. [41] Chaves F.J.P., da Silva L.F.M., de Moura M.F.S.F., Dillard D.A., Esteves V.H.C., Fracture Mechanics Tests in Adhesively Bonded Joints: A Literature Review, J Adhes 2014; 90: 955-92 https://doi.org/10.1080/00218464.2013.859075. [42] Dugdale D.S., Yielding of steel sheets containing slits, J Mech Phys Solids 1960; 8: 100-4 https://doi.org/10.1016/0022-5096(60)90013-2. [43] Barenblatt G.I., The Mathematical Theory of Equilibrium Cracks in Brittle Fracture, Adv Appl Mech 1962; 7: 55-129 https://doi.org/10.1016/S0065-2156(08)70121-2. [44] ABAQUS 6.13-Software Package Dassault Systems 2013. [45] Castagnetti D., Dragoni E., Standard finite element techniques for efficient stress analysis of adhesive joints, Int J Adhes Adhes 2009; 29: 125-35 https://doi.org/10.1016/j.ijadhadh.2008.01.005. [46] Irwin G.R., Analysis of stresses and strains near the end of a crack traversing a plate, J. appl. Mech. 1957. [47] Shokrieh M.M., Rajabpour-Shirazi H., Heidari-Rarani M., Haghpanahi M., Simulation of mode I delamination propagation in multidirectional composites with R-curve effects using VCCT method, Comp Mater Sci 2012; 65: 66-73 https://doi.org/10.1016/j.commatsci.2012.06.025. [48] da Silva L.F.M., Rodrigues T.N.S.S., Figueiredo M.A.V., de Moura M.F.S.F., Chousal J.A.G., Effect of Adhesive Type and Thickness on the Lap Shear Strength, J Adhes 2006; 82: 1091-115 https://doi.org/10.1080/00218460600948511. [49] Banea M.D., da Silva L.F.M., Campilho R.D.S.G., The Effect of Adhesive Thickness on the Mechanical Behavior of a Structural Polyurethane Adhesive, J Adhes 2014; 91: 331-46 10.1080/00218464.2014.903802. [50] Rocha R.J.B., Campilho R.D.S.G., Evaluation of different modelling conditions in the cohesive zone analysis of single-lap bonded joints, J Adhes 2018; 94: 562-82 https://doi.org/10.1080/00218464.2017.1307107. [51] Sousa C.d., Campilho R., Marques E., Costa M., Silva L.d., Overview of different strength prediction techniques for single-lap bonded joints, Proc Inst Mech Eng Pt. L J Mater Des Appl 2017; 231: 210-23 https://doi.org/10.1177/1464420716675746. [52] Fernandes T.A., Campilho R.D., Banea M.D., da Silva L.F., Adhesive selection for single lap bonded joints: experimentation and advanced techniques for strength prediction, J Adhes 2015; 91: 841-62 https://doi.org/10.1080/00218464.2014.994703.

31