Progressive failure mechanism of laminated composites under fatigue loading

Abstract

A cohesive zone model for delamination propagation in laminated composites under static and fatigue loading has been derived and validated with experimental data under different mode conditions. This study presents a new approach to quantify fatigue delamination degradation based on damage mechanics to evaluate the rate of fatigue damage ( $\partial \mathrm{D}/\partial \mathrm{N}$ ). The static damage evaluation and fatigue damage degradation are derived from damage surface concept. Both static and fatigue damage linked each other to establish fatigue crack growth formula in the laminated composites. A user-defined subroutine, UMAT, has been employed to develop and implement a damage model in ABAQUS. Two different specimens; a double cantilever beam and a single lap joint were used to investigate the effectiveness of the new method. The simulation results revealed that the developed model had good agreement with experimental data available in literature.

Introduction

Fatigue-driven delamination in laminated composites is one of the general failure modes in composite structure e.g. airplanes, automobiles and wind turbines, all of which undergo cyclic loading. ¹ Fatigue is also a major contributor to failures in rotating machinery and has been addressed in many research work.^{2, 3} Life estimation of composite materials is considered a serious challenge for engineers because of the uncertainties associated with damage initiation sites and propagation direction in these composites. There is a real and urgent need for reliable analytical methods to reduce the time and cost of experimental testing, and to enable optimally designed structures. ⁴

In composite structures, it is important to detect the damage at an early stage of failure and to know how the damage will grow during the service life of structures. A relationship between fatigue crack growth and stress intensity factor is generally written as a power law (i.e. the Paris law) for metals and composite materials. ⁵

Different approaches have been employed to analyse fatigue delamination growth, such as linear elastic fracture mechanics e.g. the virtual crack closure technique (VCCT), and cohesive zone models. Although the VCCT technique is widely used to simulate fatigue crack propagation, it has some limitations regarding crack path and re-meshing requirements. An alternative to the VCCT to simulate interface crack propagation, such as delamination in composite materials, is the cohesive zone model.

Over the last twenty years many researchers have used finite element methods, including the cohesive zone model, rather than attempt analytical solutions to simulate crack growth

Maiti and Geubelle ⁶ introduced the instantaneous interface stiffness degradation law under cyclical fatigue loading. The rate of change of cohesive stiffness was defined as a function of the number of cycles. They proposed a power law for the cohesive model that can be expressed in terms of discretized time steps. While the instantaneous cohesive stiffness was assumed by Serebrinsky and Ortiz ⁷ to degrade exponentially with each unloading-reloading cycle.

Load envelope methods have been suggested by many researcher^8–14 to simulate fatigue damage growth based on experimental Paris law curves.

Gornet and Ijaz ¹⁵ investigated high-cyclic elastic fatigue damage model for carbon fiber epoxy matrix laminates. A classical interface damage was used to predict the fatigue damage parameters for different mode-mixtures.

Nojavan et al. ⁴ proposed a non-Paris law based fatigue cohesive zone model. The authors assumed a simple power law for fatigue delamination growth under pure Mode I and Mode II loading. Although their proposed model can predict mixed mode delamination growth from simple parameters, their research demonstrated that the suggested cohesive zone model can not be used for arbitrary damage onset and propagation. Therefore, a crack path must be pre-defined in order to arrange the elements along the crack path in the finite element model.

Fatigue loading

Composite structures can be subject to cyclic loads with different wave shapes such as sinusoidal, square, etc. The load fluctuates with a certain frequency (fn) with a constant amplitude and produces a periodical stress or strain. Figure 1 shows typical sinusoidal and square loading curves with relevant parameters i.e. average stress is ((σ_max+σ_min)/2), stress amplitude is (((σ_max−σ_min))/2), and R-ratio (σ_min/σ_max). High amplitudes of fatigue loading reduce the remaining life of the structure. ¹⁶

Figure 1.

Sinusoidal and square stress functions.

In fracture mechanics, the crack growth rate depends on the stress intensity factor or the energy release rate, and propagates according to the Paris law. In the last few years many researchers have introduced algorithms into the cohesive zone model which include fatigue damage degradation. Belnoue et al. ¹⁷ employed the Paris law to represent fatigue in an interface element based on change in strain energy release rate, ΔG. The crack growth rate, $\partial \mathrm{a}/\partial \mathrm{N}$ , was assumed to be a function of change in strain energy release rate, ΔG, in the crack tip within each cyclic loading and can be written as:

