A review on failure theories and simulation models for adhesive joints

ABSTRACT

In the framework of the Cost Action CERTBOND (Reliable roadmap for certification of bonded primary structures), a wide group of researchers from 27 European Countries have had the opportunity to work on the topic of certification of bonded joints for primary structural applications from different engineering sectors such as the aerospace, automotive, civil engineering, wind energy and marine sectors. Since virtual testing and optimization are basic tools in the certification process, one of the key objectives of CERTBOND is to critically review some of the available models and failure theories for adhesive joints. The present paper summarizes the outcome of this task. Nine different models/theories are described in detail. Specifically, reviewed are the Classical Analytical Methods, the Process Zone Methods, Linear Elastic Fracture Mechanics (LEFM), the Virtual Crack Closure Technique (VCCT), the Stress Singularity Approach, Finite Fracture Mechanics (FFM), the Cohesive Zone Method (CZM), the Progressive Damage Modeling method and the Probabilistic methods. Also, at the end of the paper, the modeling of temperature effects on adhesive joints have been addressed. For each model/theory, information on the methodology, the required input, the main results, the advantages and disadvantages and the applications are given.

KEYWORDS:

• Adhesive bonding
• analytical modeling
• numerical modeling
• debonding
• fracture mechanics
• cohesive zone modeling
• probabilistic modeling

1. Introduction

Adhesive joining technology finds an increasing use in assembly and repair applications in the aerospace, automotive, marine and civil structures sectors due to its many advantages over mechanical joining techniques .^[1–3] The adhesive joining applications are restricted to secondary structures due to stringent certification guidelines. The main reasons for the lack of certification of primary adhesively bonded structures are the uncertainty of the quality of the bondline, combined with the inability of the existing non-destructive testing methods to fully detect defects and characterize surface quality, and the inability to control by design the crack growth rate and length in the bondline. For the aerospace sector, the proposed means of compliance that describe the methodologies/measures that need to be undertaken are the following ones^[4]: (i) The maximum disbonds of each bonded joint consistent with the capability to withstand the loads […] must be determined by analysis, tests, or both. Disbonds of each bonded joint greater than this must be prevented by design features. (ii) Proof testing must be conducted on each production article that will apply the critical limit design load to each critical bonded joint; or (iii) Repeatable and reliable non-destructive inspection techniques must be established that ensure the strength of each joint. Evaluating these means of compliance in the context of large-scale structures, we conclude that the only feasible measure is the limitation of the maximum disbond size, whereas proof testing is cost-prohibitive and non-destructive techniques have not yet the capability to measure or correlate bonding strength to date. Given that the mean of compliance (i) may be reached by analysis and/or tests, the role of models towards certification is crucial.

In 2019, the COST Action CERTBOND (Reliable roadmap for certification of bonded primary structures)^[5] has started. The action aims to deliver a reliable roadmap for enabling certification of primary bonded structures from various industrial sectors. One of the specific tasks of CERTBOND is critically review and report the existing failure theories and simulation models for adhesive joints. The current paper is the outcome of this task.

In general, the modeling and failure simulation of adhesive joints is a difficult task due to the complicated phenomena that take place .^[6] These include the different failure mechanisms that extend from the adhesive to the adhesive/adherent interface and to the adherent, the effect of surface treatment and the presence of defects in the bondline, the local variation of stiffness and strength and the different materials behavior (elastic, plastic, brittle and ductile). The complexity becomes even higher if the adherents are made from composite material and if dynamic loads are applied. The requirements for the models are the computation of stress and strain fields in the adhesive, the simulation of crack initiation and growth in the bondline (debonding) and the simulation of damage in the adherents and especially the composite adherents. More advanced requirements are the modeling of defects in the bondline and the modeling of environmental effects.

In the literature, there have been published only a few reviews on the modeling of adhesive joints. Da Silva et al. have performed in a two-paper series^[7,8] an extensive literature review on existing analytical models for both single and double-lap joints and a comparative study. Ramalho et al.^[9] have reviewed static prediction methods for adhesive joints. They have included analytical and numerical methods and the largest part of the review is devoted to numerical methods. The review also contains a very useful critical discussion section. From the above short overview, it becomes evident that there is still plenty of room in the literature for critical reviews of the models of adhesive joints and especially for debonding simulation models and for models of complicated phenomena such as the behaviour under cyclic loadings and the temperature effects.

In the present paper, we review a series of failure theories and models for adhesive joints. The reviewed theories/models include Classical Analytical Methods, Process Zone Methods, Linear Elastic Fracture Mechanics (LEFM), the Virtual Crack Closure Technique (VCCT), Stress Singularity Approach, Finite Fracture Mechanics (FFM), the Cohesive Zone Method (CZM), the Progressive Damage Modeling method and Probabilistic methods. Also, at the end of the paper, the modeling of temperature effects on adhesive joints have been addressed. Some of the reviewed theories/models do not appear in any of the published reviews. The present review is the outcome of an interdisciplinary approach as the contributors work on different industrial sectors such as the aerospace, the automotive, civil engineering, wind turbines and marine sectors. Therefore, each method is approached through a different viewpoint (different objectives, requirements, specifications and applications), thus enabling the authors to present a more integrated evaluation of the selected failure theories and models. The authors realize that not all theories/models are included in the review, but this was not the scope of the work. The main topics covered in each section are the following: description of the theory/model, required parameters, main results, advantages-disadvantages, and applications (industrial sector, loads, etc.). At the end of the paper, a critical comparison between the different approaches is reported, highlighting their advantages, limitations and application fields.

2. Classical analytical methods

The quantitative evaluation of the strength of an adhesively bonded joint demands adequate knowledge of the stress distribution within the joint and a suitable failure criterion. The stress distribution can be obtained from the static or dynamic equilibrium equations for the joint, which might be solved analytically or by means of numerical methods such as the Finite Element Analysis (FEA). While FEA is the method of choice for complex geometries and nonlinear material models, analytical solutions are more appropriate to build up a general understanding of the problem in short time. Knowing that general-purpose finite element packages are too complex for everyday practical application a number of intermediate complexity calculation tools for widespread technical needs were developed. A short history of the development of such tools is provided in .^[10] The majority of stress (analytical) models in adhesive joints have been incorporated in the Joint Designer package .^[11] On the other hand, user-friendly JointCalc software^[10] provides accurate solutions for most of the joint configurations encountered in practice, while KISPEO software^[12] allows the calculation of stress distributions and failure loads for a SLJ configuration. Adhesive joints have been intensively investigated over the past 80 years and numerous analytical models have been proposed.

Most of the analytical models for adhesively bonded joints are two-dimensional. Disregarding the class of tubular joints loaded in torsion, in which the stresses run along the width of the joint, this statement can be acknowledged as sufficient because the stresses in the width direction are general significantly lower than in the direction of the loading. In these analyses, it is assumed that the adhesive joints are in a state of plane stress or plane strain in the plane perpendicular to the width direction, neglecting the stresses across the width direction caused by Poisson’s ratio strains in the adherends and the anticlastic bending of the adherends. However, there are some exceptions such as^[13] and^[14] that consider three-dimensional effects. Nonlinear material behaviour is difficult to incorporate because the analysis becomes very complex. Thus, most of the analyses are linear elastic for both adherends and adhesive.

As the degree of complexity and the number of stress components in the adhesive and the adherends increase, the initial analytical problem must be solved numerically. Table 1 in Volkersen’s seminal paper^[15] gives a detailed summary of the available analytical models, indicating the conditions of applicability the stresses they give, and the type of solution (algebraic or numerical). For example, if the joint bending is not severe and the adhesive is brittle, Volkersen’s analysis is sufficient. However, if the adhesive and/or the adherends undergo plastic deformation and substantial peeling is present, a more comprehensive model is necessary, which is typically less amenable to be solved via analytical methods.

Table 1. Summary of both linear and nonlinear two-dimensional analytical models available in the literature^[15].

Available analytical models	Material linearity				Adherends					Adhesive stresses
	Adhesive		Adherend		I	C	S	DS		${\sigma }_x$	${\sigma }_y$	${\tau }_{xy}$
	L	NL	L	NL				T	M
Volkersen^[Citation15]	X		X		X		X	X				X
Goland and Reissner^[Citation16]	X		X		X		X				X	X
Wah^[Citation17]	X		X		X	X	X	X				X
Hart-Smith^[Citation18]	X	X	X		X		X				X	X
Pirvics^[Citation19]	X		X		X		X	X	X	X	X	X
Grimes and Greimann^[Citation20]	X	X	X	X	X	X	X	X	X		X	X
Renton and Vinson^[Citation21,Citation22]	X		X		X	X	X	X	X		X	X
Srinivas^[Citation23]	X		X		X	X	X	X	X		X	X
Allman^[Citation24]	X		X		X	X	X				X	X
Ojalvo and Eidinoff^[Citation25]	X		X		X		X		X		X	X
Delale et al.^[Citation26]	X	X	X		X	X	X	X	X		X	X
Bigwood and Crocombe^[Citation27]	X		X		X		X	X	X		X	X
Bigwood and Crocombe^[Citation28]	X	X	X		X		X	X	X		X
Cheng et al.^[Citation29]	X		X		X		X	X	X	X	X	X
Crocombe and Bigwood^[Citation30]	X	X	X	X	X		X	X	X		X
Adams and Mallick^[Citation31]	X	X	X		X	X	X	X	X	X	X	X
Tong^[Citation32]	X	X	X		X		X				X	X
Yang and Pang^[Citation33]	X		X		X	X	X	X	X		X	X
Frostig et al.^[Citation34]	X		X		X	X	X	X	X		X	X
Sawa et al.^[Citation35]	X		X		X		X	X	X	X	X	X
Mortensen and Thomsen^[Citation36]	X	X	X		X	X	X	X	X			X
Adams et al.^[Citation37]	X	X	X	X	X		X
Wang et al.^[Citation38]	X	X	X	X	X		X	X	X	X	X	X
Smeltzer and Klang^[Citation39]	X	X	X	X	X	X	X	X	X		X	X

Notations: L – Linear, NL – Nonlinear I– Isotropic, C – Composite, S – Similar, DS – Dissimilar, T – Thickness, M – Material

The simplest linear elastic analysis considers one of the most common joints that can be found in practice, the single-lap joint (SLJ). In this analysis, the adhesive is assumed to deform only in shear, while the adherends stay rigid, as illustrated in Figure 1(a). The adhesive shear stress τ_xy is constant over the overlap length and is given by:

(1) ${ }_{xy}=\frac{F}{bL}$(1)

Figure 1. Deformations and shear stresses in loaded single-lap joint (SLJ). (a) Rigid adherends and (b) Elastic adherends.

where F is the applied load, b is the joint width and L is the overlap length. The value of the shear stress computed via equation 1(1) ${ }_{xy}=\frac{F}{bL}$(1) can be interpreted as the average shear stress acting on the adhesive layer. Although the modelling assumptions leading to equation 1(1) ${ }_{xy}=\frac{F}{bL}$(1) represent major simplifications of the actual system, they still provide the main guidance to evaluate the shear strength of adhesive joints in conformity to ASTM and ISO standards.

Volkersen’s analysis^[15] introduced the concept of differential shear, which is schematically illustrated in Figure 1(b). The assumption of pure shear deformation of the adhesive is maintained as in the simplest analysis, the adherends can deform in tension because they are considered elastic and not rigid. An effect that is not taken into account in this analysis is the bending moment due to the eccentric load path of SLJs. The solution is more representative of a Double Lap Joint (DLJs) than SLJs, since bending of the adherends in DLJs is not as significant as in SLJs.

The eccentric load path in a SLJ generates a bending moment and a transverse force at the ends of the joint, in addition to the applied tensile load per unit width. Because of this bending moment, the joint will rotate, altering the direction of the load line with the tendency of the applied tensile forces to come into line. As the joint rotates, the bending moment will decrease, giving rise to a nonlinear geometric problem where the effects of the large deflections of the adherends must be accounted for. The first to consider these effects were Goland and Reissner, ^[16] whose analysis was later extended by Oplinger, ^[40] who considered large deflections both outside and inside the overlap, as well as the individual deformation of the upper and lower adherends in the overlap. The results reported by Oplinger^[40] were close to those of Goland and Reissner^[16] in the case of large adherend-to-adhesive thickness ratio, whereas they differed significantly in the opposite limit of adherends thinner than the adhesive layer. Zhao^[41] developed a simpler form of the bending moment factor that is accurate for thick and stiff adherends but has limitations for short overlaps.

After the determination of the loads at the ends of the overlap, Goland and Reissner^[16] calculated the shear and peel stresses in the adhesive layer, solving a plane strain problem. Instead of solving a nonlinear geometric problem due to the eccentric load path, they solved a linear problem in the overlap with the loads applied at the ends. In this way, they avoided a more complex problem with the consideration of the geometric nonlinearity effect.

After the so-called classical works, some authors tried to obtain more general closed-form solutions by including, for example, adherends with dissimilar thickness and material properties or composite adherends .^{[9–23−25–53]} However, as the models get more general, their mathematical structure become increasingly complicated, thus making the need for numerical methods more prominent, either to evaluate stress values from closed-form solutions, ^[15,16] or to compute an approximate solution of the governing differential equations.

Models of adequate fidelity are particularly necessary when the joints include laminated composite adherends, since through-thickness (or transverse) shear and normal deformations in the adherends can be hardly neglected. The most important of the earlier analyses to account for these deformations were done by Renton and Vinson, ^[22] Srinivas^[23] and Allman .^[24] Due to the increase in use of composite materials at that time, Renton and Vinson^[21,22] suggested that the analysis should consider not only the anisotropic properties of composites, but also the laminated construction (anisotropic properties of each lamina and lamina fibre orientation). Using composite laminated plate theory, they developed a linear elastic analysis between two similar or dissimilar laminated or isotropic adherends for a SLJ.

Wah^[17] was the first to consider laminated composite adherends that were symmetrical about their midsurface. The adhesive shear stress was assumed to be constant through the thickness whereas the adhesive peel stress could vary. The case of asymmetric composite adherends in balanced or unbalanced joints was considered more recently by Yang and Pang^[33] and Mortensen and Thomsen^[36] where the coupling effect of the external tensile loading and the bending moment due to the asymmetry of the composite laminates was considered.

Most analytical methods indicate that the strength of adhesive joints is enhanced by thicker bondlines (i.e. the stress decreases as the bondline thickness increases) as reported, for instance, in the study by Srinivas .^[23] However, in practice, the adhesive lap joint strength decreases as the glueline gets thicker. Adams and Peppiatt^[44] proposed that this seemingly contradictory behaviour could be understood as an effect of the presence of defects, such as voids and micro cracks, which are presumably more numerous per unit volume in thicker bondlines. Other possible reasons for decreasing strength in thicker bondlines might be identified in the more pronounced adherend bending due to load misalignment, or interface adhesive-adherend stresses .^[45,46] Therefore, in order to realistically predict failure loads, analytical models should include the variation of stress through the adhesive thickness, including the interface stresses.

When adhesives having a large plastic strain to failure are used, such as rubber-modified epoxies, the adhesive plasticity must be included in order to correctly simulate the stress and strain distributions when the adhesive yields. Adherends can yield too, and the analysis needs to account for this behaviour if realistic failure loads are to be predicted. One of the most important works considering adhesive plasticity was done by Hart-Smith for SLJs^[18] and DLJs .^[19] Material nonlinearity due to plastic behaviour is not often included because of the increased complexity in the mathematical formulation and the need to resort to numerical methods (typically, FEA), for a solution. Only three analyses were found in the literature on analytical models that considered both adherend and adhesive nonlinear behaviour .^[20,30,38]

3. Process zone methods

In engineering applications, the presence of stress concentrations (also known as hot-spots) in design of adhesively bonded structures is almost unavoidable .^[47] These hot-spots can be generated from geometrical features (e.g. notches, holes, corners, edges) that are inherent characteristics on the construction of components .^[48] However, in the modelling of adhesive joints the consideration of only the peak (i.e. maximum) value of stress (or strain) in such hot spots, have been shown to produce over-conservative predictions .^[49]

Therefore, based on the concept that failure takes place within a process zone, process zone models (i.e. averaging methods) have been developed to address the effect of stress concentration on the mechanical behaviour of materials .^[50] This approach takes into consideration that fracture involving crack initiation and propagation is strongly influenced by aspects of the stress (or strain) field within the process zone, such as the gradient of stress, or the absolute volume of material under high stress.

The increasing use of process zone models is closely related to the evolution of computer-aided simulation (e.g. Finite Element Analysis – FEA), which made the calculation of stress very straightforward even for complex structures .^[51] In this regard, the use of standard FEA techniques, which do not require special purpose elements, is especially interesting for large scale models in industrial applications .^[52] These techniques employ standard elements (e.g. continuum elements) to accurately calculate the distribution of stress in bonded joints by considering, for instance, kinematic constraints between adhesive and adherends .^[53]

Another advantage of process zone methods is their ability to circumvent issues related to the high dependency of the mesh size when considering maximum value criteria, as well as capability to deal with singularities (e.g. overlap ends of single-lap joints) .^[12] Process zone methods have been investigated since the 50′s (e.g. Neuber^[54] and Peterson^[55]) and were brought back in the recent years by Taylor .^[56]

Process zone methods rely on a length parameter (i.e. critical distance) which defines the size of the zone (as shown in Figure 2). This process length can be defined using calibration techniques, as done by Hoey et al., ^[57] a technique that has been already validated for adhesives in fatigue conditions by Beber et al.^[58]

Figure 2. Theory of critical distances: (a) Point method, (b) Line method.

Then, the defined failure criterion (stress or strain) can be considered as an average value in a point, in a line, in an area or in a volume .^[59] For instance, by averaging an equivalent stress ${\sigma }_{eq}$ within a critical distance L, an effective stress ${\sigma }_{eff}$ is obtained which can be used for prediction purposes, as illustrated in Figure 2.

The TCD comprises several methods that consider the averaged stress in a point, in a line, in an area or in a volume. The TCD has been applied to material such as metals, ceramics, polymers (including adhesives) considering both linear-elastic, as well as plastic behaviour .^[59] In the Point Method (${\sigma }_{eff,PM}$) the effective stress is postulated to be the equivalent stress at a point determined by a critical distance ($L_{PM})$ from the singularity (or notch), that is:

(2) ${\sigma }_{eff,PM}={\sigma }_{eq}\left(L_{PM}\right)$(2)

In the Line Method (${\sigma }_{eff,LM}$), the effective stress is considered to be the equivalent stress averaged in a line with a length determined by a critical distance ($L_{PM}$) from the notch:

(3) ${\sigma }_{eff,LM}=\frac{1}{L_{LM}}\underset{0}{\overset{L_{LM}}{\int }}{\sigma }_{eq}\left(x\right).dx$(3)

On the other hand, Zhao et al.^[60] considered the adhesive thickness as the distance over which the maximum stresses are averaged and then compared to the adhesive yield strength. Clark and Mcgregor^[61] established that for the failure to take place the maximum principal stress must exceed the ultimate tensile stress of the adhesive within a finite process zone. Three joints types, namely single-lap, double strap and T-peel joints were used to demonstrate that length of such zone was independent of the joint geometry. Critical strain at a distance was used by Towse et al.^[62] as failure criterion the critical strain at a distance to predict the behaviour of double-lap joints. Authors have applied a nonlinear analysis, which included the effect of residual thermal stresses.

Regarding the static failure of adhesively bonded joints, Crocombe^[63] applied a stress criterion based on an averaging method in order to predict the failure load of single-lap joints considering a linear-elastic material behaviour obtaining good predictions for ductile adhesives. Critical longitudinal strain method has been employed for single lap joints to successfully predict the effect of overlap length and substrate thickness (Khoramishad et al.^[64]), bondline thickness (Akhavan-Safar et al.^[65]), and dissimilar adherends (Cruz et al.^[66]).

For the fatigue of adhesives, Schneider et al.^[67] used the maximum principal stress to predict the lifetime of scarf, thick adherend and single lap joints at different temperatures using a linear-elastic material model. Beber et al. extended this analysis to include elasto-plastic material models, ^[68] as well as to consider multiaxial stress response .^[69] These considerations were able to improve the quality of predictions, especially for SLJs. Schmidt et al.^[70] performed a reliable fatigue strength assessment of adhesively bonded thin‐walled steel structures taking into account both point and line averaging methods.

A practical example of a lifetime prediction carried out using a process zone method based on the critical distance is presented in Figure 3. The objective is to predict the lifetime associated with an operational load $F_{APP}$. A requirement that should be previously defined is the equivalent stress (e.g. maximum principal) and the critical distance L. Based on the operational load, the effective stress is calculated ${\sigma }_{a,eff}^{APP}$. Considering a reference sample, which SN curve is be used for definition of the lifetime prediction, a reference force $F_{REF}$ that causes the same effective stress is obtained from a simulation. From the reference force a stress amplitude is calculated, which is then finally correlated with a predicted lifetime that corresponds the lifetime on the sample to be predicted.

Figure 3. Application of the process zone method based on the use of the critical distance for lifetime prediction of adhesive bonded joints.

Beber et al.^[58] succesfully carried out critical distance-based fatigue lifetime predictions of three notched bulk specimens and a single-lap joint considering the Drucker-Prager equivalent stress as failure criterion. More recently, Sousa et al.^[71] have employed a process zone method based on an equivalent notch approach to consider mixed mode loading for the prediction of single-lap joints.

4. Linear elastic fracture mechanics

Linear Elastic Fracture Mechanics (LEFM) is a basic theory of fracture, originally proposed by Griffith^[72] and further developed by Irwin^[73] and Rice .^[74] This method is used to predict the critical load that leads to the propagation of an existing crack under a static load or the crack growth rate under cyclic loadings. It can be applied when the material behaves in a linear-elastic manner and the fracture process zone is included in the singularity region around the crack tip. In other cases, elastic-plastic fracture mechanics or Cohesive Zone Models can be used to describe the fracture behaviour.

The fundamental problem of the Griffith’s theory is shown in Figure 4. In this case, the critical remote stress leading to the unstable growth of a crack with a length 2a within an infinite plate can be calculated as:

(4) ${\sigma }_c=\sqrt{\frac{E{G}_{Ic}}{\pi a}}$(4)

Figure 4. Fundamental problem of the Griffith’s theory – a crack in a plate under tension.

where G_Ic is the mode I critical Energy Release Rate (ERR) and E is the Young’s modulus.

In a more general formulation, the crack propagation under a quasi-static load occurs when the ERR, G, is higher than a critical value G_c:

(5) $G=-\frac{d{\mathrm{\Pi }}_p}{dA}=G_c$(5)

where Π_p is the total potential energy of the body and $A$ the crack area. In general, the total ERR is the sum of mode I, II and III contributions, and its critical value depends on the mode mixity (ratio between single components).

After the work by Williams, ^[75] the LEFM was extended to interface cracks, so that is could be adopted for predicting, for instance, the propagation of cracks in bonded connections. The complete stress fields in the neighbourhood of an interface crack tip were presented by Rice, ^[74] showing that the singular behaviour with an exponent equal to −0.5 is preserved, even if a logarithmic term appears, thus generating oscillatory stress fields in the close vicinity of the tip. This leads to difficulties in the determination of the mode I and II components of the Stress Intensity Factors or the ERR (Rice,^[74]). Indeed, the mode decomposition, and therefore the values of G_I and G_II, are not univocally defined but depend on the reference length adopted for their calculation. Therefore, special care has to be taken when the ERR components are calculated through the Virtual Crack Closure Technique, VCCT (see the next section), as they depend on the adopted element size.