$$\partial \mathrm{a}/\partial \mathrm{N}=\mathrm{C}(\Delta \mathrm{G}⁄\mathrm{Gc})\mathrm{ }^\mathrm{m}$$

(1)

where G_C is the critical fracture energy, C and m are fitting parameters determined experimentally.

When the specimen is under mixed mode loading, i.e. Mode I and Mode II, the parameters C and m can be determined depending on linear or parabolic equations based on pure mode cases as in the following expressions ¹²

Linear rule

$$\mathrm{C}=\frac{{\mathrm{G}}_{\mathrm{I}}}{{\mathrm{G}}_{\mathrm{I}}+{\mathrm{G}}_{\mathrm{II}}}{\mathrm{C}}_{\mathrm{I}}+\frac{{\mathrm{G}}_{\mathrm{II}}}{{\mathrm{G}}_{\mathrm{I}}+{\mathrm{G}}_{\mathrm{II}}}{\mathrm{C}}_{\mathrm{II}}$$

(2)

$$\mathrm{m}=\frac{{\mathrm{G}}_{\mathrm{I}}}{{\mathrm{G}}_{\mathrm{I}}+{\mathrm{G}}_{\mathrm{II}}}{\mathrm{m}}_{\mathrm{I}}+\frac{{\mathrm{G}}_{\mathrm{II}}}{{\mathrm{G}}_{\mathrm{I}}+{\mathrm{G}}_{\mathrm{II}}}{\mathrm{m}}_{\mathrm{II}}$$

(3)

Parabolic rule

$$\log \mathrm{C}=\log {\mathrm{C}}_{\mathrm{I}}+\left(\frac{{\mathrm{G}}_{\mathrm{II}}}{{\mathrm{G}}_{\mathrm{T}}}\right)\log {\mathrm{C}}_{\mathrm{m}}+{\left(\frac{{\mathrm{G}}_{\mathrm{II}}}{{\mathrm{G}}_{\mathrm{T}}}\right)}^2\log \frac{{\mathrm{C}}_{\mathrm{II}}}{{\mathrm{C}}_{\mathrm{m}}{\mathrm{C}}_{\mathrm{I}}}$$

(4)

$$\mathrm{m}={\mathrm{m}}_{\mathrm{I}}+{\mathrm{m}}_{\mathrm{m}}\left(\frac{{\mathrm{G}}_{\mathrm{II}}}{{\mathrm{G}}_{\mathrm{I}}+{\mathrm{G}}_{\mathrm{II}}}\right)+({\mathrm{m}}_{\mathrm{II}}-{\mathrm{m}}_{\mathrm{I}}-{\mathrm{m}}_{\mathrm{m}}){\left(\frac{{\mathrm{G}}_{\mathrm{II}}}{{\mathrm{G}}_{\mathrm{I}}+{\mathrm{G}}_{\mathrm{II}}}\right)}^2$$

(5)

where the subscripts I and II indicates Mode I, Mode II, and mode ratio. m_m is a mixed mode parameter to be experimentally determined.

Fatigue delamination does not propagate if the strain energy rate does not reach the threshold value G_th. In this region no measurable fatigue crack growth occurs and it is called subcritical growth region, as shown in Figure 2 (Region I). The delamination grows linearly when the energy release rate during the cyclic loading exceeds the threshold value. This is known as the region of stable growth (Region II) when plotted on a log–log scale, as shown in Figure 2. The nonlinear or unstable growth (Region III) occurs if the strain energy release rate is near to the static fracture energy G_c, see Figure 2. ¹³

Figure 2.

Typical fatigue delamination growth curve. ¹³

The change in strain energy rate during fatigue loading can be calculated from the instantaneous maximum strain energy and the load ratio R as ¹⁷

$$\Delta \mathrm{G}=\left(1-{\mathrm{R}}^2\right){\mathrm{G}}_{\max }$$

(6)

where R is the load ratio and is defined by the user input in the numerical model. G_max is calculated from the traction–separation curve of the cohesive zone model (CZM) for each fatigue cycle.

The interface element is easily simulated in a numerical model based on traction-separation responses, the fracture process zone ahead of a crack tip can be modelled using a CZM. Many different types of traction-separation relations have been investigated in the literature e.g. polynomial law, ¹⁸ trapezoidal model, ¹⁹ perfectly plastic relationship, ²⁰ exponential law, ²¹ and bilinear cohesive law. ²² Among these the bilinear law gives good convergence of the finite element model and numerically is the simplest formulation.^13,23