In addition, the theoretical displacement fields are characterised by an unrealistic crack face compenetration in the close vicinity of the tip. It was shown that in several conditions such feature is confined within a region with a negligible extension, so that the LEFM framework proposed by Rice^[76] remains applicable. On the other hand, when the compenetration zone is more extended, the LEFM can still be adopted, according to the framework proposed by Comninou .^[77]

As a further complication, the critical value of the ERR and the Paris-like law for the fatigue crack propagation are strongly dependant on the ratio between the mode I and II components, generally with a higher resistance to crack propagation as a pure mode II condition is approached (see, for instance, Wang^[78]; Ducept et al.^[79]; Carraro et al.^[80]; Adamos et al.^[81]). Several models were proposed in the literature for assessing the critical ERR or the Paris-like law in mixed mode conditions. For instance, the critical ERR can be estimated from the simple phenomenological expression proposed by Hutchinson and Suo^[82] or Ducept et al., ^[79] or with the damage-based model by Wang, ^[78] or with the analytical model by Adamos et al. .^[81]

Under fatigue loading, the crack propagation in bonded joints was sometimes shown to be well predicted by using only the mode I component of ERR and the relevant Paris-like law obtained through a pure mode I test (Kinloch and Osiyemi^[83]; Abdel Wahab et al.^[84]; Rocha et al.^[85]). Other authors made use of empirical expressions for defining an “effective” or “equivalent” ERR, incorporating the mode mixity, to be used in the Paris-like law for predicting the fatigue crack growth (Cheuk et al.^[86]; Quaresimin and Ricotta^[87]). A model, valid both for static and fatigue loading, was proposed by Carraro et al.^[88] based on the crack propagation mechanisms observed in composite bonded joints. The proposed criterion uses the only mode I ERR when the loading condition is mode I-dominated, and local stress parameter, derived from the maximum principal stress field in the adhesive, when the loading mode is shear-dominated. This last criterion, however, cannot be fully ascribed to the LEFM framework, as the stress parameter is calculated within a process zone that may extend beyond the singularity region.

Despite the mentioned complexities and difficulties, the LEFM was successfully applied to predict the crack propagation in bonded connections (see, for instance, some of the applications listed in the next section, as the VCCT technique is often adopted for the calculation of the ERR).

Special care must be taken when using the LEFM approach in bonded joints. Indeed, the fracture process zone can overcome the singularity region, leading to a significant dependence of the critical ERR on the adhesive thickness, under mode I and, mainly, mode II loadings (Ji et al., ^[89,90] Sarrado et al.^[91]). In these conditions, other approaches, such as the Cohesive Zone Model described later in the paper, can be applied, even if it requires a higher computational effort (Rabinovitch^[92]).

As a final remark, it can be said that LEFM can be adopted to predict the static or fatigue crack propagation in bonded joints, provided that the fracture process zone is small, which is typically the case when the adhesive is sufficiently brittle (Sun et al.^[93]). Critical ERR or Paris-like curves are preferably obtained on joints with the same adhesive thickness as the joint to be designed or verified. Crack initiation falls outside the scope of LEFM, hence that phenomenon cannot be accounted for by this method.

It should be emphasized that, despite its important advantage – simplicity and fast calculations, LEFM presents also some significant disadvantages. The idealizations inherent to the method limit its potential to analyse problems encountered in industrial practice. The method is only suitable for materials with linear elastic characteristics. Furthermore, the assumption of pre-existing cracks is crucial for LEFM, which excludes some types of cracks from the analysis (e.g., very small cracks).

The main challenges of the LEFM methods can be summarized as follows:

• to develop numerical tools for the fast computation of the SIF and ERR for interface cracks, with free and coarse mesh suitable to be applied to big structures;

• to define standardised testing procedures for the mode III properties, useful in practical applications;

• to model the influence of the adhesive thickness on the fracture toughness.

5. Virtual Crack Closure Technique (VCCT)

Virtual Crack Closure Technique (VCCT) is a linear elastic fracture mechanics-based method used to compute the strain energy release rate based on 2D and solid 3D FEA which provide the mode separation required when using the mixed-mode fracture criterion. The method has been developed by Rybicki and Kanninen^[94] and detailed by Krueger .^[95]

The VCCT method is based on the assumption that, when a crack grows, the energy released balances the work needed to close to crack to its initial length before propagation. The key input parameters of the method are the critical mode-I strain energy release rate G_IC and critical mode-II strain energy release rate G_IIC, which have to be determined experimentally. The VCCT method is usually applied through the FE method.

In the three-dimensional VCCT model, three stress intensity factors (SIF) are calculated using the following equation:

(6) $G_i=\frac{K_i^2}{E}\beta \left(i=1,2,3\right)$(6)

where $G_i$ is the energy release rate in mode i, $K_i$ the stress intensity factor for mode i,$E$ the elastic modulus, $\beta =1$ for plane stress and $\beta =1-{\nu }^2$ for plane strain whereas $\nu$ is the Poisson ratio. Thus, either for 2-dimensional or 3-dimensional problems, the previous equation for a crack extension Δc may be expressed as:

(7) $G_I=-\frac{1}{2\mathrm{\Delta }c}{F}_{22}\mathrm{\Delta }y$(7)

(8) $G_{II}=-\frac{1}{2\mathrm{\Delta }c}{F}_{11}\mathrm{\Delta }x$(8)

(9) $G_{III}=-\frac{1}{2\mathrm{\Delta }c}{F}_{33}\mathrm{\Delta }z$(9)

where G_I, G_II, G_III the energy release rates in mode I, II and III respectively, F₁₁, F₂₂ and F₃₃ the corresponding reaction forces on the crack tip and Δx, Δy and Δz the relative displacements between the upper and lower nodes of the crack as shown in Figure 5. In the particular case of the absence of one of the three modes, the corresponding equation may be neglected. Nevertheless, by using the VCCT method the corresponding G_i values should be imposed.

Figure 5. Schematic illustration of the concept of the VCCT method.

In the last decade, VCCT method has been implemented into commercial FE codes such as ANSYS and Abaqus. Initially, the method has been widely used to simulate delamination growth in CFRP composites^[95] and in the last decade it has found many applications to adhesive joints together with the Cohesive Zone Modelling method.

A short description of some representative works showing the evolution of VCCT in bonded joints follows.

Jokkinen et al.^[96] have studied the applicability of VCCT for crack growth analysis of an adhesively bonded joint with a ductile adhesive and self-similar crack growth. The analysis was performed on a Double Cantilever Beam (DCB) specimen with epoxy adhesive. First, a linear elastic behaviour was assumed for the adhesive and in the next step, a nonlinear analysis was performed by using a linear elastic – ideally plastic material model. For both approaches, the numerical results match well the experimental curves.

Senthil et al.^[97] used the VCCT method to ascertain the effective load carrying capacity of the adhesively bonded composite joints in the presence of closed debonds (defects). For the study, they used two adhesively bonded symmetric flat rectangular laminates of identical size in length and width having closed (embedded) debond of square shape at the specimen centre subjected to a uniform in-plane end-shortening. The authors investigated the influence of parameters such as laminate sequence, debond location, size and its shapes (square and circular).

Quaresimin and Ricotta^[87,98] developed a model for predicting the fatigue life of bonded joints, based on the distinction between a crack initiation and propagation phase. They used the VCCT, combined with a Paris-like curve linking the ERR to the crack growth rate, for predicting the life spent for the debond propagation in composite single lap joints, resulting in a sound agreement with experimental results.

Pirondi et al.^[99] developed a fatigue VCCT model based on the Direct Cyclic procedure of the Abaqus FE code. The model has been successfully applied to the DCB, ENF, and SLJ specimens and the results are in the form of crack length vs. number of cycles and G vs. crack length. The authors have also compared the fatigue VCCT model with a fatigue CZM model and have found that (a) the two models are in overall good agreement with each other, with the exception of the first instants of propagation where the CZM process zone has to shape up, whereas the VCCT starts with a sharp crack and (b) while the modelling effort is a bit higher (need of introducing a layer of cohesive elements), the CZM is easier to use (no need to identify the proper number of Fourier terms and time increment to represent cyclic loading). At the same time, it results more efficient as the computation is lower up to one order of magnitude, despite the less performing hardware used to run the analyses.

Many authors have combined the VCCT method with the CZM method. Jokinen et al.^[100] developed a method for modelling the fracture process with separate nucleation and propagation phases by combining the VCCT with the CZM on the finite element basis to take into account the development of fracture toughness. The method was applied to simulate a double cantilever beam (DCB) test as an example. The analysis focused on the physical validity of the VCCT-CZM coupling and on the determination of simulation parameter values. By using experimental data as a reference, the simulation results were compared to the results of traditional CZM and VCCT simulations. The comparison indicated that the combined CZM-VCCT method reproduced the DCB test cycles more accurately than the separate application of CZM and VCCT models.

Finally, recently, De Carvalho et al.^[101] have proposed a modelling approach to simulate delamination propagation in fatigue that combines VCCT with a progressive nodal release strategy. The progressive nodal release alleviates the artificial stress concentrations found when using the VCCT with instantaneous release to model 3D delamination without re-meshing. This enables crack shapes that do not conform to the underlying mesh to be readily simulated. Contrarily to previous implementations, the local delamination propagation direction is not assumed to follow mesh lines, but it is instead computed as part of the iterative procedure proposed to determine the maximum energy release along the crack front. The results suggest that VCCT can accurately simulate Mode I, Mixed Mode I/II and Mode II fatigue delamination growth and be used to simulate growth of arbitrarily shaped cracks that do not conform to the underlying mesh. The new approach can be potentially used to simulate fatigue debonding propagation in bonded joints.

The existing literature suggests that the VCCT is an efficient method to simulate debonding growth in bonded joints under static and fatigue loads. Of course, it has certain advantages and disadvantages and its implementation on the FE method is constantly improving. Overall, it stands below the CZM method.

In, ^[102] Floros and Tserpes, have compared the VCCT method with the CZM method in simulating debonding growth in a crack lap shear composite bonded joint containing a bolt as a crack stopper (Figure 6). They have found that the VCCT method is more capable than the CZM method in simulating sudden debonding.

Figure 6. (a) Geometry and boundary conditions of the crack-lap shear bonded/bolted specimen, ^[102] (b) Comparison of load-displacement curves predicted using the VCCT and CZM methods, ^[102] (c) Comparison of debonding growth simulated using the VCCT and CZM methods^[102]

6. Stress singularity approach

Adhesive bonding between dissimilar (or equal) materials gives rise to points where geometry and material properties change abruptly, which makes the stress state at these points unbounded from a linear elastic point of view. Examples of critical points of stress singularities in a double-lap shear joint are shown in Figure 7. Therefore, the use of failure initiation criteria based on nominal stresses is not a viable strategy to predict failure, whereas the use of singularity parameters to characterize the stress state is a feasible alternative, which is consistent with the principles of Linear Elastic Fracture Mechanics.

Figure 7. Examples of critical points where stress singularities appear at an adhesive joint.

A comprehensive review of the existing literature on stress singularities and prediction of failure initiation in adhesive joints falls beyond the scope of the present work. The main features of the method will be illustrated here as well as an historical perspective on its development.

From the pioneering works appeared in the early 20^th century, where the analytical solutions of 2-dimensional elasticity problems with stress singularities were first derived, (e.g. Wieghardt, ^[103] Williams^[104]), many advances have been made and published in this area.

One of the first classical attempts to use the singular stress state as the driving force to study the propagation of a crack running along the adhesive interface is due to Malyshev and Salganik, ^[105] who assumed an initial crack emanating from the corner at one end of the overlap length, and did not address the initiation of failure.

Two classical works, covering the analysis of the singular stress state without assuming a pre-existing crack or failure in bimaterial corners were studied by Bogy and Wang^[106] and Hein and Erdogan .^[107] It is important to notice that there are two clearly separated states, nucleation or initiation of the crack/failure, and the crack progression until the catastrophic failure.

Dempsey and Sinclair^[108] described in detail all types of stress singularities in 2D elasticity problems while a rigorous mathematical analysis of these solutions can be found in Nicaise and Sändig .^[109] In this scheme, and taking a polar coordinate system centred at the corner tip, the stresses admit, in the vast majority of cases, variable separation in the polar coordinates (r,$\theta$), which results in the following decomposition of the stress tensor:

(10) ${\sigma }_{ij}\left(r,\theta \right)\cong \sum _{n=1}^N{K}_n{r}^{{\lambda }_{n-1}}{f}_{ij}^{\left(n\right)}\left(\theta \right)$(10)

where, K_n (n= 1, …, N) are the Generalized Stress Intensity Factors (GSIFs), λ_n are the order of stress singularities, and f_ij($\theta$) is the characteristic angular function for stresses. The prefix “Generalized” refers to a Stress Intensity Factor which is associated to a singularity stress field different from that of a crack. While the order of stress singularities λ_n and the characteristic angular function f_ij($\theta$) only depend on the local geometry, local boundary conditions and the mechanical properties of the materials in the corner, the Generalized Stress Intensity Factors K_n depend on the global geometry and far field loading conditions of the problem. The quantities λ_n and f_ij($\theta$) can be obtained using analytical or semianalytical procedures, while the calculation of K_n requires detailed numerical methods, such as the Finite Element (or Boundary Element) Analysis. Alternatively, GSIFs can also be estimated from experimental data.

Very few particular configurations give rise to pure logarithmic stress singularities as pointed out by Sinclair .^[110]

One of the first proposals to predict the failure initiation in adhesive joints, without assuming a pre-existing crack is due to Gradin and Groth, ^[111] who defined a generalized stress intensity factor and compared it with a critical allowable value (the so-called generalized toughness) for the failure to initiate, as proposed in an analogous way by Hattori et al .^[112]

The growing adoption of composite materials, mainly in aeronautic lightweight structures, and the use of adhesive bonding with other materials complicate the analytical modelling of these singularity stress fields. Ting^[113] made a relevant contribution to the analysis of multi-material corners with the presence of anisotropic materials, which enabled the detailed study of joints with composite materials.

A recent comprehensive review of the mathematical aspects and practical applications of the singularity stress characterization in joints with anisotropic materials (including composites), has been presented by Mantič et al.^[114] Besides the overview of more than 120 reference on the topic, the authors introduced a failure criterion based on the singularity parameters of the stress state, together with corroborating experimental evidence.

The stress singularity approach has also been experimentally validated in the study of fatigue crack initiation (Lefebvre et al.^[115]), and fatigue strength in metallic bonded lap joints, Imanaka et al.^[116] Fatigue crack initiation in composite single lap joints was investigated by Quaresimin and Ricotta^[87] and Meneghetti et al., ^[117] thus proving the capability of the method to predict the effect of the overlap length.

Depending on the particular geometry of the corner and the elastic properties of the materials, the values of the orders of stress singularities λ_n vary between 0 (to make the stress singular) and 1 (to remain the strain energy with a finite value), and consequently, the dimensions of the GSIF vary accordingly. In many studies, a proper definition of the structure of the stress field is done by taking a characteristic length (L₀) to the power of the order of stress singularity in the denominator, equal than of the radial distance, so that the GSIF always has the same units as the stress. In any case, an appropriate definition of the GSIF is needed in the different corner configurations, to standardize its definition and to allow comparisons between different studies. A good example of this, can be found in Hwu and Kuo.^[118]

The use of singularity stresses to predict the failure of adhesively bonded joints is still the subject of intense research efforts. Examples of recent developments include the extension of analytical solutions to three-dimensional problems with anisotropic materials, (Ren et al ^[119] , Zappalorto et al.^[120]), which is quite useful to analyse joints with thick composite laminates as adherents. There are also applications in joints where adhesives with dissimilar stiffness are used along the overlap length in order to make the shear stress distribution more uniform and maximize the failure load of the joint (Breto et al.^[121]), and procedures to determine the critical value of the stress intensity factors^[122] in adhesive joints with stress singularities.

Singularity stress states influence the occurrence of premature failure in many problems, as well as in many standard test configurations, as pointed out recently by Barroso et al, ^[123] where at least three typical tests in composite materials with adhesives show significant higher failure loads after removing the local configuration which gives rise to the stress singularity. The excellent correlation between the predictions and the experimental evidence establishes confidence in the method.

The main advantage of this approach is the detailed and accurate knowledge of the local stress field at the corner tip, with analytical, or semi-analytical calculations for the key singular parameters of the stress representation, and with the knowledge of the contribution of each term in the stress decomposition defined in equation 10(10) ${\sigma }_{ij}\left(r,\theta \right)\cong \sum _{n=1}^N{K}_n{r}^{{\lambda }_{n-1}}{f}_{ij}^{\left(n\right)}\left(\theta \right)$(10) .

The main drawback of the stress singularity approach is the inherent complexity of the calculations required to derive the singularity parameters of the stress field, that is the order of stress singularities and the angular shape functions. The good accuracy achieved with existing computational tools in problems which involve isotropic materials, is no longer preserved in presence of anisotropy. Another inconvenience of this approach is the lack of standardized tests for the experimental determination of the critical GSIF value (Generalized Toughness) for a general corner configuration.

Additional aspects, such as the influence of local plasticity, local three-dimensional effects, or stress relaxation, have also been addressed in literature, although, as mentioned before, it is not the aim of this work to cover all these issues in detail, and just to make an overview of the approach.

Due to its nature, with the presence of predicted unbounded stresses, the Stress Singularity approach shows its full potential mainly in problems with no plasticity effects (fatigue problems with low load levels, or simply brittle materials). The great advantage is the explicit knowledge of the local stress field, whereas its main drawback is the complexity in the evaluation of the singularity parameters defining the stress state.

7. Finite fracture mechanics

Finite Fracture Mechanics (FFM) is a mathematical framework designed to predict crack initiation in presence of stress concentrations. It was first introduced in the pioneering work of Leguillon et al.^[124] for homogeneous and isotropic notched plates, and it has received considerable attention in the last two decades. FFM is based on the idea that failure does not initiate at a single point through the generation of a crack of infinitesimal length, but rather with the formation of a finite-size crack appears, characterized by a length a_c and a critical load, P_c, which are both unknowns of the problem. In the context of bonded joints, the principles of FFM can be illustrated by considering a bi-material corner as in Figure 8(a), with a remote applied load P. Crack initiation occurs with a finite crack length, a_c, when the load reaches its critical value P_c, as shown in Figure 8(b).

Figure 8. Bi-material corner a) before and b) after the initiation of a finite crack.

These two unknowns (a_c,P_c) can be calculated by imposing the following two conditions: a stress-based failure criterion, which requires that the stress in the pristine configuration, along the crack initiation length a_c, is at least equal to the material strength, and the energy criterion, which dictates that the released potential energy balances the energy required for the formation of new crack surfaces. These two constraints can be formally expressed as the following system of equations:

(11) $f\left(\sigma ,\tau ,{\sigma }_R,{\tau }_R\right)=1$(11)

(12) $-\frac{\mathrm{\Delta }{\mathrm{\Pi }}_p}{\mathrm{\Delta }\mathrm{A}}=G_c\left(W_p\right)$(12)

In equations 11(11) $f\left(\sigma ,\tau ,{\sigma }_R,{\tau }_R\right)=1$(11) (stress-based failure criterion) and 12 (energy criterion), $\sigma ,\tau$ are the normal and shear interface stress components, ${\sigma }_R,{\tau }_R$ are the normal and shear interface strengths, $W_p$ is the elastic potential energy, A is the crack area and $G_c$ is the critical total Energy Rate Release (ERR), which is typically a function of the mode mixity, $\mathrm{\Psi }$. It is worth highlighting that the symbol $\mathrm{\Delta }$ is adopted in the energy criterion, instead of δ as in the classical fracture mechanics, meaning that a finite energy amount is released in correspondence of a finite crack area increment $\mathrm{\Delta }\mathrm{A}$. A typical qualitative plot of the critical load estimated by the two criteria against the initial crack length is shown in Figure 9. The intersection of the two curves represents the solution of the system and identifies the crack initiation condition.

Figure 9. Qualitative plot of the stress and energy criteria.

The FFM approach was first applied to bonded joints by Leguillon et al.^[125] The case study consisted of a joint with steel adherends, bonded through an epoxy adhesive. Two strategies were adopted by the authors to predict the load at crack initiation: a macro-scale and a micro-scale approach. In the macro-scale approach, the existence of the bond-line was neglected, and the specimen was treated as a homogeneous steel plate with a sharp notch with an opening angle of 90° (square edge joint). In the micro-scale approach, instead, the steel/epoxy bi-material corner was considered. In both cases, the singular local stress fields and the ERR were calculated theoretically through an asymptotic expansion. A similar accuracy in the prediction of crack initiation load was obtained with the two methods, even if the micro-scale approach was the only one that could, in principle, describe the influence of the local geometry, such as the actual corner angle and the adhesive thickness. The micro-scale approach is also more complicated due to the complex nature of the stress fields in the presence of bi-material interface cracks (Rice^[74,76]).

Later, Garcia and Leguillon^[126] formalised the FFM approach for bonded connections between isotropic adherends made of the same material, with a generic corner angle. They used a multiaxial failure criterion (in the stress criterion) and a mode-dependency formulation for the fracture toughness (in the energy criterion). The macro-scale methodology was implicitly adopted, and the asymptotic expansion approach was used for the calculation of the singular stress fields.

Tran et al.^[127] adopted a similar macro-scale approach to predict the failure load of PC/epoxy/PC and PMMA/epoxy/PMMA joints, thus quantifying the influence of the corner angle with good accuracy.

Cornetti et al.^[128] and Weißgraeber and Becker^[129] used the FFM together with the Weak Interface (WI) approach for predicting crack initiation in concrete/composite and metallic bonded joints, respectively. In the WI approach, the bond-line is described as a bed of normal and shear springs, which makes the calculation of the pristine stress fields and the ERR particularly simple from an analytical point of view, even if the stress singularity in the presence of sharp corners is not accounted for. The proposed methodology was shown to be capable of describing the effect of the overlap length and the adhesive thickness in lap joints (Weißgraeber and Becker^[129]). As a further note, the effect of the local corner geometry (opening angle) cannot not be correctly captured with this approach due to the limitations of the simplified stress analysis.

Moradi et al.^[130] proposed a numerical procedure based on the FFM, in which the stress fields and the total ERR were calculated through FE analyses. The case of steel/epoxy/steel single lap joints was considered. The stress criterion involved the interface peel stress only and no mode dependency was considered for the critical ERR in the energy criterion. The approach captured the effect of the adhesive thickness and was also proved to capture the detrimental effect of the increasing adhesive Young modulus on the crack initiation load.

A similar approach was adopted by Mendoza-Navarro et al., ^[131] with the difference that both the peel and shear stresses were included in the stress criterion. The effects of the overlap length and multiaxial loads was correctly captured for steel/polyester resin/steel joints.

Later, Hell et al^[132] proposed a similar numerical procedure including also the possibility of crack initiation inside the adhesive, not at the interface. They applied the model to steel single lap joints, obtaining a good agreement with experimental results with different overlap length and adhesive thickness. Both the peel and the shear stresses were involved in the stress criterion and the mode dependency of the critical ERR was also considered.

Le Pavic et al.^[133] applied the FFM approach to aluminium joints tested through an Arcan fixture with different loading angles. They used both a simple analytical expression (as those adopted by Weißgraeber and Becker^[129]) and FE computations for the stress fields and the ERR. FEA showed very good correlation with experimental data, while the analytical method turned out to be more conservative, albeit much faster in its application.