It is worthwhile to mention that the static damage evolution should be included in modelling the delamination growth under fatigue conditions. The total damage parameter can be calculated during cyclic loading as:

$${\mathrm{D}}_{\mathrm{tot}}=\mathrm{D}+{\mathrm{D}}_{\mathrm{f}}$$

(7)

where D is static damage evolution and ${\mathrm{D}}_{\mathrm{f}}$ is fatigue damage parameter.

The total damage ${\mathrm{D}}_{\mathrm{tot}}$ is used to update the cohesive traction at each element after each time-step. Complete failure of an element takes place once the total damage reaches unity and ²⁴

$$\tau =\left(1-D\right)k\delta \mathrm{ }\mathit{static}\mathrm{ }\mathit{loading}\mathrm{ }$$

(8a)

$$\tau =\left(1-D_{\mathit{tot}}\right)k\delta \mathrm{ }\mathit{fatigue}\mathrm{ }\mathit{loading}$$

(8b)

where $\mathrm{ }\mathrm{\tau }$ is cohesive traction, k is interface stiffness, and δ is relative displacement

Energy release rate at the crack tip can be calculated by integrating the traction-separation curve under cyclic loading, unless complete failure occurs. The fatigue loading causes stiffness degradation in the interface element, this leads to increase separation as total damage evolution accumulates. The highest strain energy release rate will be in those elements that lie in the process zone near to the numerical crack tip. The response of traction separation will follow a non-vertical path which matches the bi-linear law for static loading according to Turon et al. ⁹ and Harper and Hallett. ²⁴

Developed fatigue damage degradation

The quadratic nominal stress criterion to predict the onset of delamination under mixed mode loading is one of the most frequently adopted failure criteria. ²⁵ This criterion has been successfully utilized by many researchers e.g. Hameed et al., ²⁶ Zou and Hameed, ²⁷ Belnoue et al. ¹⁷ and Kawashita and Hallett, ¹² and is written as;

$${F_i=\mathrm{ }\left(\frac{〈{\tau }_n〉}{{\tau }_{nc}}\right)}^2+{\left(\frac{{\tau }_s}{{\tau }_{sc}}\right)}^2+{\left(\frac{{\tau }_t}{{\tau }_{tc}}\right)}^2$$

(9)

where ${\tau }_n,\mathrm{ }{\tau }_s,\mathrm{ }{\tau }_t$ are the normal, shear and tear stresses respectively and where the operator $〈x〉\mathrm{ }$ is defined as;

$$〈x〉=\left\{\begin{array}{c}x, x>0\\ 0, x\le 0\end{array}\right.$$

Damage growth under pure-mode loading is predicted using the energy release rate during loading and unloading, and fracture toughness:

$$F_p=\frac{G}{G_c}$$

(10)

The static damage evaluation is derived from damage surface concept and is written as:

Damage initiation criterion ( $F_i$ ) + Damage propagation criterion ( $F_p$ ) = 1

From this condition, the static damage evaluation can be expressed as

$$D=1-\sqrt{\frac{1-F_p}{\left[\frac{k_1^2{\delta }_1^2}{{\tau }_{1c}^2}+\frac{k_2^2{\delta }_2^2}{{\tau }_{2c}^2}+\frac{k_3^2{\delta }_3^2}{{\tau }_{3c}^2}\right]}}$$

(11)

The ratio between damaged ( $L_d$ ) and undamaged( $L_u$ ) parts was assumed by Turon et al. ⁹ to equal the ratio of the dissipated energy over the representative interface part. Therefore, the damage growth criterion for pure mode can be written as:

$$F_p=\frac{G}{G_c}=\frac{L_d}{L_u}$$

(12)

Where $L_u$ is mesh size and $L_d$ is the length of damaged elements.

Substitute equation (12) into the static damage evaluation equation (11), to obtain the relation between damage evaluation and damaged length ratio;

$$D=1-\sqrt{\frac{1-\frac{L_d}{L_u}}{\left[\frac{k_1^2{\delta }_1^2}{{\tau }_{1c}^2}+\frac{k_2^2{\delta }_2^2}{{\tau }_{2c}^2}+\frac{k_3^2{\delta }_3^2}{{\tau }_{3c}^2}\right]}}$$

(13)

By differentiating the above equation with respect to the damaged length, we obtain a novel term $\frac{\partial D}{{\partial L}_d}$ , not found elsewhere in the literature and which can be written as:

$$\frac{\partial D}{{\partial L}_d}=\frac{1}{2}\mathrm{*}\frac{1}{L_u}\mathrm{*}\frac{1}{\sqrt{\frac{1-\frac{L_d}{L_u}}{\left[\frac{k_1^2{\delta }_1^2}{{\tau }_{1c}^2}+\frac{k_2^2{\delta }_2^2}{{\tau }_{2c}^2}+\frac{k_3^2{\delta }_3^2}{{\tau }_{3c}^2}\right]}}}\mathrm{*}\frac{1}{\left[\frac{k_1^2{\delta }_1^2}{{\tau }_{1c}^2}+\frac{k_2^2{\delta }_2^2}{{\tau }_{2c}^2}+\frac{k_3^2{\delta }_3^2}{{\tau }_{3c}^2}\right]}$$

(14)

The evaluation of the rate of fatigue damage for each cyclic loading ( $\partial D/\partial N$ ) on the traction–separation curve is a function of the novel term ( $\partial D/{\partial L}_d$ ) and the damage growth rate ( ${\partial L}_d/\partial N$ ) over the process zone. The expression for the required rate can be written as:

$$\frac{\partial D}{\partial N}=\frac{\partial D}{{\partial L}_d}\mathrm{*}\frac{{\partial L}_d}{\partial N}$$

(15)

So,

$$\frac{\partial D}{\partial N}=\frac{1}{2}\mathrm{*}\frac{1}{L_u}\mathrm{*}\frac{1}{\sqrt{\frac{1-\frac{L_d}{L_u}}{\left[\frac{k_1^2{\delta }_1^2}{{\tau }_{1c}^2}+\frac{k_2^2{\delta }_2^2}{{\tau }_{2c}^2}+\frac{k_3^2{\delta }_3^2}{{\tau }_{3c}^2}\right]}}}\mathrm{*}\frac{1}{\left[\frac{k_1^2{\delta }_1^2}{{\tau }_{1c}^2}+\frac{k_2^2{\delta }_2^2}{{\tau }_{2c}^2}+\frac{k_3^2{\delta }_3^2}{{\tau }_{3c}^2}\right]}\mathrm{*}\frac{{\partial L}_d}{\partial N}$$

(16)

The damage accumulates ahead of the crack tip in the fracture process zone to give additional delamination length. This new length of delamination is calculated based on the damage fatigue crack lengths, $L_d$ , of all elements which lie in the cohesive zone area. It is written as a summation function in the direction of delamination growth as:

$$a=\sum _{e\in L_{CZ}}{L}_d^e$$

(17a)

where $L_{CZ}$ is the length of the cohesive zone and calculated as ⁹

$$L_{CZ}=\frac{9\pi }{32}\frac{G_c}{{\left({\tau }_c\right)}^2}\mathrm{ }E$$

(17b)

E is the Young’s modulus of the material.

The equation (17a) is derived with respect to the number of cycles to obtain the crack growth rate, $\partial a/\partial N$ , as function of ( ${\partial L}_d/\partial N)$ .

$$\frac{\partial a}{\partial N}=\sum _{e\in L_{CZ}}\frac{{\partial L}_d^e}{\partial N}$$

(18)

The incremental rate of delamination is assumed to be equal for all elements within the cohesive zone area. Therefore, the summation sign can be rewritten in terms of the number of elements ( ${EN}_{CZ})$ in the process zone, $L_{C\mathrm{Z}}$ . ²⁸ The fatigue delamination growth rate can then be written as:

$$\frac{\partial a}{\partial N}=\mathrm{ }{EN}_{CZ}\mathrm{*}\frac{{\partial L}_d}{\partial N}\mathrm{ }$$

(19)

where

$${EN}_{CZ}=\mathrm{ }\mathit{length}\mathrm{ }of\mathrm{ }\mathit{cohesive}\mathrm{ }\mathit{zone}\mathrm{ }\left(L_{CZ}\right)/\mathit{element}\mathrm{ }\mathit{length}\mathrm{ }\left(L_u\right).$$

By substituting equation (19) into equation (16), the rate of fatigue damage ( $\partial D/\partial N$ ) is related to Paris’ law as:

$$\frac{\partial D}{\partial N}=\frac{1}{2}\mathrm{*}\frac{1}{L_u}\mathrm{*}\frac{1}{\sqrt{\frac{1-\frac{L_d}{L_u}}{\left[\frac{k_1^2{\delta }_1^2}{{\tau }_{1c}^2}+\frac{k_2^2{\delta }_2^2}{{\tau }_{2c}^2}+\frac{k_3^2{\delta }_3^2}{{\tau }_{3c}^2}\right]}}}\mathrm{*}\frac{1}{\left[\frac{k_1^2{\delta }_1^2}{{\tau }_{1c}^2}+\frac{k_2^2{\delta }_2^2}{{\tau }_{2c}^2}+\frac{k_3^2{\delta }_3^2}{{\tau }_{3c}^2}\right]}\mathrm{*}\frac{1}{{EN}_{CZ}}\mathrm{*}C{\left(\frac{\Delta G}{G_c}\right)}^m$$

(20)

In the UMAT subroutine, the fatigue damage degradation is updated in the model for each time-step as;

$$D_f^{\mathit{new}}=D_f^{\mathit{old}}+\Delta N\frac{\partial D}{\partial N}$$

(21)

$\Delta N$ is calculated as follows:

$$\Delta N=\mathit{frequency}\mathrm{*}\mathit{step}\mathrm{ }\mathit{time}$$

(22)

Double cantilever beam (DCB) model

Double cantilever beam (DCB) specimens as reported by Landry and LaPlante ¹³ were adopted for this numerical simulation. The length and width of the specimen was 150 mm and 20 mm respectively with pre-crack length of 50 mm. The specimens have two unidirectional 16 ply-thick (3 mm thick) layers made from carbon/epoxy with material properties shown in Table 1. All simulations in this study were subjected to Mode I static and fatigue loading, and carried out using commercial finite element software ABAQUS. A two dimensional finite element model of DCB was implemented to reduce computational effort and save computer time. The layers were meshed using 4-node plane strain elements, whereas the cohesive elements were meshed with element type COH2D4 and placed between layers where the crack was expected to propagate.

Table 1.

Material properties of the carbon/epoxy layers used. ¹³

E11 (GPa)	E22 (GPa)	E33 (GPa)	$v_{12}$	$G_{12}$ (GPa)	$G_{IC}$ (N/m)
155	10.5	10.5	0.32	4.83	422

First, the static load was numerically applied on the DCB by specifying an opening displacement of each cantilever arm. The numerical results for energy as a function of displacement have been validated with experimental data reported by Landry and LaPlante. ¹³ The response of strain energy with each initial displacement is in good agreement with experimental results, as shown in Figure 3.

Figure 3.

Comparison between developed model and experimental results of fracture toughness.

The experimental procedure was performed by Landry and LaPlante ¹³ as follows; the specimens were cracked approximately 5 mm and then subsequently subjected to the quasi-static load. When opening displacement reached to 9.5 mm, the static load was stopped and then, sinusoidal displacement cycles with an amplitude of 8 mm were performed.

The novel approach to fatigue damage degradation was implement in ABAQUS via the UMAT subroutine. It was necessary to define the cohesive zone length in the model which was estimated based on traction values, fracture toughness and Young’s modulus of the cohesive material. Fatigue delamination growth is determined by relating the damage accumulation of cohesive elements to the Paris growth law, equation (20). The damage accumulates with each cyclic loading based on equation (21) to produce a new crack length. Figure 4 illustrates damage accumulation for the DCB under different numbers of cycles. It is clear that the delamination length increases when the number of cycles increases from 250 to 3500. The difference between the experimental results and developed model is determined as;

$$\mathit{difference}=\frac{\mathit{experimental}\mathrm{ }\mathit{value}-\mathit{simulation}\mathrm{ }\mathit{value}}{\mathit{experimental}\mathrm{ }\mathit{value}}\mathrm{ }\mathrm{\%}$$

(23)

Figure 4.

Delamination growth in DCB under different cyclic loading.

Figure 5.

Comparison between developed model predictions and experimental results for delamination growth.

For example, the difference at 250 cycle is calculated as;

$$\mathit{difference}=\frac{61-58}{61}\mathrm{ }\mathrm{\%}\cong 5\mathrm{\%}$$

In the same way, the difference at 3500 cycles can be determined and it is about 3% as shown in Figure 5.

The fatigue crack growth rate using developed can be represented by Paris plot as shown in Figure 6. The simulation results of mixed mode were validated against experimental data available in literature. ²⁹

Figure 6.

Comparison of simulated fatigue crack growth rates and experimental data available in literature. ²⁹

Single lap joint (SLJ) model

Adhesively bonded joints have many advantages e.g. lower weight, and better fatigue performance compared with classical mechanically fastened joints. Therefore, they are used in many structural applications in industries such as automotive and aerospace. The increase in demand for adhesively bonded joints has been accompanied by investigating and developing numerical models to predict the response of the structure, and to describe the behaviour of damage in the adhesive joints. ³⁰ The cohesive zone model is a very effective way to model adhesively bonded joints using finite element analysis

Khoramishad et al. ³¹ carried out static and fatigue loading of single lap joints (SLJ). The specimen consisted of two aluminium 2024-T3 layers bonded with FM 73 M OST toughened epoxy film adhesive. The substrates were pre-treated prior to bonding. This pre-treatment consisted of a chromic acid etch (CAE) and phosphoric acid anodising (PAA) followed by the application of BR 127 corrosion inhibiting primer, to maximise environmental resistance and bonding durability. The joints were cured at 120 °C, under 0.28 MPa pressure for 60 min. The dimensional details of the SLJ are shown in Figure 7. The overlap length, the width and the thickness of the bond line were 30, 12.5 and 0.2 mm, respectively.

Figure 7.

Single lap joint geometries.

The static tests were carried out by Khoramishad et al. ³¹ on the single lap joints (SLJ) to find cohesive zone properties e.g. fracture energy. The displacement control was applied to the specimen at a rate of 0.1 mm/min. Average static strengths of 10.0 kN was observed for the SLJ. After that, the single lap joints (SLJ) specimen was tested under fatigue load condition at 5 Hz with a load ratio of 0.1.

The interface element length used over the cohesive zone depends on mesh size in the finite element model. It effects the sensitivity of the fatigue response according to Naghipour et al. ³² Vassilopoulos, ⁵ in previous studies investigated the effect of mesh element size on convergence of the numerical model and found that using a mesh size less than 1 mm for the interface element provided better solution convergence. Therefore, in the present simulation the element size adopted for the adhesive was 0.2 mm.

Fatigue modelling was performed at maximum fatigue load of 40% of the static final failure load which was measured experimentally at about 10 kN. When the specimen was loaded statically at 4 kN, as shown Figure 8, it can be seen that, initially, the damage parameter is less than unity for each of the three elements. When the fatigue loading was activated, the first and second elements needed 200 cycles for the rate of increase of damage to begin to accelerate rapidly whereas in the third element the rate of increase of damage began to accelerate after about 500 cycles. This was because the value of the damage parameter for the third element at the end of the static stage was less than for the first and second elements.

Figure 8.

Damage propagation of the first, second and third element.

The delamination growth was simulated based on the model developed here. It should be observed that the damage grew symmetrically on both side of the adhesively bonded joint as shown in Figure 9. The variation of the damage parameter under different numbers of cycles, as predicted by the present model is also illustrated in Figure 9. The elements having state variable (SDV3) equal to unity in the damage contours represent complete failure of the cohesive zone, signifying no fatigue resistance. Based on model proposed here the crack needed 2000 cycles to grow to 0.30 mm and propagated to 0.45 mm when the number of cycles reached 4000. The measured delamination lengths by Khoramishad et al. ³¹ were 0.285 mm and 0.43 mm for 2000 and 4000 cycles respectively. Thus, the difference between numerical model and test was 5.2% at 2000 cycles and 4.6% at 4000 cycles.

Figure 9.

Delamination growth under different cyclic loading (a) developed model (b) experimental data by Khoramishad et al. ³¹

Conclusions

The objective of this paper was to present a novel fatigue delamination accumulation model. The novel approach was tested numerically under two different load conditions; Mode I loading which was performed on a DCB specimen, and Mode II loading which was carried out on an adhesively bonded joint. In addition, the simulation was performed of the mixed-mode model and compared with experimental data available in the literature ²⁹ as shown in Figure 6. The simulation findings were compared against experimental results available in literature, the maximum difference between them was 5.0% for the DCB specimen and 5.2% for the SLJ specimen. This is taken as confirmation that the proposed model converges with experimental data when tested using numerical simulation.