Carrere et al.^[134] made a very interesting comparison between the FFM approach and a Continuum Damage Model (CDM), applied to steel single lap joints. Both the FFM and the CDM can be classified as coupled criteria, in the sense that both the interface strength and fracture toughness determine the structural failure. The difference is that, in the FFM the stress and energy criteria need to be simultaneously satisfied for the initiation of a finite size crack. Differently, in the CDM the satisfaction of the stress criterion triggers the degradation of a finite element, culminating in its suppression driven by the ERR. The two approaches provided quantitatively similar results when the CDM was used in combination with small-displacement FE simulations, which led to an underestimation of the experimental failure load. Although the two methods show a similar capacity to predict the occurrence of failure in bonded connections, FFM is more straightforward and more conservative than CDM. A drawback of FFM with respect to CDM is that the crack location has to be known a priori. The greater flexibility of CDM has, however, to be weighed against the higher computational cost of this method, due to the finer discretization of the adhesive layer (at least 20 elements).

All the examples reported so far considered 2D problems, which can be considered acceptable for simple geometries such as single or double lap joints. An extension of the FFM to 3D problems was formulated by Doitran and Leguillon .^[135,136] The stress fields and the ERR were calculated by FE analyses and the shape of the interface crack was chosen based on the iso-stress contour lines obtained from the analysis of the joint in the pristine condition. The 3D formulation was shown to slightly improve the failure predictions for scarf joints (Doitran and Leguillon^[136]).

To conclude, the FFM approach, initially conceived for homogeneous notched plates, proved to be a very useful tool for predicting the failure (meant as crack initiation) of bonded joints. The following main conclusions can be drawn:

• in principle, FFM makes possible to model the effect of the joint geometry, adhesive thickness, corner angle, material elastic properties and multiaxial loads in the prediction of the failure load;

• a macro-scale (i.e., neglecting the adhesive layer) or a micro-scale (considering the adhesive layer) approach can be adopted, the latter one being the only method to assess the influence of the adhesive thickness and tensile modulus, as well as the local corner shape;

• FFM requires the knowledge of the stress fields in the pristine configuration and the ERR for a crack initiated at the bi-material corner. They can be calculated analytically (through an asymptotic expansion of the singular stress fields or through simplified formulations), or numerically. Using the asymptotic expansion of the stress field or FEA with very fine local meshes are the only possibilities to account for the local corner geometry;

• the interface strength and fracture toughness need to be known. Standardised procedures can be found for measuring the interface toughness, whereas characterising the strength is more complicated because of the singular nature of the stress fields at a bi-material interface. A procedure was however proposed by Lauke and Barroso^[137] and Barroso et al.^[138];

• multiaxial failure criteria and mode-dependency formulations should be adopted in the stress and energy criteria, respectively;

• FFM can be applied to quasi-static loads only and it enables the evaluation of the crack initiation load, whereas it cannot predict the ultimate joint failure.

On the basis of these considerations, FFM results to be a very useful tool to make comparative and preliminary analysis based on plausible assumptions on interface properties, and to define optimal solutions in terms of geometry, constituents and loading conditions. The approach can be implemented in FE codes, even if the very fine meshes required for the computation of the stress fields and the ERR hinders its application to large structures, a limitation which can be circumvented by means of sub-modelling techniques.

As a last observation, to the best of the authors’ knowledge, no application of the FFM to composite bonded joints has been reported in the literature. Therefore, the potential of FFM to analyse composite bonded joints remains to be explored.

The main challenges of the FFM method can be summarized as follows:

• to assess the experimental standard procedures to obtain the input interface properties, and particularly, the interface strength;

• to investigate possible extensions of the method to the case of fatigue loading.

8. The cohesive zone method

8.1. Theory and application

Cohesive Zone Modelling (CZM) is a damage mechanics-based numerical approach suitable for modelling crack initiation and propagation. CZM is based on the works of Barenblatt^[139,140] and Dugdale^[141] in the late 1950s/early 1960s, when the idea of a cohesion zone in front of the top of a crack was proposed. Later, in 1976, it was numerically implemented by Hilleborg et al.^[142] and a numerical model which defines a function between traction and crack separation was proposed. Needleman^[143] proposed various polynomial and exponential functions for the definition of traction-separation law. Typical bi-linear traction-separation laws for mixed-mode load analyses are shown in Figure 10, where ${\delta }_I,{\delta }_{II}$, are the ultimate relative displacements, ${\delta }_I^0,{\delta }_{II}^0$ are the damage initiation relative displacements, ${\delta }_I^F,{\delta }_{II}^F$ are the ultimate relative displacements and $\beta ={\delta }_{II}/{\delta }_I$ is the “mode mixity”, which correlates ultimate normal and tangential relative displacements as can be seen in the small hatched triangle in Figure 10. $I$ and $II$ subscripts stand for normal and tangential direction, respectively Additionally, $T,S$ are the peak tractions in normal and tangential direction. Additionally, the area under the triangle for Mode I and Mode II is the critical energy release rate of the joint.

Figure 10. The bi-linear traction separation laws for a Mixed-Mode loading analysis^[144].

During the last decades CZM has evolved and it has eventually achieved the status of method of choice for simulating delamination of composite materials and onset and debonding growth of adhesive joints. The method can be implemented either through spring elements^[145,146] or by using cohesive elements between the adherends in 2-D^[147] or 3-D problems .^[148] The main reason for the wide acceptance of this method is its accessibility through the most common commercial Finite Element (FE) platforms and its capacity to be used in complex geometries.

CZM is based on the assumption that one or multiple interfaces can be artificially introduced in structures. Damage growth at these interfaces is modelled by the introduction of a discontinuity in the displacement field .^[149] As mentioned, the method is applied through traction-separation laws which can have various forms such as triangular, ^[150] linear-parabolic, ^[151] polynomial, ^[152] exponential^[153] and trapezoidal .^[154]

CZM approach has been successfully used by numerous authors for the simulation of delamination of CFRP materials as well as the debonding initiation and growth of joints with metal and/or composite adherents. Dias et al.^[155] applied CZM for simulating Mode I loading in a DCB specimen. De Moura et al.^[156] studied Mode II loading by means of experiments and numerical models, by applying a trapezoidal traction-separation law. Campilho et al.^[157] investigated the influence of the input parameters and the form of the traction-separation law. Alvarez et al.^[158] applied CZM approach for simulating DCB configuration and examine the CZM parameters’ sensitivity. As can be seen in Figure 11(a), numerical, analytical, and experimental results are in good agreement. Furthermore, Floros et al.^[144] performed an experimental study and numerical simulations using CZM approach to simulate Mode I, Mode II and Mixed Mode I+ II loading of composite bonded joints. Figure 11(b) depicts the comparison of the load-displacement curve of the DCB specimen predicted by the CZM method with the respective experimental curves. As can be seen, the CZM method predicts accurately the debonding behavior of the DCB composite bonded joints, except for the regions where unstable debonding growth occurred during the tests.

Figure 11. Comparison between numerical (CZM) and experimental load-displacement curves obtained for DCB specimens (a)^[158] (b)^[144].

For the most widely used bilinear traction-separation law, the main required inputs for the CZM approach are the initial stiffness of the joint, the maximum peak tractions and the critical energy release rate for normal and tangential directions. By some authors, ^[159] initial stiffness is treated as “penalty stiffness” while other authors^[160] suggest its definition by using interface thickness, $t$, and elastic moduli of the interface to calculate it. Peak tractions are measured by conducting tests of the bulk adhesive material and critical energy release rates are measured by performing fracture toughness tests under Mode I and Mode II loading. In most commercial FE codes, the ultimate relative displacements are required as input. Given that the area under the bilinear law for Mode I and Mode II is the critical energy release rate under Mode I and Mode II, respectively, ${\delta }_I^F$ and ${\delta }_{II}^F$ are calculated using the following equations:

(13) ${\delta }_I^F=\frac{2G_{IC}}{T}$(13)

(14) ${\delta }_{II}^F=\frac{2G_{IIC}}{S}$(14)

Apart from simulating delamination and debonding under pseudo-static loading, CZM has also been used to simulate fatigue crack growth, which has been the subject of intense research in the last decade. The simulating of mechanical behavior as well as debonding growth under fatigue loading is very important for the design of adhesive joints and the evaluation of their load bearing capacity.

The general idea which underpins the simulation of fatigue debonding growth is the modification of the bilinear traction-separation law by degrading the strength, as well as the stiffness and fracture energy of the cohesive elements as a function of the applied load cycles .^[161,162] For this purpose, a fatigue damage parameter ($d_f$) is introduced which accounts exclusively for the damage due to fatigue loading. The concept of the fatigue damage parameter and the way it is used to modify the bi-linear traction-separation law are illustrated in Figure 12. More specifically, ${\sigma }_{max}$ is the maximum traction, $d_s$ and $d_f$ are the damage parameters due to static and fatigue loading, $D_{tot}=d_s+d_f$ is the sum of the damage parameters. When $D_{tot}$ is equal to 1, due to the accumulated fatigue damage, debonding occurs and the cohesive element is deleted

Figure 12. The modified bi-linear traction-separation law for fatigue load^[163].

The fatigue load application in FE models is carried out mainly in two ways, either by applying the cycle by cycle or by applying a loading index strategy based on the maximum fatigue load. It is obvious that the first method is much more computationally expensive. Therefore, the second method of load application is preferable in simulating the debonding growth of adhesive joints as it is possible to simulate large number of loading cycles avoiding long-lasting analyses.

Concerning the degradation models, several authors proposed various approaches to perform simulations under fatigue loading, based on CZM principles. Turon et al.^[164] implemented a model in terms of damaged area $\left(dA/dN\right)$ for delamination simulation under high-cycle fatigue. Pirondi and Moroni, ^[165] studied Mode I and Mode II loading cases using a subroutine in which the strain energy release rate has been calculated and updated automatically. De Moura and Gonçalves implemented a formulation with a unique damage parameter accounting for cumulative damage resulting from static and fatigue loading for Mode I^[166] and Mode II .^[167]

Numerical simulation of adhesive joints under fatigue loading takes place by means of simulating the debonding growth or creating a Paris-like law. In most cases the numerical models are validated by comparing the numerical to the experimental results. Fernandez et al.^[168] simulated Mode I fatigue loading, Eklind et al.^[169] simulated high cycle fatigue growth and, Pirondi and Moroni^[170] proposed an improvement of a Cohesive Zone Model in order to account for fatigue delamination rate. Monteiro et al.^[171] studied exclusive Mode II fatigue loading. Additionally, several authors have been simulated fatigue debonding under Mixed-Mode I+ II loading using Compact Tension-Shear (CTS) specimen, ^[172] Short Beam Shear (SBS) specimen, ^[173] a Mixed-Mode I+ II bending apparatus^[174] and Cracked Lap Shear (CLS) specimen .^[175] Furthermore, more recently, Al-Azzawi et al.^[176] performed numerical and experimental investigations in Fiber-Metal Laminates under fatigue loading using a novel trapezoidal traction-separation law. Rocha et al.^[177] and Sousa et al.^[178] investigated Mixed Mode fatigue fracture behavior for different adhesive systems and studied the effects of mode mixity. Concerning joining of dissimilar materials, Choi and Kim^[179] studied the fatigue crack growth specimens with straight and penny-shaped cracks. Indicative results from the application of the CZM method to predict fatigue crack growth in the CLS specimen are given in Figure 13.

Figure 13. Predicted vs measured fatigue crack growth in a CLS specimen^[180].

CZM approach requires numerous input parameters which have to be determined experimentally. More specifically, critical energy release rate for Mode I and Mode II, peak traction in normal and tangential direction are necessary for the FE model, as well as ultimate displacement in normal and tangential direction. For the common case of bilinear traction-separation law ultimate displacements can be calculated by means of the area of the triangle by making use of the critical energy release rate and peak traction. Depending on the approach chosen for fatigue simulation, several input parameters may be required, such as the debonding growth as a function of elapsed cycles $\left(da/dN\right)$ and the energy release rate for various known number of cycles for Mode I and Mode II fatigue loading. Furthermore, a user-defined “frequency” is necessary for the correlation of the accumulation of fatigue damage with the fatigue cycles.

As a method for failure simulation of adhesively bonded joints, CZM has many advantages comparing to other techniques. The implementation of CZM in commercial FE codes is straightforward and it has the capacity to model complex geometries from coupon level to bigger structural parts in automotive and aerospace field. However, the method is mesh-dependent and fine mesh is necessary to simulate and predict the phenomena which take place during debonding growth. This drawback could lead to erroneous predictions in case of coarse mesh selection or large analyses time and computationally expensive FE simulations in case of very fine mesh size. For this reason, the user has to perform mesh-convergence study to optimize between solution accuracy and analysis time. Moreover, although the input parameters are not too many, most of them must be determined experimentally and there is a lack of standardized procedures for that purpose. A further limitation of the method is that, using the standard traction-separation law, sudden debonding growth cannot be predicted efficiently.

In spite of its popularity, algorithms for the simulation of fatigue debonding growth simulation codes have not yet been incorporated in most of the commercial FE tools. The user has to program his own routines in order to modify the traction-separation law and take into account the fatigue damage accumulation.

8.2. Application of the CZM method to glued-in rods in cross laminated timber

Literature efforts that take advantage of the CZM technique for the building construction field are the most general and spread. Furthermore, a growing trend is represented by the innovative application to timber joints for constructions .^[181–183] As an example, reference^[181] and Figure 14, illustrate how a standard pull-pull test configuration can be used for the mechanical analysis and characterization of the CZM input parameters. The key components for the reference specimen and FE model in Figure 14 are represented by a threaded rod that is glued into a CLT log (cross section centre).

Figure 14. Reference numerical model (ABAQUS) for the analysis and mechanical characterization of glued-in-rods in CLT applications (reproduced from^[181] with permission from Elsevier).

The basic modelling approach follows earlier applications of the CZM strategy to timber joints and components, ^[182,183] and it is further adapted to include the bonding effect of the interposed glue. A basic feature of the overall modelling approach in Figure 14 is that the separate constituent layers of CLT are described and connected by a rigid bond, while considering different material properties for each layer (crossed moduli). The CZM interaction takes place at the interface of the solid steel rod and the adjacent wooden elements, that are all described in the form of brick elements from the ABAQUS library. Such a modelling strategy is typically characterized by high computational cost, which makes convenient to take advantage from symmetry (as in Figure 14). The intrinsic risk for similar applications is the potential occurrence and propagation of local damage mechanisms (especially at the symmetry plane of rods and restrained CZM nodes) that could result in misleading structural predictions for the specimens to verify. In this regard, it is often convenient to analyze the full nominal geometry and avoid potential numerical issues in the region of rods and contacts .^[182]

Finally, the study in^[181] gives evidence of the sensitivity of numerical predictions to different arrangements of the rods, which could be glued parallel or perpendicular to the grain of the hosting wooden layers.

In general, the numerical investigation reported in^[181] proves that the use of CZM for the structural assessment of building structural elements can be efficient and accurate. Typical results and comparisons to earlier pull-out experiments can be found in Figure 15.

Figure 15. Typical numerical and experimental predictions for glued-in-rods in CLT applications: rods arranged (a) parallel or (b) perpendicular to the grain of timber (reproduced from^[181] with permission from Elsevier).

Besides, as a common aspect of the literature applications of the CZM technique to wooden joints and components (see for example^[181–183] for a selection of different joint typologies) major uncertainties and technical issues are still represented by the realistic and robust calibration of the key input parameters that are responsible of the CZM initiation and propagation. This criticality derives both from the high variability of metal joint solutions that can be used in timber structures, and also to the typically high sensitivity of the CZM damage parameters to failure mechanisms that are typical of an imperfect, orthotropic building materials like wood.

9. Progressive damage modelling

Progressive damage modeling (PDM) is a widely used technique for predicting the fatigue/breakage behavior and strength of bonded joints based on their damage state’s evolution. Several approaches for the simulation of damage progression in bonded joints have been introduced over the last decades. These approaches mainly focus on the prediction of the static strength of bonded joints based on different concepts, such as stress/strain-based failure criteria, ^{[148,184,185]} continuum mechanics, ^[186,187] cohesive zone models, ^[188–190] and fracture mechanics .^[191–195]

Harris et al.^[184] used a nonlinear FE technique to predict the failure mode and failure load of SLJs. They used a failure criterion based on the uniaxial tensile properties. Depending on the adhesive toughness, maximum stress or maximum strain criteria were used. The former was used for two uncured adhesives, while the latter was used for two cured adhesives. Apalak et al.^[185] investigated the initiation and propagation of damage zones in a bonded interlining joint’s composite panels. Panigrahi and Pradhan^[148] developed a 3D FE analysis to calculate the out-of-plane normal and shear stresses in a bonded liner joint (SLJ) with laminated FRP composite panels. This method is used to calculate local 3D stress fields in the most critical area that differ along the overlap length.

De Moura et. al.^[186] discussed the suitability of continuum damage models for the prediction of bonded joints’ mechanical behavior. They developed a continuum mixing damage model based on fracture characterization to accurately simulate the cases where the bond thickness plays an important role. Chen et al.^[187] proposed a failure criterion for progressive damage modeling of SLJ based on the criterion of strain energy. The criterion was used to model crack initiation and propagation in SLJ joints using a brittle adhesive and a ductile adhesive.

Xie et.al.^[192] presented a prediction methodology for hat-stiffened structures made of steel and CFRP composite material through an experiment to analyze and fail bonded structures. The virtual crack closure technique (VCCT) was used to calculate the crack tip’s strain energy release rates.

Cameselle-Molares et. al.^[195] developed a quasi-static PDM to predict the fracture behavior and strength of bonded FRP joints based on four steps, see Figure 16. The model was based on a mixed bearing fracture criterion^[196,197] as a function of a Master R curve resulting from experimental data from standard fracture mechanics connections. The developed failure criterion was dependent on the crack length and mode mix and considered the fiber bridge effect’s contribution. The energy release rate values for bonded DLJ joints were determined using the VCCT method in an FE model, and the numerically determined SERR was compared with the critical SERR given by the mixed-mode failure criterion .^[196]

Figure 16. Flow chart of PDM strength prediction methodology by Cameselle-Molares et. al.^[195].

One of the main challenges of PDM is its dependency on the experimental input. The second main challenge is its need for proper calibration of the method.

The algorithmic implementation of PDM is straightforward, and it usually gives accurate results. Once the method is verified and validated, it can be calibrated for different joint configurations with minimum effort to predict such joints’ behavior.

10. Probabilistic methods

Quality assessment of adhesive bonds is primarily based on physical testing as well as process control. Physical testing provides limited insights on the long-term behaviour of adhesive bonds because manufacturing processes and specimen preparation introduce several sources of variability which require careful consideration to be controlled, such imperfections in the bondline thickness and local concentrations of pores, kissing bonds and other defects. Furthermore, the loads applied to the bonded structures in the test environment might differ significantly from the loading scenarios found in service.

Material and component testing provide the necessary data to determine the parameters (i.e., to calibrate) of physics-based numerical models, which can then be used to simulate the occurrence of failure mechanisms and estimate the strength of adhesive bonds in virtually unlimited combinations of materials, geometric configurations and loading conditions. However, computational models are approximations of reality and the uncertainty in design variables, parameters and model structure must be properly accounted for to make reliable predictions. Although deterministic analysis might suggest that performance requirements are fulfilled (e.g. stress in the adhesive lower than the maximum allowable value), the inclusion of uncertainty in input variables and modelling assumptions might reveal that the probability of failure is yet significant (see^[198] for an example). The limitations of experimental and deterministic numerical methods motivate the introduction of probabilistic methods for the assessment of the structural performance of adhesively bonded joints.

The problem of quantifying the effects of uncertainty in the analysis of the mechanical behaviour of adhesive joints has been addressed in various forms. We identified three categories of approaches to facilitate the navigation throughout the existing literature, schematically illustrated in Figure 17:

Figure 17. Illustration of different approaches to model the effects of uncertain properties in the structural performance of adhesive joints: statistical modelling (a), structural reliability (b) and stochastic structural mechanics (c).

1. statistical modelling (primarily, data-driven)

2. structural reliability

3. stochastic structural mechanics (primarily, physics-based).

These categories should not be intended too rigidly, since hybrid combinations are also possible, as shown for example in, ^[198] where a structural model with adhesive thickness modelled as a stochastic field was combined with a reliability model to compute the probability of failure for adhesively bonded composite patches.

In the statistical modelling approach, the failure state of adhesive joints is characterized by a stochastically varying physical quantity such as load, stress or time which is modelled as a random variable $X$ with a continuous probability density function $p_X$. The parameters of the probability density function are estimated from experimental data generated from tests where the adhesive joints are loaded until failure. In the attempt to account for the limited control and information over real materials, geometric dimensions and loads, the target of the quantitative analysis shifts from the deterministic conditions for failure to the probability that such an event occurs. For an arbitrarily small increment of the failure-driving variable $dx$ (e.g. load) applied to a large population of $N$ nominally identical (e.g. in terms of manufacturing, age and loading conditions) joints, a certain number of failures will be observed. The resulting fractional change in the size of the population of functioning joints $dN/N$ can be identified (in the limit of arbitrarily large $N$) with the probability of failure for an individual joint randomly drawn from the population, that is $dN/N\approx p_{X}dx$, which represents an informal definition for the notion of probability density function $p_X$. The probability of failure at a specific value of $x$ is obtained by integration of the probability density function over the range of possible realizations of $X$ up to $x$, which defines the cumulative distribution function of the random variable $X$, denoted as $CDF\left(x\right)$. The most widely used statistical model in the literature on adhesive joints is the Weibull distribution, ^[199–206] which is defined by the probability density function:

(15) $p_X\left(x;\alpha ,\gamma ,x_0\right)=\left\{\begin{array}{cc}0& x\le x_0\\ \frac{\gamma }{\alpha }{\left(\frac{x-x_0}{\alpha }\right)}^{\gamma -1}{e}^{-{\left(\frac{x-x_0}{\alpha }\right)}^{\gamma }}& x>x_0\end{array}\right.,$(15)

where $x$ is a nonnegative variable which characterizes the failure state of the joint, such as critical load, ^{[199,205,206]} or a measure of stress in the adhesive at failure under tensile tests, shear tests^[200–203] or combinations thereof^[204]; $\alpha ,\gamma >0$ denote the scale and shape parameters, respectively, and $x_0\ge 0$ a shift factor which is sometimes omitted (see Figure 18 to visualize the shape of Weibull PDF and CDF for different values of the parameters).

Figure 18. Weibull probability density (a) and cumulative distribution function (b) for different values of the scale and shape parameters $\alpha ,\gamma$.

Weibull distribution finds broad application in the analysis of product reliability for a wide range of materials (e.g. metals, composites), failure mechanisms (e.g. fatigue, yielding, wear) and generally all the systems where the occurrence of failure localized at a small volume or component determines failure for the whole system (i.e. the weakest link theory). The reliability of a device such as a mechanical joint is defined as the probability that it performs without failure under stated operating conditions for a specified period of time .^[207] In general, the probability of safe performance (i.e. not failure), it is simply related to the probability of failure by $P_{nf}\left(\mathrm{x}\right)=1-P_f\left(x\right)=1-\mathrm{C}\mathrm{D}\mathrm{F}\left(x\right)$, which in the case of Weibull distribution takes the form

(16) $P_{nf}^W\left(x;\alpha ,\gamma ,x_0\right)=e^{-{\left(\frac{x-x_0}{\alpha }\right)}^{\gamma }}$(16)

All the parameters in equations 15(15) $p_X\left(x;\alpha ,\gamma ,x_0\right)=\left\{\begin{array}{cc}0& x\le x_0\\ \frac{\gamma }{\alpha }{\left(\frac{x-x_0}{\alpha }\right)}^{\gamma -1}{e}^{-{\left(\frac{x-x_0}{\alpha }\right)}^{\gamma }}& x>x_0\end{array}\right.,$(15) and 16(16) $P_{nf}^W\left(x;\alpha ,\gamma ,x_0\right)=e^{-{\left(\frac{x-x_0}{\alpha }\right)}^{\gamma }}$(16) are estimated from experimental data by means of standard methods such as maximum likelihood, Q-Q plots, linear regression. Experimental methods reported in the literature to investigate the reliability of adhesively bonded joints include the tensile and shear tests, ^[200] End-Notched Flexure test, ^[206] and shear-tensile interaction test .^[204]

The sensitivity of the estimated probability of failure to geometric parameters such as adhesive thickness and overlap length has been investigated theoretically in^[199] and experimentally in^[205] (for a single lap aluminium joint with acrylic adhesive) and^[204] (for double lap joints of pultruded Glass Fibre Reinforced Polymer profiles). The onset of failure is related to the number of randomly distributed defects inside the adhesive, such as voids, inclusions, and kissing defects. As the volume of the adhesive increases, more defects are likely to be included in the adhesive and therefore the strength decreases, while the stress concentration near the ends of the adhesive layer decreases. The competition between the volume effect (which favours lesser thickness) and the concentration of stress at the ends (which favours greater thickness) results in the optimal thickness for which the strength of the adhesive is maximum. Assuming that the failure stress within elementary volumes of adhesive follows Weibull distribution, the scale effect can be theoretically derived for an arbitrary volume of adhesive .^[199]

Fitting probabilistic models to failure data of adhesive joints is a versatile approach to estimate the probability of failure for virtually any type of material combinations and bond geometry. However, these models provide little understanding about the conditions for the occurrence of failure. Modelling strategies which do not rely solely on empirical models have to be introduced to provide deeper insights into the physics of failure mechanisms and the effects of design parameters on the resistance of adhesive joints to failure.

As an alternative to direct estimation from experimental data, the framework of structural reliability^[207] offers the possibility to compute the failure probability of adhesive joints on the basis of their geometry, material properties and loads. The failure condition is characterized by the definition of a limit state function $g\left(X\right)$, where $X$ denotes the set of design variables (e.g. material parameters, thickness, etc.), which might be all, or in part, specified as random variables with assigned probability density functions. A simple example of limit state function can be devised in terms of the stress state of the joint (expressed by some scalar measure $\sigma$, to not overcomplicate the notation):

(17) $g\left(\sigma ,{\sigma }_{max}\right)=\sigma -{\sigma }_{max}$(17)

where the strength ${\sigma }_{max}$, is the maximum stress that the joint is capable to withstand without failure. Equation 17(17) $g\left(\sigma ,{\sigma }_{max}\right)=\sigma -{\sigma }_{max}$(17) exemplifies a general implication of limit state functions, that is the partitioning of the design space into regions of safe and unsafe performance, determined by the inequalities $g\le 0$ and $g>0$, respectively. The stress state (or, equivalently, the external load) as well as the strength of the joint can be modelled as random variables with probability density functions estimated from failure data, as in .^[208] Alternatively, a model-based characterization of the stress state is also possible .^{[198,209–211]} For example, the analytical model by Volkersen was considered in^[198] to quantify the shear stress inside the adhesive in composite-aluminium patches with randomly varying adhesive thickness. The maximum shear stress was then extracted from the numerical solution of the model and inserted in equation 17(17) $g\left(\sigma ,{\sigma }_{max}\right)=\sigma -{\sigma }_{max}$(17) as representative measure of the stress driving the occurrence of failure in the adhesive joint.

For a given failure mechanism, the associated limit state function is amenable of several parameterizations, which might include linear or nonlinear combinations of stochastic and deterministic variables. Independently of the specific form of the limit state function, the probability of failure in the framework of structural reliability is evaluated as

(18) $P_f=\underset{g\left(X\right)\le 0}{\int }{f}_X\left(X_1,X_2\dots \right)dX,$(18)

where $f_X\left(X_1,X_2\dots \right)$ is the joint probability density function of all the uncertain input variables and the integral is computed over the failure region of the design space. The multi-dimensional integral in equation 18(18) $P_f=\underset{g\left(X\right)\le 0}{\int }{f}_X\left(X_1,X_2\dots \right)dX,$(18) can be evaluated numerically via sampling techniques such as Monte Carlo simulations, ^[207] or via semi-analytical methods such as FORM, which approximates the failure probability with $P_f\cong \mathrm{\Phi }\left(-\beta \right)$, where $\mathrm{\Phi }$ is the standard normal distribution and the reliability index $\beta$ is the solution of a constrained optimization problem in the design space .^[207]

Both statistical and reliability-based methods target the probability of failure as the main quantity of interest to characterize the structural performance of adhesive joints. Reliability and quality are closely related. Normally, quality focuses on the prevention of defects during the warranty phase, whereas reliability looks at preventing failures throughout the life cycle of the joint, that is from manufacturing to the end of service life. The derivation of safety margins for the design variables is crucial for the development of robust procedures for quality control and certification of load-carrying bonded structures. Safety margins quantify the range within which design variables may vary without failing to fulfil requirements such as “probability of failure less than 1%”. For example, the problem of random variations of adhesive thickness along the bondline was addressed in^[198] from a reliability perspective. In order to derive the safety margin $C_s$ which could bound the probability of failure for the joint below 1%, the adhesive thickness $e_a^{C_s}\left(\mathrm{x}\right)$ was expressed as

(19) $e_a^{C_s}\left(\mathrm{x}\right)=C_s{e}_a\left(\mathrm{x}\right),$(19)

where $e_a^{C_s}\left(\mathrm{x}\right)$ is local coordinate along the bondline and $e_a\left(x\right)$ the nominal design thickness. A parametric study was conducted by evaluating the right-hand-side of equation 17(17) $g\left(\sigma ,{\sigma }_{max}\right)=\sigma -{\sigma }_{max}$(17) for different safety margins, which enabled to determine the value of $C_s$ corresponding to a maximum probability of failure of 1% (i.e. $C_s\cong 1.54$). A practical indication that the authors of^[198] could draw from this analysis was that, in order to satisfy the safety requirement of $P_f<1\mathrm{\%}$, the manufacturing process of the joint should be designed to make the adhesive thickness as close as possible to $e_a^{C_s}\left(\mathrm{x}\right)$ rather than to the nominal value $e_a\left(\mathrm{x}\right)$.

As long as failure criteria can be expressed in terms of a limit state function which partitions the design space in safe and unsafe regions, the reliability of any type of adhesive joint can be quantified within the framework of structural reliability. The impact of the uncertain design variables on the probability of failure can be assessed through the evaluation of “importance factors”, which can often be done efficiently within the same computational framework used for $P_f$, thus providing a quantitative criterion to reduce the number of random variables to be considered in the analysis. The probability of failure as computed through equation 19(19) $e_a^{C_s}\left(\mathrm{x}\right)=C_s{e}_a\left(\mathrm{x}\right),$(19) incorporates the available information from experiments as well as theoretical models.

The role of experimental data differs between statistical (i.e. phenomenological) and structural (i.e. theoretical) models. While in statistical analysis the data drive the selection of the model to estimate the failure probability, structural models are derived from fundamental physical principles and data are used for calibration and validation purposes. Deterministic models acquire a stochastic character by introducing random variables to model parameters, loads or boundary conditions which are expected to be significantly affected by uncertainty due to intrinsic variability or lack of information. Mechanistic or physics-based modelling of adhesive joints under uncertainty is based on stochastic partial differential equations, which are discretized and solved numerically in space, time, and probability domains .^[212] For example, the following stochastic partial differential equation was derived in^[213] to assess the load-bearing capacity of adhesive anchors embedded in concrete under the effect of material degradation due to creep deformation, which makes the adhesive shear stiffness and strength uncertain:

(20) $\frac{\pi \left(D+t\right)}{t}\mathrm{G}\left(z,\omega \right)\mathrm{w}\left(z,\omega \right)-\mathrm{E}\mathrm{A}\frac{d^{2}w\left(z,\omega \right)}{d{z}^2}=0,$(20)

where $\mathrm{w}\left(z,\omega \right)$ is the longitudinal displacement of the anchor (modelled as a rod with diameter $D$, Young modulus $E$ and cross section area $A$), $z$ the coordinate along the rod axis, $\omega \in \mathrm{\Omega }$ sample space, which is the set of all possible realizations of the random variables, $t$ the adhesive thickness and $\mathrm{G}\left(z,\omega \right)$ the position-dependent adhesive shear modulus. The specification of $\mathrm{G}\left(z,\omega \right)$ and boundary conditions for equation 20(20) $\frac{\pi \left(D+t\right)}{t}\mathrm{G}\left(z,\omega \right)\mathrm{w}\left(z,\omega \right)-\mathrm{E}\mathrm{A}\frac{d^{2}w\left(z,\omega \right)}{d{z}^2}=0,$(20) determine the displacement $\mathrm{w}\left(z,\omega \right)$ and, therefore, the shear stress random field $\tau \left(z,\omega \right)$, which is compared to the allowed maximum shear stress to identify regions of safe and unsafe performance in the design space.

Further applications of stochastic structural mechanics to adhesive joints include the study of crack propagation in a Double Cantilever Beam setup with random interface stiffness and elongation at break, ^[214] and the sensitivity analysis of Finite Element Models of adhesively bonded metals with respect to uncertain cohesive zone parameters .^[215]

Probabilistic methods can be efficiently combined with analytical and numerical techniques to predict the occurrence of failure in adhesive joints, thus providing an assessment of the robustness of the conclusions drawn from the analysis, which is particularly relevant for the purpose of certification of novel structures. Reliability-based design offers a quantitative framework where the feasibility and safety of adhesive joints can be assessed at the very early phases of product development, opening wide opportunities for performance and cost optimization.

11. Additional considerations

Adhesive joints between dissimilar materials, made with adhesives curing at temperature have, at least, two identified sources of residual stresses associated with temperature. One source is associated to the chemical shrinkage during the curing stage, and the second source is associated with the different Thermal Expansion Coefficient (TEC) of the materials being joined. A comprehensive review of the different sources of stresses in adhesive joints is detailed in .^[216]

While some authors^[217] indicate that it is neither easy nor realistic to separate the stress changes due to the curing process from purely thermal effects, others^[218] argue instead that it is possible to clearly identify each of the two phenomena separately.

In the case of unidirectional long fibre reinforced materials, it is important to understand that for a general stacking sequence, residual stresses will appear inside the laminate, due to the different thermomechanical properties associated to the different orthotropic directions of the unidirectional reinforced material. Nevertheless, considering the laminate as an equivalent homogeneous orthotropic material, thermal effects will also appear if the equivalent TEC of the laminate is different than the TEC of the other material of the joint. Therefore, it is of paramount importance to calculate the equivalent TEC of the laminate from the thermomechanical properties of the individual plies of the laminate .^[219,220] An excellent study of the influence of temperature on the strength reduction of adhesive joints with composite materials is the classical work by Hart-Smith, ^[221] in which the main mechanisms contributing to the strength reduction of the joint are clearly identified and quantified: a) difference in the TEC of the materials, and b) values of the curing temperature of the adhesives.

Other important aspects when analysing temperature effects in adhesive joints are the thermal properties to take into account for the adhesive and the reference temperature value from which these thermal effects account from.

Characteristic temperatures polymeric materials such as adhesives are^[34] the glass transition temperature (T_g), which is particularly important for structural adhesives, the melting temperature (T_m), which has significance for partially crystalline polymers and negligible effect on noncrystalline or amorphous polymers (e.g. epoxides, phenolics, and acrylics), the decomposition temperature (T_d) at which the adhesive undergoes chemical degradation, and the curing (T_c) and room temperatures (T_amb), which are typical reference values for the thermal effects in the cooling stage of the curing process, when residual thermal stresses are generated (see^[222] for further information).

A still debated question regards the temperature jump to take into account for the thermal residual stresses during the cooling stage of the curing process, that is whether to consider the glass transition temperature^[223] (ΔT = T_g-T_amb) or the curing temperature^[224] (ΔT = T_c-T_amb). The difference between Tg and Tc in some adhesives can be relatively large, ^[222] which makes some authors^[225] argue that the reference temperature should be experimentally determined by means of a bi-material laminate which presents curvature after the curing stage; in that case, a gradual heating until the recovery of the flat original configuration would determine the more appropriate reference temperature.

Considering the adhesive as an isotropic material, the stress-strain law with temperature dependent properties is:

(21) ${\sigma }_{ij}=2\mu {\epsilon }_{ij}+\lambda {\epsilon }_{kk}{\delta }_{ij}-\left[\left(3\lambda +2\mu \right)\underset{T_0}{\overset{T}{\int }}\alpha \left(T{ }^{\mathrm{\prime }}\right)d{T}^{\prime }\right]{\delta }_{ij}$(21)

Where μ is the shear modulus, λ is the Lame constant and $\alpha$ is the TEC.

An accurate evaluation of the stresses developed during the cooling stage of the curing process should take into account the detailed temperature dependency of the TEC $\alpha$(T), but if only the final situation (at T) is of interest, equation (21)(21) ${\sigma }_{ij}=2\mu {\epsilon }_{ij}+\lambda {\epsilon }_{kk}{\delta }_{ij}-\left[\left(3\lambda +2\mu \right)\underset{T_0}{\overset{T}{\int }}\alpha \left(T{ }^{\mathrm{\prime }}\right)d{T}^{\prime }\right]{\delta }_{ij}$(21) can be reduced to the simpler relation:

(22) ${\sigma }_{ij}=2\mu {\epsilon }_{ij}+\lambda {\epsilon }_{kk}{\delta }_{ij}-\left(3\lambda +2\mu \right){\alpha }_{avg}\mathrm{\Delta }T{\delta }_{ij}$(22)

Where α_avg is the average value of the TEC in the temperature range (ΔT = T-T₀).The approximation that leads to equation (22)(22) ${\sigma }_{ij}=2\mu {\epsilon }_{ij}+\lambda {\epsilon }_{kk}{\delta }_{ij}-\left(3\lambda +2\mu \right){\alpha }_{avg}\mathrm{\Delta }T{\delta }_{ij}$(22) from equation (21)(21) ${\sigma }_{ij}=2\mu {\epsilon }_{ij}+\lambda {\epsilon }_{kk}{\delta }_{ij}-\left[\left(3\lambda +2\mu \right)\underset{T_0}{\overset{T}{\int }}\alpha \left(T{ }^{\mathrm{\prime }}\right)d{T}^{\prime }\right]{\delta }_{ij}$(21) tends to overestimate the residual stress values, as pointed out by^[226] and also observed experimentally in .^[227]

In^[227] two simple experiments were conducted to check the potential residual stress field developed in the cooling stage of the curing process. A composite-metal single-lap joint was bonded using an adhesive curing at T_c = 115°C. The thermal stresses in the critical part of the joint (where failure starts) are numerically calculated by means of FEA for different ΔT. The stress fields are completely different when considering the expected ΔT = T_c-T_amb = 115–25 = 90°C than those obtained when considering lower ΔT values. This lower ΔT are completely equivalent to testing the joint at temperatures above the room temperature, and those test were carried out in^[227] finding absolutely no significant influence in the failure load or failure path (considering detailed failure path observations at the critical part of the joint). These experimental findings indicate that the expected residual stress profile has not fully developed, maybe due to the relaxation effects mentioned in^[225] in presence of high stress gradients and to the viscoelastic behaviour of the adhesive.

Another experiment carried out in^[227] investigated time-dependent relaxation effects in 24 samples of the same single-lap joint configuration which were bonded and tested at equal intervals of time (in a logarithmic scale) from just 15 minutes after the curing stage (the necessary time to prepare the sample) to a 7 months later. The two bonded materials had different TEC and the adhesive was cured at temperature, letting the expected residual stress state to develop. The test results did not show any significant variation, see^[227,228] for further details and information. No relaxation was observed after 15 minutes from the end of the curing process. If any relaxation effect takes place in the adhesive layer, it might occur during the cooling stage of the curing process.

In summary, it seems that thermal residual stresses that could appear in the cooling stage of the curing when using an adhesive curing at temperature are overestimated if no relaxation effects are taken into account.

Different techniques have been used in bonded joints to measure (direct or indirectly) thermal stress fields. Photoelasticity^[229] was successfully used in bonded structures, as well as Moire interferometry^[230,231] giving full displacement fields, and strain gauges, ^[232,233] which provide only local information. Indirect measurements, using nonsymmetric laminates^[216,224] were also used to determine stress free reference temperature in bonded joints, or embedded optical fiber Bragg grating sensors^[232] for measurements inside the adhesive layer. More sophisticated methods such as Neutron diffraction^[234] and X-ray diffraction^[235] also allows internal and surface (or near surface) strain measurement, respectively. Also, Digital Image Correlation has been used to experimentally measure the residual strain/stress fields .^[236]

Theory and experiments show that adhesive joints between dissimilar materials generate residual thermal stresses, which increase with the thermal expansion coefficient mismatch and with the curing temperature value. Although it is important to take into account these thermal stresses for failure initiation prediction, the complexity of different relaxation mechanisms in the adhesive complicate the exact determination of the residual thermal stresses, especially in the long-term.

12. Concluding remarks

In the present paper, a review of some of the most widely used failure theories and models for adhesive joints has been performed. The evaluation has been performed by experts from different scientific areas, so it represents a more generalized viewpoint. Table 2 summarizes the evaluation findings for the different methods that have been reviewed. Thus, it may be used as a guide for a quick selection of the method.

Table 2. Evaluation summary of the different failure theories and models.

Modelling approach	Material input data	Results/application		Advantages	Limitations
		CI-CP-FF*	S-F-D**
Classical Analytical Methods	- Elastic or elastic/plastic properties - Critical peel and shear stress values	CI	S+F	- Simple and fast calculations - Estimation of the effect of the macro-scale geometry - Account for material non-linearity (only some models)	- No through-the-thickness stress distribution - Neglects the effects of local corner geometry
Process Zone Methods	- Materials’ constitutive laws. - Adhesive/interface strength - Adhesive/interface S-N curve (fatigue) - Critical length	CI or CP	S+F	- Simple application through FE analysis - No mesh convergence issues - Valid also for non-LE behaviour	- Non-trivial definition of the critical length - No damage evolution prediction when conceived for CI
LEFM	- Elastic properties - Critical ERR/SIF - Paris-like law (fatigue)	CP+FF	S+F+D	- Simple and fast calculations - Simple evaluation of the structural load bearing capability in the presence of detected damage - Once coupled with a crack initiation criterion, allows the development of a comprehensive predictive framework	- Applicable for LE behaviour only, with small process zones - Requires the calibration of the fracture parameters for a specific adhesive thickness
VCCT	- Elastic properties - Critical ERR/SIF - Paris-like law (fatigue)	CP+FF	S+F+D	- It does not require too fine meshes - Same advantages as LEFM	- Same as LEFM
Stress Singularity Approach	- Elastic properties - Critical GSIF or GSIF-life curves	CI	S+F	- Accounts for corner stress singularity and, thus, the effect of local corner geometry and adhesive thickness	- Different critical GSIFs for different corner geometry - Complex stress fields calculation - Limited to linear elastic behaviour
Finite Fracture Mechanics	- Elastic properties - Interface strength and fracture toughness	CI	S	- It predicts the influence of the global and corner local geometry, adhesive thickness and elastic properties - Useful in a preliminary design phase and for comparing different solutions. - It requires LE analyses only, though with significant local mesh refinements	- LE behaviour only - The input interface strength is of non-trivial characterisation - Labour intensive, mainly for complex components/geometries
CZM	- Materials’ constitutive laws - Mode I and II cohesive laws - Mode I and II Paris-like laws	CI+CP+FF	S+F+D	- Applicable in case of non-LE behaviour and in case of large process zone - Available in most commercial FE codes, at least for static loadings	- Challenging identification of the cohesive laws - Dependence on the cohesive law type - Mesh dependence issues - Computationally expensive
Statistical modelling	Parameters of probability distribution, e.g. for Weibull model: $\alpha$, $\gamma$ and $x_0$	FF	S + F + D	Flexible to fit different types and moderately large of failure datasets Safety factors for design variables can be derived from the reliability function of the joint	Limited insights about the physics of failure
Structural reliability	Limit state function $\mathrm{g}\left(\mathrm{S},\mathrm{L}/\sigma \right)$, where $\mathrm{S}$ is the capacity of an adhesively bonded joint to resist the external load $\mathrm{L}$/internal stress $\sigma$. The $\mathrm{g}$ function depends on material and geometric properties which might be deterministic or random with assigned probability distributions	FF	S + F + D	No assumptions about the statistical properties of key variables for failure Safety factors for design variables can be derived from the reliability function of the joint	Trade-off between accuracy and computational resources to evaluate the failure probability integral (e.g. via Monte Carlo simulations)
Stochastic structural mechanics	Varying with specific models. All or part of material and geometric properties (bulk and interface) are modelled as random variables with assigned probability distributions	CP, FF	S+F+D	Detailed insights into failure mechanisms Potentially high predictive power	Substantial modelling expertise, computational resources, model calibration and validation to secure reliability of the predictions

* CI = Crack Initiation, CP = Crack Propagation, FF = Final Failure

** S = Static, F = Fatigue, D = Dynamic loading

Nomenclature

τ_xy:adhesive shear stress (SLS specimen)

F:Axial tensile force (SLS specimen)

b:Adherent’s width (SLS specimen)

L:Adherent’s length (SLS specimen)

σ_eq:Equivalent stress

σ_eff:Effective stress

L_PM:Critical distance

σ_eff,PM:Effective stress at the critical distance

F_APP:Operational force

F_REF:Reference forcé

a:Half crack length

E:Young’s modulus

G_IC:Mode I critical energy reléase rate

σ_c:Critical remote stress

Π_p:Total potential energy

$A$:Crack area

K_iStress intensity factor for mode i

G_i:Energy release rate

G_I:Mode I energy release rate

G_II:Mode II energy release rate

G_III:Mode III energy release rate

F₁₁, F₂₂, F₃₃: Corresponding reaction forces on the crack tip

Δx, Δy, Δz: Relative displacements at the crack tip

K_n (n = 1, …, N):Generalized Stress Intensity Factors (GSIFs)

f_ij($\theta$):Characteristic angular function for stresses

λ_n:Stress singularities

W_p:Elastic potential energy

σ_R:Normal interface strength

τ_R:Shear interface strength

Ψ:Mode mixity

X:Random variable

p_x:Probability density function

N:Population

α, γ>0:Scale and shape parameters

x₀:Shift factor

P_nf(x):Probability of failure

$P_{nf}^W$:Weibull distribution

g(x):Limit state function

f_X:Joint probability density function

C_s:Safety margin

$e_a^{C_s}$:Local coordinate along the bondline

w(z,ω):Longitudinal displacement

G:Adhesive shear modulus

K:Generalized Stress Intensity Factor (GSIFs)

λ_n:Order of stress singularity

f_ij($\theta$):Characteristic angular function for stresses

L₀:Characteristic length

T_g:Glass transition temperature

T_m:Melting temperature

T_d:Decomposition temperature

T_c:Curing temperature

T_amb:Room temperature

$\alpha$:Thermal expansion coefficient (TEC)

$\alpha$_avg:Average value of the thermal expansion coefficient

$\mu$:Shear modulus

λ:Lamé constant

Acknowledgements

This publication is the outcome of a collaborative work that has been performed in the frame of the COST Action CA18120 (CERTBOND - https://certbond.eu/), supported by COST (European Cooperation in Science and Technology - https://www.cost.eu/).

Correction Statement

This article has been republished with minor changes. These changes do not impact the academic content of the article.

Additional information

Funding

This work was supported by the European Commission [COST Action CA18120-CERTBOND].