Engineering Fracture Mechanics 266 (2022) 108353

Contents lists available at ScienceDirect

Engineering Fracture Mechanics

journal homepage: www.elsevier.com/locate/engfracmech

Extrinsic cohesive zone modeling for interface crack growth:

Numerical and experimental studies
P.J. Saikia, N. Muthu *
Indian Institute of Technology Guwahati, Guwahati, Assam 781039, India

A R T I C L E I N F O A B S T R A C T

Keywords: The present work is concerned with the study of interface crack using the extrinsic cohesive zone
Cohesive zone modeling model (CZM) that nullifies the crack tip stress singularity. A semi-analytical procedure combining
Interaction integral the generalized crack opening displacement solution and the finite element (FE) result is shown to
Interface fracture
estimate the traction-separation law (TSL) parameters. The stress intensity factor (SIF), which
Stress intensity factor
symbolizes the strength of the singularity, is computed using the domain-independent interaction
integral that considers crack-face traction. The SIF obtained from the FE simulations is plotted
along the cohesive zone length to demonstrate the reduction in the strength of the crack tip
singularity and found in good agreement with the analytical solutions. Finally, the effectiveness of
the proposed procedure is validated for an interface crack propagation between isotropic parent
materials.

1. Introduction

In recent times, material interfaces have been part of numerous structural and engineering applications, ranging from adhesive
joining to laminated composites. Often, interfacial stresses cause debonding followed by crack propagation leading to failure, thus
weakening the structure. With increasing demands on multifunctional needs in mechanical, aerospace, and biomedical applications,
understanding the mechanical behavior of such interfaces is very essential for optimizing the design. Therefore, accurate modeling of
the interfacial crack propagation phenomenon is vital to avoid material failure during service loading conditions.
The cohesive zone model (CZM) is a widely used technique for analyzing the interface fracture phenomenon within computational
fracture mechanics studies. It allows a smooth transition between the strength and energy-based fracture mechanics approach pro­
posed by Inglis [1] and Griffith [2], respectively. The CZM concept originated from Barenblatt [3] and Dugdale’s [4] strip yield model.
The main idea behind the CZM is to eliminate the stress singularity of the asymptotic solution inherent to the linear elastic fracture
mechanics (LEFM). The constitutive response of the fractured surfaces that form the material interface follows a traction-separation
law (TSL) characterized by the interface strength, stiffness, and fracture toughness. The damage in the cohesive zone will initiate
whenever the interfacial stresses reach the failure strength of the material. The damage subsequently evolves, resulting in crack
propagation when the expended energy becomes equal to the fracture energy or the final separation of the cohesive zone reaches a
critical separation value. Various TSLs were proposed to correlate the experimental load vs. crack growth data of cracked structures.
Some of the popular TSLs include – the linear softening cohesive law by Camacho and Ortiz [5], the exponential cohesive law by
Needleman [6,7], and Xu and Needleman [8], the trapezoidal and the polynomial cohesive law by Tvergaard and Hutchinson [9]. The

* Corresponding author.
E-mail addresses: [email protected] (P.J. Saikia), [email protected] (N. Muthu).

https://doi.org/10.1016/j.engfracmech.2022.108353
Received 3 January 2022; Received in revised form 19 February 2022; Accepted 22 February 2022
Available online 3 March 2022
0013-7944/© 2022 Elsevier Ltd. All rights reserved.
P.J. Saikia and N. Muthu Engineering Fracture Mechanics 266 (2022) 108353

Nomenclature

2a Original crack length

2b Actual crack length considering the cohesive zone
An , Bn Coefficients of complex Eigenvalue expansion
c n , dn Coefficients used to define near tip stress field
E Young’s modulus
GF Cohesive fracture energy
KI , KII Pure mode I and mode II stress intensity factor, respectively
Kapplied , Kcohesive Stress intensity factor generated by the far-field loading and the presence of cohesive zone, respectively
l, b, t Geometrical dimensions of the specimen
k Kolosov constant
nx , ny Direction cosines along x and y-direction, respectively
N Exponent defining the cohesive law profile
SE Brittleness number at the interface
Tx , Ty Traction field applied at the periphery of the K-dominated region
u, v Near tip displacement field in the cohesive zone along x and y-direction, respectively
z Complex variable used to represent near tip fields
δn Normal separation along the cohesive zone
δ̂n Normalized normal separation
δs In-plane shear separation along the cohesive zone
δ̂s Normalized shear separation
δcn , δct Critical normal and shear separation along the cohesive zone, respectively
εaux
jk , εjk Auxiliary and total strain, respectively
η(x) Arbitrary location along the crack axis
μ Shear modulus
ν Poisson’s ratio
ρ Cohesive zone length
σn Normal cohesive traction
σ xx Near tip stress field along the cohesive zone in x- direction
σ yy Near tip stress field along the cohesive zone in y-direction
σ∞ Far-field normal stress
σ max , τmax Peak normal and shear traction along the cohesive zone
τs Shear cohesive traction
τxy Near tip stress field along xy-plane
τ∞ Far-field shear stress
φ(z), χ (z) Complex analytical functions used to represent the near tip field

trapezoidal cohesive law can be used for crack propagation in elastic-plastic solids and recently used to evaluate the force on the
bridging platelet and investigate the toughening mechanism of nacre [10]. These TSLs are popularly known as potential-based intrinsic
CZMs, where the cohesive traction is zero at its initial separation displacement. The constitutive behavior of the fracture surface is
defined by a potential function that represents the energy distribution along the fracture surface. However, none of those above-
mentioned TSLs can nullify the stress singularity, which is one of the main objectives of the CZM. The existing literature showed
that a TSL with non-zero traction at the onset of the separation could only remove the stress singularity [11]; such CZMs are cate­
gorized as extrinsic cohesive zone models, which is the focus of the present work. In an extrinsic CZMs, the cohesive traction reaches a
finite value at zero separation displacement, decreases as the cohesive zone widens up, and reaches zero as the separation reaches the
critical separation value. The simplest extrinsic CZMs to nullify the stress singularity are the constant traction and linear softening
cohesive laws. In the case of extrinsic CZMs, the cohesive elements are adaptively inserted while carrying out the simulation [12],
unlike the intrinsic CZMs, where the cohesive elements are initially introduced along the potential failure path [12–14].
The stress intensity factor (SIF) parameter, pertinent to the LEFM, characterizes the stress singularity associated with a crack tip.
There are various computational techniques [15–21] to compute the SIF within the framework of computational fracture mechanics.
The determination of the SIF for a crack embedded with a cohesive zone was also emphasized in some of the earlier established work to
estimate the effectiveness of the cohesive model [22,23]. It was observed that the cohesive zone’s characteristic length and cohesive
traction play a vital role in nullifying the stress singularity. While for a constant TSL, closed-form solutions are available to evaluate the
cohesive zone’s size [24], there exists no approach to determine the cohesive zone parameters for other TSLs, including linear and non-
linear softening laws; only complicated expression containing integrals are available for nullifying the strength of singularity at the tip
of the cohesive zone [11,25]. This has motivated the author for the present work, and accordingly, it aims to address the following
questions explicitly. For any generalized TSLs, what is the influence of the characteristic length of the cohesive zone on the SIF? At

2
P.J. Saikia and N. Muthu Engineering Fracture Mechanics 266 (2022) 108353

what will be the values of the cohesive zone parameters for which the SIF drops to zero?
Rice and Sih [26] gave the expression in the integral form to determine the SIF due to cohesive traction. To evaluate the SIFs, the
expressions of the cohesive stress variation in the deformed configuration are needed. But, the analytical expression of the normal and
shear tractions along the cohesive zone for linear or higher-order TSLs does not exist due to the complex nature of the problem.
Previously, Jin and Sun [11] and Wang [25] proposed a computationally involved iterative numerical methodology for evaluating the
SIF’s in the presence of the cohesive zone. A semi-analytical solution of the normal and shear tractions for any general extrinsic TSL
was obtained in this present work. A general expression of the asymptotic stress and displacement field for mode-I cohesive crack was
derived by Xiao and Karihaloo [27], which we compare with the FE solution to determine the actual solution. This stress and
displacement field solution corresponding to a particular TSL was used to evaluate the size of the cohesive zone for nullifying the SIF.
This approach was extended to mode II cohesive crack studies. To ascertain the null singularity, the SIF has been computed using the
domain-independent interaction integral procedure [28–30]. The crack-face cohesive tractions were taken into consideration while
evaluating the SIFs by interaction integral technique. The comparison was made between the normal contour integral devoid of crack
face traction and the domain integral to demonstrate the domain independency of the technique.
Finally, the present work synthesizes the proposed semi-analytical procedure and mode-I experimental interface crack propagation
studies. Concerning interface crack propagation experiments, Kinloch and Shaw [31] earlier used a double cantilever beam epoxy
glued steel sandwich specimen to determine the interface properties where the crack is allowed to propagate in the pre-cracked ad­
hesive layer. Banks-sills et al. [32–34] proposed a procedure for evaluating the interface fracture properties of isotropic and anisotropic
interfaces and validated it by the experimentation of glass/epoxy and ceramic clay material interfaces. As per the author’s knowledge,
there is no direct application of extrinsic CZMs involving null SIFs in experimental studies involving interface crack propagation.
Therefore, this work conducted a pure mode I fracture study on a non-dimensional glued disc-shaped CT (compact tension) specimen.
The load-displacement behavior was validated by selecting a linear-softening TSL with suitable cohesive zone parameters for the
numerical studies. By comparing the load-displacement behavior of the numerical studies with the experimental results, it was evident
that the proposed procedure for extrinsic CZMs can be used for modeling interface crack propagation problems under varieties of TSLs
with non-zero traction at the onset of opening displacement.

2. Cohesive zone model for the removal of crack tip stress singularity

Consider a plate with a crack of length 2a at the center, as shown in Fig. 1. The plate is subjected to a far-field loading consisting of
tensile stress σ ∞ and shear stressτ∞ . Due to the external loading, a cohesive zone of length ρ = b − a developed ahead of the crack tip.
Assuming ρ ≪ a, we can consider a ≈ a + ρ. The stress intensity factor (SIF) [24] generated due to the applied load is given by.
√̅̅̅̅̅
Kapplied = (σ∞ + iτ∞ ) πa (1)

By considering Green’s function approach, Rice and Sih [26] evaluated the SIF generated due to the cohesive traction at the tip of
the cohesive zone.
√̅̅̅ ∫ ρ
2
(2)
1
Kcohesive = − [σn (η) + iτs (η)]η− 2 dη
π 0

where σn (η) and τs (η) are the normal and shear cohesive traction at any arbitrary location η along the crack axis.

Fig. 1. Infinite central crack plate with the K-dominated region used for the analysis.

3
P.J. Saikia and N. Muthu Engineering Fracture Mechanics 266 (2022) 108353

For the removal of stress singularity, we can use

Kapplied + Kcohesive = 0 (3)
√̅̅̅ ∫ ρ
√̅̅̅̅̅ 2
(4)
1
(σ∞ + iτ∞ ) πa = [σ n (η) + iτs (η)]η− 2 dη
π 0

For a cohesive zone with constant normal and shear traction σn (η) = σ n and τs (η) = τs , Eq. (4) leads to
√̅̅̅ ∫ρ
√̅̅̅̅̅ 2
(5)
1
(σ∞ + iτ∞ ) πa = (σ n + iτs ) η− 2 d η
π 0

√̅̅̅
√̅̅̅̅̅ 2 √̅̅̅
(σ∞ + iτ∞ ) πa = (σ n + iτs )2 ρ (6)
π
Assuming a pure mode I or mode II case, separating the real and imaginary part of Eq. (6), we get
√̅̅̅
√̅̅̅̅̅ 2 √̅̅̅
σ∞ πa = σn 2 ρ
π
√̅̅̅ (7)
√̅̅̅̅̅ 2 √̅̅̅
τ∞ π a = τs 2 ρ
π
Therefore, for a given constant normal and shear cohesive traction σ n , τs , and a prescribed loading σ∞ and τ∞ , the cohesive zone
length ρ can be evaluated from Eq. (7).
However, for any generalized TSL, the cohesive stress distribution σ n (η) or τs (η) mentioned in Eq. (4) along the cohesive zone is not
readily available to determine the cohesive zone parameters for the removal of stress singularity. Due to the complexity of the problem,
a semi-analytical approach is adopted to nullify the stress singularity.

3. Analytical formulations of the cohesive zone model

The asymptotic pure mode II stress and displacement field (assuming no friction between the fractured surfaces) ahead of the
cohesive zone tip was obtained by adopting the similar analytical procedure of Xiao and Karihaloo [27] for their pure mode I analysis.
The above-mentioned formulation is an extension of Muskhelishvili’s problem of planer elasticity. According to Muskhelishvili [35],
the convenient representation of stress and displacement in cartesian coordinate is

σ xx + σ yy = 2[φ’(z) + φ’(z)]
σ yy − σ xx + 2iτxy = 2[zφ˝(z) + χ ˝(z) ] (8)
2μ(u + iv) = kφ(z) − zφ’(z) − χ ’(z)

k and μ are the Kolosov constant and the shear modulus, respectively. For elasticity problems, k = 3 − 4ν and k = 3−1+ν for plane strain
ν

and stress conditions, respectively. E and ν are Young’s modulus and the Poisson’s ratio respectively. z is a complex variable in polar
coordinate with z = reiθ , as shown in Fig. 2. The prime and overbar are used for the differentiation and the complex conjugate,
respectively. Sih and Liebowitz [36] expressed the analytical functions φ(z) and χ (z) as a series of complex Goursat functions:
∑ ∑
φ(z) = An zλn = An rλn eiλn θ
(9)
n=0 n=0
∑ ∑
χ (z) = Bn zλn +1 = Bn rλn +1 ei(λn +1)θ
n=0 n=0

Fig. 2. Asymptotic near tip field in the presence of a cohesive zone.

4
P.J. Saikia and N. Muthu Engineering Fracture Mechanics 266 (2022) 108353

An = a1n +ia2n and Bn = b1n +ib2n are the coefficients of complex eigenvalue expansion. As the crack propagation occurs, the
cohesive zone opens up, and the fracture process zone (FPZ) is developed ahead of the crack tip. The end coordinate of FPZ is defined
mainly by two critical parameters; a physical tip where cohesive stress vanishes and a mathematical tip with zero crack opening and
maximum traction, as shown in Fig. 2. Xiao and Karihaloo [27] obtained the asymptomatic stress and displacement field near the
cohesive crack tip by substituting the complex function Eq. (9) into Eq. (8). The generalized expression for stress and displacement in
the presence of the cohesive zone was provided in the appendix section. A comparison between the work of Xiao and Karihaloo [27]
and the obtained analytical formulations was summarized in Table 1.

4. Solution methodology

The K-dominated region covering the cohesive zone tip at its center, as shown in Fig. 1, was considered for the analysis. This
simulation approach is identical to Zavattieri et al. [37] for their crack propagation along with the rough interfaces. A singularity-
dominated circular region was considered with known Linear elastic fracture mechanics (LEFM) pure mode I and mode II stresses
were applied to the outer periphery, as shown in Fig. 3(a).

4.1. Model geometry and boundary conditions

A central cohesive crack was assumed, and its opening was observed independently under a pure mode I and mode II far-field
loading. The mesh generation was shown in Fig. 3(b), with the coordinate system at the center of the circular domain.
The commercial FE package ABAQUS/standard interface was used to carry out the numerical analysis. The model geometry was of
2D shell features. The parent elements were represented by a four-node bilinear plane stress quadrilateral hourglass control with four-
noded 2D cohesive elements were placed ahead of the physical tip. The geometrical dimensions of the problem were
r = 100 mm, a = r − ρ, t=1

where t is the plane stress thickness, a and r are the original crack length and the radius of the circular domain, respectively. ρ rep­
resents the cohesive zone length used for the simulation. The material properties of the parent material were.
√̅̅̅̅̅̅̅̅
E = 70GPa, ν = 0.3, K = KI = KII = 100 MPa mm

where E is Young’s modulus, ν is the Poisson’s ratio.KI and KII are the prescribed SIF for pure mode I and II, respectively. The traction
field applied at the periphery of the circular domain can be expressed as
[ ] [ ][ ]
Tx σ xx τxy nx
= (10)
Ty τxy σ yy ny

where for pure mode I,

( )[ ( ) ( )]
KI θ θ 3θ
σ xx = √̅̅̅̅̅̅̅ cos 1 − sin sin
2π r 2 2 2
( )[ ( ) ( )]
KI θ θ 3θ
σ yy = √̅̅̅̅̅̅̅ cos 1 + sin sin (11)
2πr 2 2 2
[ ( ) ( ) ( )]
KI θ θ 3θ
τxy = √̅̅̅̅̅̅̅ sin cos cos
2π r 2 2 2

For pure mode II,

Table 1
Comparative study between the present formulation with Xiao and Karihaloo [27].
Parameters Xiao and Karihaloo [27] (mode I) Present analysis (mode II)
⃒ ⃒ ⃒ ⃒
Boundary conditions =0
σyy ⃒θ=π = σyy ⃒θ=− π ∕ σyy ⃒θ=π = σyy ⃒θ=− π = 0
⃒ ⃒ ⃒ ⃒
τxy ⃒θ=π = τyy ⃒θ=− π = 0 =0
τxy ⃒θ=π = τyy ⃒θ=− π ∕
⃒ ⃒
τxy ⃒ θ=0 = 0, v|θ=0 = 0 σyy = σxx ⃒ θ=0 = 0, u|θ=0 = 0
Asymptotic displacement field 2n + 3 2n + 3
) )
∑ r 2 ( 2n + 3 ∑ r 2 ( 2n + 3
δn = n=0 a1n (k +1)sin π δs = n=0 b2n (k +1)sin π
μ 2 μ 2
Effective normalized displacement 2n + 3 2n + 3
∑ ∑
δn = n=0 cn r μ
̂ δ̂s = n=0 dn r μ
( ) ( )
(k + 1) 2n + 3 (k + 1) 2n + 3
cn = a sin π dn = b sin π
μδcn 1n 2 μδcs 2n 2

5
P.J. Saikia and N. Muthu Engineering Fracture Mechanics 266 (2022) 108353

Fig. 3a. Prescribed boundary condition.

Fig. 3b. Finite element mesh of the geometry.

( )[ ( ) ( )]
KII θ θ 3θ
σ xx = − √̅̅̅̅̅̅̅ sin 2 + cos cos
2π r 2 2 2
[ ( ) ( ) ( )]
KII θ θ 3θ
σ yy = √̅̅̅̅̅̅̅ sin cos cos (12)
2πr 2 2 2
( )[ ( ) ( )]
KII θ θ 3θ
τxy = √̅̅̅̅̅̅̅ cos 1 − sin sin
2π r 2 2 2

σxx , σyy , and τxy were near tip stress fields in the presence of the cohesive zone. nx = cosθ, ny = sinθ was the direction cosines along the x
and y-axis, respectively.
For pure mode I and mode II, displacement ahead of the cohesive tip along the crack axis was constrained in the y and x-direction,
respectively, as mentioned in Table 1. Two different regions were considered in the model geometry. One was the square region
surrounding the cohesive tip at its center, with the other one circular region covering the square domain. Each side of the square was
twice that of the value of cohesive zone length, ρ considered. A square domain (consisting of eight elements, excluding the cohesive
one) with the cohesive zone tip at its center was considered for evaluating the interaction integral procedure.

6
P.J. Saikia and N. Muthu Engineering Fracture Mechanics 266 (2022) 108353

4.2. Cohesive interface model

The cohesive behavior for negligible small thickness interface elements was taken care of by a user-defined subroutine (UMAT). A
linear-elastic constitutive equation was used to define the parent material behavior. The constitutive behavior of the cohesive was
defined as
For the pure mode I,
⎧ ( )N
⎪
⎨σ δn
max 1 − , δ < δcn
σ(δ) = δcn (13)
⎪
⎩
0, δ ≥ δcn

For the pure mode II,

⎧ ( )N
⎪
⎨τ δs
max 1 − , δ < δcs
(14)
c
τ(δ) = δs
⎪
⎩
0, δ ≥ δcs

where σmax and τmax are the peak cohesive normal and shear traction. N is an exponent defining the shape of the cohesive law. For N = 0
andN = 1, the above Eq. (13) and Eq. (14) results in constant traction and linear softening cohesive law, respectively.
The damage evolves as soon as the cohesive zone opens up for the pure mode I, followed by the element deletion once it reaches its
critical separation δcn . However, for pure mode II, negative separation generates anti-shear traction in the cohesive.

4.2.1. Constant traction cohesive law

Throughout the simulation, the cohesive stress was assumed to be constant and equal to σmax and τmax for pure mode I and mode II,
respectively. Once the critical separation is achieved, the cohesive energy approaches its critical fracture energy GF , and stress
drastically drops to zero.
The assumed material properties were

K2 N
σmax = τmax = 40 MPa, GF = = 0.1428 .
E mm
The cohesive zone length ρ is varied from 1 mm to 3 mm with an increment of 0.25 mm, and Kcohesive for a particular value of ρ can be
obtained directly from Eq. (7). Five contours of arbitrary shape and size were assumed in the ABAQUS, and their average value of
contour integral was presented. A comparison between the overall SIF obtained from contour integral, interaction-integral, and
analytical formulations were shown in Table 2 and 3 for pure mode I and mode II, respectively. The results obtained by the contour
integral became inaccurate and deviated more from the analytical results with larger values of cohesive zone length. However, the
interaction integral results were closely matched with the analytical results.

4.2.2. Linear-softening cohesive law

The cohesive stress was σmax and τmax for pure mode I and mode II, respectively, at the start of the simulation. The cohesive traction
decreases as the cohesive zone opens up and approaches zero at its final separation displacement. The constitutive behavior of the
cohesive can be represented as
For pure mode I,
⎧ ( )
⎨ σmax 1 − δn , δ < δc
⎪
(15)
c n
σ(δ) = δn
⎪
⎩
0, δ ≥ δcn

For pure mode II,

⎧ ( )
⎨ τmax 1 − δs ,
⎪
δ < δcs
(16)
c
τ(δ) = δs
⎪
⎩
0, δ ≥ δcs

The assumed material properties are.

K2 N
σmax = τmax = 51.2 MPa, GF = = 0.1428 .
E mm
As mentioned for the constant traction cohesive law, the cohesive zone length was assumed from ρ = 1.5 mm to 3 mm with an
identical step size of .25 mm.
We used the pure mode I and II asymptotic displacement field mentioned in Table 1 to approximate the stress field ahead of the
cohesive tip. After obtaining the displacement field numerically around the cohesive tip for a specified cohesive law, unknown co­
efficients were determined by comparing the numerically obtained field with the analytical field mentioned in Table 1. Coefficient up

7
P.J. Saikia and N. Muthu Engineering Fracture Mechanics 266 (2022) 108353

Table 2
Overall SIF obtained for constant TSL (pure mode I).
ρ (mm) Contour integral Interaction-integral Analytical formulation

1 41 38.55 36.15
1.25 23.80 30.83 28.61
1.5 20.45 23.39 21.80
1.75 17.35 16.28 15.53
2 14.48 9.888 9.706
2.25 11.71 3.756 4.23
2.5 9.132 − 2.279 − 0.9508
2.75 6.598 − 7.996 − 5.87
3 4.158 − 13.65 − 10.58

Table 3
Overall SIF obtained for constant TSL (pure mode II).
ρ (mm) Contour integral Interaction-integral Analytical formulation

1 43.56 36.05 36.15

1.25 22.85 30.14 28.61
1.5 20.20 23.44 21.80
1.75 17.80 17.48 15.53
2 15.41 11.66 9.706
2.25 13.18 6.210 4.23
2.5 10.98 1.036 − 0.9508
2.75 8.856 − 3.938 − 5.87
3 6.786 − 8.840 − 10.58

to their fifth-order was considered to get desired computational accuracy, as shown in Fig. 4. For a linear softening cohesive law, the
obtained value of the coefficients and their variations along with the cohesive zone length was shown in Fig. 5 and Fig. 6 for pure mode
I and mode II, respectively. All the coefficients were monotonically increased or decreased with the change in cohesive zone length ρ.
After obtaining the displacement field, the stress field in the cohesive crack for a particular cohesive law can be easily evaluated.
Table 4 and 5 summarized the overall SIF obtained by contour, M-integral, and semi-analytical formulations for mode I and mode II,
respectively. For both the constant and linear softening TSLs, the overall SIF’s were getting decreased with the increasing values of
cohesive zone length. Additionally, a constant cohesive law generated more SIF’s compared to a cohesive law with linear softening
behavior with increasing values of cohesive zone length. Figs. 7(a)-7(b) and Figs. 8(a)-8(b) compared the pure mode I and mode II
displacement and stress behavior for a linear softening TSL respectively for a specific value ofρ = 1.75mm. The agreement was
excellent, with a minimum error less than a reasonable limit. It is important to mention that the present methodology is only developed
for the crack in an interface embedded in isotropic media. Deriving a generalized expression of displacement in the presence of a
cohesive zone for a crack in a bimaterial interface is another major challenge and is still an open problem for the scientific community.
One of the essential aspects of the study was to estimate the SIF by the domain-independent interaction integral and contrast with
the contour integral in ABAQUS. In the case of interaction integral, three square domains with the mathematical tip at its center were
used to calculate the SIFs, as shown in Fig. 9. For contour integral, five regions were considered for evaluating the SIF. Abaqus
automatically locates the elements that make up each ring from the crack tip regions. [38]. The SIFs obtained using the domain-
independent interaction integral are consistent for all three domains. The main reason for this is attributed to the inclusion of
crack face traction terms in the integral. This is given by Shih et al. [39].
∫ [ ) ] ∫
1( ( )
IS = σ ij uaux aux
j,1 + σ ij uj,i − σ jk εaux aux
jk + σ jk εjk δ1,i q,i dA − tj uaux
j,1 qdC (17)
A 2 +
C +C −

σaux aux
ij ,uj and, σij , uj represent the stress and displacement vector due to the auxiliary field and actual field, respectively. tj is the
cohesive traction on the crack face C+ and C− mentioned in Fig. 10. uj,1 denotes a change in displacement with x1 coordinate. A is a
simply connected domain enclosed by the contourC.q is a smoothly increasing function enclosed by the contour C with zero value on C1
and unity inΓ. IS represents the interaction integral term. The detailed interaction integral formulation is provided in the appendix
section.
The in-built contour integral in ABAQUS lacks crack face traction term, and therefore, the SIF obtained is dependent on the contour.
Inclusion of crack face traction would have solved the issue, but mode separation would still be an issue for contour integral. Therefore,
the domain-dependent interaction integral was used for the subsequent analysis. The SIFs obtained by both the interaction integral and
contour integral for a particular cohesive zone parameter ρ = 1.75 mm and σ max = 40 MPa for a constant cohesive law are shown in
Fig. 11.

8
P.J. Saikia and N. Muthu Engineering Fracture Mechanics 266 (2022) 108353

Fig. 4. Polynomial order considered to estimate the displacement field (for pure mode I, ρ = 1.75 mm).

Fig. 5. Variation of the coefficients with the cohesive zone length (pure mode I).

Fig. 6. Variation of coefficients with the cohesive zone length (pure mode II).

9
P.J. Saikia and N. Muthu Engineering Fracture Mechanics 266 (2022) 108353

Table 4
Overall SIF obtained for linear-softening TSL (pure mode I).
ρ (mm) Contour integral Interaction integral Semi-analytical formulation

1.5 25.15 26.97 25.37

1.75 22.43 20.54 19.07
2 19.11 13.05 11.93
2.25 13.76 1.87 1.40
2.5 4.47 − 15.99 − 15.83
2.75 − 2.35 − 29.51 − 28.66
3 − 7.44 − 39.89 − 37.71

Table 5
Overall SIF obtained for linear-softening TSL (pure mode II).
ρ (mm) Contour integral Interaction-integral Semi-analytical formulation

1.5 25.37 28.65 27.61

1.75 23.95 24.49 22.70
2 22.29 20.05 17.74
2.25 20.89 15.82 13.30
2.5 17.05 7.057 4.773
2.75 11.33 − 5.318 − 8.061
3 − 0.2774 − 28.348 − 30.95

5. Experimental procedure for predicting crack growth

The material considered for the experimental studies was colorless and amorphous thermoplastic polymer PMMA (Polymethyl
methacrylate), popularly known as acrylic or Plexiglas. PMMA is a compatible structural polymer and can be joined with most
commercial adhesives to study the interface behavior. A commercially available 10mm thickness acrylic sheet with the brand name
BIGMALL (India) was considered to conduct the experimental investigation.

5.1. Strength and toughness measurement

As there is no standard method for evaluating the interfacial fracture toughness, a procedure identical to Alam et al. [40] for
measuring the deflection-penetration behavior of the PMMA interface is adopted. Before proceeding ahead, accurate material prop­
erties (Young’s modulus and Poisson’s ratio) were essential to perform interface crack propagation simulation; five specimens were
tested as per ASTM D638-14 test method is to evaluate the requisite elastic properties. A schematic diagram of a dog-bone-shaped
tensile specimen was shown in Fig. 12. All the specimens were machined in a CNC (computer numeric control) laser cutting ma­
chine to confirm good dimensional accuracy. The tests were conducted in a 250 kN BISS servo-hydraulic Universal Testing Machine
with a 1 mm/min strain rate. The obtained elastic properties after the experimentation are shown in Table 6.
A non-standard fracture test was conducted to evaluate the Mode-I fracture toughness at the interface. The schematic disc-shaped
CT specimen of thickness t = 10 mm used for the experimentation was shown in Fig. 13. A low strength, brittle Araldite resin (Epofine

Fig. 7a. Mode-I displacement profile of the crack.

10
P.J. Saikia and N. Muthu Engineering Fracture Mechanics 266 (2022) 108353

Fig. 7b. Mode-I stress profile of the crack.

Fig. 8a. Mode-II displacement profile of the crack.

Fig. 8b. Mode-II stress profile of the crack.

11
P.J. Saikia and N. Muthu Engineering Fracture Mechanics 266 (2022) 108353

First contour

Second contour

Third contour

Fig. 9. Domain for the interaction integral.

Fig. 10. Schematic of interaction integral under crack face traction C = C1 − Γ + C+ + C− .

Fig. 11. Domain independence for interaction-integral.

1564) is adequately mixed with an equal volume of hardener (FINEHARD-3486-2) and used to join the native PMMA with a minimal
thickness and avoid discontinuity in local properties. All the interfaces were properly cleaned by fine-grade sandpapers and cleaned
with acetone to remove the debris and impurities. Two semicircular pieces of acrylic sheet are glued by the above-mentioned adhesive
ahead of the notch tip from the center of the mid-plane. Two circular holes of 13mm diameter were machined on each specimen to
properly clamp the specimen in suitable fixtures. The toughness of the interface of all the five identical samples was tested in a 25kN

12
P.J. Saikia and N. Muthu Engineering Fracture Mechanics 266 (2022) 108353

29

R 76

57 19 115 246

Fig. 12. Schematic of the tensile specimen (all dimensions are in mm).

Table 6
Mechanical properties of acrylic.
Property Value

Young’s modulus (GPa) 2.88

Poisson’s ratio 0.38

Fig. 13. Schematic of non-dimensional CT specimen.

13
P.J. Saikia and N. Muthu Engineering Fracture Mechanics 266 (2022) 108353

INSTRON 8801 digitally controlled dynamic testing machine, as shown in Fig. 14. A quasi-static test is performed at a displacement
rate of 0.1mm/min to determine the overall load–deflection behavior of the interface accurately. The load resist by the specimen
increases linearly with the separation and drastically drops once crack propagation occurs. The instantaneous drop in the resistance
mimics the perfectly brittle failure of the interface. The average load–displacement behavior of all the test coupons with the standard
deviation is shown in Fig. 15.

6. Numerical studies

6.1. Validation of the CZM for ductile-brittle fracture transition at the interface

The numerical procedure was validated by considering a three-point bending problem, as shown in Fig. 16. A similar problem was
earlier solved by Carpinteri and Colombo [41] using the node release technique within FE analysis. The geometric dimensions of the
problem.
b = t = 150 mm, l = 600 mm, a = 45 mm

with the following material parameters

E = 36.5 GPa, σmax = 3.19 MPa, ν = 0.1

where E, σ max and ν were the Young’s modulus, cohesive strength, and the Poisson’s ratio of the material, respectively. A linear
softening cohesive law, as mentioned in Eq. (15) was incorporated in ABAQUS using the user subroutine UMAT.
The domain-independent interaction integral formulation was used to evaluate the SIF. The displacement along the x direction was
constrained at the point of load application. The parent material was meshed with 4-noded bilinear plane stress quadrilateral elements,
and 4-noded 2D cohesive element was used along the straight predefined starter crack, as shown in Fig. 17. The cohesive elements were
incrementally inserted as the crack propagates; the SIF at the cohesive tip was iteratively calculated by the procedure as presented by
Moes and Belytscko [42], however, for problems without a starter crack.
The displacement at the point of application of loading shown in Fig. 16 was monotonically increased to a level where the SIF at the
mathematical tip is zero within a reasonable limit. The interface crack propagation through the cohesive zone can be stated as.

Step 1: Initially, an arbitrary value of cohesive zone of length da = 15mm was considered ahead of the physical crack tip. A linear
softening cohesive law was assessed with an initial guess of displacement δ was applied at the point of application, as shown in
Fig. 16. The value of reaction force, P at the support point, was simultaneously calculated. While doing so, if the separation of any of
the cohesive elements was greater than the critical separation δc , eliminate the cohesive element, and the procedure is repeated.
Step 2: SIF at the cohesive tip was evaluated using the domain-indepenent interaction integral procedure.

Fig. 14. Fracture toughness testing experimental set up.

14
P.J. Saikia and N. Muthu Engineering Fracture Mechanics 266 (2022) 108353

Fig. 15. Load displacement behavior at the interface.

Fig. 16. Schematic of a three-point bending specimen.

Mathematical tip

Cohesive zone

Physical tip

Pre-crack

Fig. 17. Finite element mesh generation of the geometry (all dimensions are in mm).

Step 3: If the overall SIF was nullified within some reasonable limit, another incremental value of cohesive zone length da = 15mm
was added, step 1 and 2 was repeated. This procedure was carried out till the end of the cohesive zone length.

In most engineering structures, crack growth is usually steady and progresses slowly. In some other cases, crack formation and its
propagation is catastrophic, resulting in the abrupt loss of load-bearing capacity. The failure behavior of the structure might range
from ductile to brittle depending on their material properties, structural geometry and the loading condition. To investigate the effect

15
P.J. Saikia and N. Muthu Engineering Fracture Mechanics 266 (2022) 108353

of geometric dimensions (size effect) on the global load–displacement behavior and to compare the interface brittleness, a dimen­
sionless number (popularly known as brittleness number) was introduced by Carpinteri [43]. The brittleness number can be defined as.
GF
SE = (18)
σmax b

where GF is the fracture energy and σ max is the maximum cohesive traction. b is the depth of the specimen considered. The
load–deflection path presumed a positive slope for particularly brittle failure (specimen with high tensile strength and low toughness)
at the interfaces. This can only be possible with simultaneous degradation of load and displacement to achieve stable and controlled
crack propagation. On the other hand, a higher value of SE , which is due to a higher value of GF or a lower value of σmax, results in a
more ductile specimen. A ductile failure always proceeds with normal softening, resulting in a negative slope in the load–displacement
behavior. Three different values of brittleness number (SE ) were considered to capture the ductile to brittle failure transition behavior
at the interface and were shown in Figs. 18-20. The crack propagation procedure agreed well with the established results of Carpinteri
and Colombo [41] for all three cases. As the length of the crack propagation or, in this case, the cohesive zone length da decreases, a
smoother load–displacement curve will be obtained.

6.2. Numerical validation

The FE model of the specimen geometry is shown in Fig. 21. An initial crack of length a = r along the mid-plane was considered up
to the center point (x = 0, y = 0) of the specimen. The center of both the circular hole was considered a reference point for applying the
boundary conditions and evaluating the interface fracture toughness. The top reference point (RP-1) was constrained to move along
the x-direction, while the bottom (RP-2) was constrained in both x andy. The displacement control boundary condition was prescribed
along y-direction at the top reference point to induce crack propagation which was governed by the procedure mentioned in the earlier
section.
The ductile to brittle fracture transition behavior at the interface for three different values of brittleness number is shown in Fig. 22.
Maximum cohesive stress of σ max = 16MPa was considered for all the simulations. The geometric dimension b was taken care of by the
thickness of the specimen. It was observed that for a brittleness numberSE = 0.001467, the numerical results obtained from the present
procedure were very close to the results obtained from the experiments, as shown in Fig. 23. So, the current procedure can be effi­
ciently used to capture the ductile to brittle fracture transition behavior of the interface fracture mechanics.

7. Conclusions

A semi-analytical asymptotic stress and displacement field solution was obtained for pure mode I and mode II fracture analysis in
the presence of cohesive zone for any generalized cohesive law that nullifies the crack tip singularity, assuming no friction between the
crack surfaces. The solution coefficients were dependent on the type of cohesive law, material geometry, and the boundary conditions
of the problem. However, there was no need to evaluate the coefficients for a given cohesive zone length ρ of constant cohesive TSL.
The obtained SIFs using the domain-independent interaction integral ascertains the reduction in the strength of the singularity and
found to agree well with the analytical formulation for determining the SIF. The presented procedure could be easily extended to
higher-order extrinsic TSLs with finite traction at zero separation. The simulation of a cracked three-point bending specimen was of
considerable interest. The brittleness number related to cohesive law parameters and material geometry specifications affects the
load–displacement curve; thus, demonstrating the ability of the proposed procedure to model the ductile to brittle crack propagation
phenomenon. Finally, a pre-cracked disc-shaped CT specimen was fabricated and subjected to external loading in the UTM to obtain
the load–displacement curve due to interface crack propagation; this corresponded to brittle failure. In the FE analysis that considered

Fig. 18. Dimensionless load-displacement graph for SE = 0.0209.

16
P.J. Saikia and N. Muthu Engineering Fracture Mechanics 266 (2022) 108353

Fig. 19. Dimensionless load-displacement graph for SE = 0.000209.

Fig. 20. Dimensionless load-displacement graph for SE = 0.0000418.

Fig. 21. Finite element mesh of the CT geometry.

17
P.J. Saikia and N. Muthu Engineering Fracture Mechanics 266 (2022) 108353

Fig. 22. Comparison of results obtained for different brittleness number.

Fig. 23. Comparison of results obtained from the experiment and the present method.

a linear softening cohesive law, the load required for an incremental crack growth was determined by nullifying the SIF at the
mathematical crack tip to obtain the load–displacement curve. The FE results agreed with the experimental data demonstrating the
effectiveness of modeling interface crack propagation problem.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to
influence the work reported in this paper.

Acknowledgments

This work was supported by Science and Engineering Research Board, India under the scheme ‘Early Career Research Award’, No:
ECR/2018/001638.

Appendix I

The asymptotic displacement and the stress field near the cohesive crack tip can be defined as:.
∑rλn
u= {k(a1n cosλn θ − a2n sinλn θ) + λn [ − a1n cos(λn − 2)θ + a2n sin(λn − 2)θ ] + (λn + 1)( − b1n cosλn θ + b2n sinλn θ)} (A1.1)
n=0
2μ

∑rλn { [ ] }
v= k(a1n sinλn θ + a2n cosλn θ) + λn a1n sin(λn − 2)θ + a2n cos(λn − 2)θ + (λn + 1)(b1n sinλn θ + b2n cosλn θ) (A1.2)
n=0
2μ

18
P.J. Saikia and N. Muthu Engineering Fracture Mechanics 266 (2022) 108353

∑
σxx = rλn − 1 {2λn [a1n cos(λn − 1)θ − a2n sin(λn − 1)θ ] − λn (λn − 1)[a1n cos(λn − 3)θ − a2n sin(λn − 3)θ ] − λn (λn + 1)[b1n cos(λn − 1)θ
n=0

− b2n sin(λn − 1)θ]} (A1.3)

∑
σyy = rλn − 1 {2λn [a1n cos(λn − 1)θ − a2n sin(λn − 1)θ ] + λn (λn − 1)[a1n cos(λn − 3)θ − a2n sin(λn − 3)θ ] + λn (λn + 1)[b1n cos(λn − 1)θ
n=0

− b2n sin(λn − 1)θ]} (A1.4)

∑
τxy = rλn − 1 {λn (λn − 1)[a1n sin(λn − 3)θ + a2n cos(λn − 3)θ ] + λn (λn + 1)[b1n sin(λn − 1)θ + +b2n cos(λn − 1)θ ]} (A1.5)
n=0

Appendix II

The auxiliary displacement and stress field for interaction integral as follows.
√̅̅̅̅̅
r [ aux I ]
uaux
i = K f (θ) + KIIaux fiII (θ) (A2.1)
2π I i

1 [ ]
σaux aux I aux II
ij = √̅̅̅̅̅̅̅ KI gij (θ) + KII gij (θ) (A2.2)
2πr

k − cosθ θ 2 + k + cosθ θ
f1I (θ) = cos , f1II (θ) = sin (A2.3)
2μ 2 2μ 2

k − cosθ θ 2 − k − cosθ θ
f2I (θ) = sin , f2II (θ) = cos (A2.4)
2μ 2 2μ 2

3 θ 1 5θ 7 θ 1 5θ
gI11 (θ) = cos + cos , gII11 (θ) = − sin − sin , (A2.5)
4 2 4 2 4 2 4 2

5 θ 1 5θ 1 θ 1 5θ
gI22 (θ) = cos − cos , gII22 (θ) = − sin + sin , (A2.6)
4 2 4 2 4 2 4 2

gI12 (θ) = gI21 (θ) = gII22 (θ), gII12 (θ) = gII21 (θ) = gI11 (θ) (A2.7)

References

[1] Inglis CE. Stresses in a plate due to the presence of cracks and sharp corners. Trans Inst Nav Arch 1913;55:219–41.
[2] Griffith AAVI. The phenomena of rupture and flow in solids. Philos Trans R Soc London Ser A, Contain Pap a Math or Phys Character 1921;221:163–98.
[3] Barenblatt GI. The mathematical theory of equilibrium cracks in brittle fracture. Adv Appl Mech., vol. 7, Elsevier; 1962, p. 55–129.
[4] Dugdale DS. Yielding of steel sheets containing slits. J Mech Phys Solids 1960;8(2):100–4.
[5] Camacho GT, Ortiz M. Computational modelling of impact damage in brittle materials. Int J Solids Struct 1996;33(20-22):2899–938. https://doi.org/10.1016/
0020-7683(95)00255-3.
[6] Needleman A. A continuum model for void nucleation by inclusion debonding. J Appl Mech Trans ASME 1987;54:525–31. https://doi.org/10.1115/1.3173064.
[7] Needleman A. An analysis of decohesion along an imperfect interface. Int J Fract 1990;42(1):21–40. https://doi.org/10.1007/BF00018611.
[8] Xu X-P, Needleman A. Void nucleation by inclusion debonding in a crystal matrix. Model Simul Mater Sci Eng 1993;1(2):111–32. https://doi.org/10.1088/
0965-0393/1/2/001.
[9] Tvergaard V, Hutchinson JW. The relation between crack growth resistance and fracture process parameters in elastic-plastic solids. J Mech Phys Solids 1992;40
(6):1377–97.
[10] Shao Y, Zhao H-P, Feng X-Q, Gao H. Discontinuous crack-bridging model for fracture toughness analysis of nacre. J Mech Phys Solids 2012;60(8):1400–19.
https://doi.org/10.1016/j.jmps.2012.04.011.
[11] Jin Z-H, Sun CT. Cohesive fracture model based on necking. Int J Fract 2005;134(2):91–108. https://doi.org/10.1007/s10704-005-7864-1.
[12] Zhang Z(, Paulino GH, Celes W. Extrinsic cohesive modelling of dynamic fracture and microbranching instability in brittle materials. Int J Numer Methods Eng
2007;72(8):893–923.
[13] Jin Z-H, Paulino GH, Dodds RH. Cohesive fracture modeling of elastic–plastic crack growth in functionally graded materials. Eng Fract Mech 2003;70(14):
1885–912. https://doi.org/10.1016/S0013-7944(03)00130-9.
[14] Paulino GH, Zhang ZY. Dynamic fracture of functionally graded composites using an intrinsic cohesive zone model. Mater Sci Forum 2005;492–493:447–52.
https://doi.org/10.4028/www.scientific.net/msf.492-493.447.
[15] Budiansky B, Rice JR. Conservation laws and energy-release rates 1973.
[16] Stern M, Becker EB, Dunham RS. A contour integral computation of mixed-mode stress intensity factors. Int J Fract 1976;12(3):359–68. https://doi.org/
10.1007/BF00032831.
[17] Wang SS, Yau JF, Corten HT. A mixed-mode crack analysis of rectilinear anisotropic solids using conservation laws of elasticity. Int J Fract 1980;16(3):247–59.
https://doi.org/10.1007/BF00013381.

19
P.J. Saikia and N. Muthu Engineering Fracture Mechanics 266 (2022) 108353

[18] Yu H, Wu L, Guo L, Du S, He Q. Investigation of mixed-mode stress intensity factors for nonhomogeneous materials using an interaction integral method. Int J
Solids Struct 2009;46(20):3710–24. https://doi.org/10.1016/j.ijsolstr.2009.06.019.
[19] Walters MC, Paulino GH, Dodds RH. Interaction integral procedures for 3-D curved cracks including surface tractions. Eng Fract Mech 2005;72(11):1635–63.
https://doi.org/10.1016/j.engfracmech.2005.01.002.
[20] Daimon R, Okada H. Mixed-mode stress intensity factor evaluation by interaction integral method for quadratic tetrahedral finite element with correction terms.
Eng Fract Mech 2014;115:22–42. https://doi.org/10.1016/j.engfracmech.2013.11.009.
[21] Yu H, Kuna M. Interaction integral method for computation of crack parameters K-T – A review. Eng Fract Mech 2021;249:107722. https://doi.org/10.1016/j.
engfracmech.2021.107722.
[22] Ingraffea AR, Gerstk WH, Gergely P, Saouma V. Fracture mechanics of bond in reinforced concrete. J Struct Eng 1984;110(4):871–90. https://doi.org/10.1061/
(ASCE)0733-9445(1984)110:4(871).
[23] Cendón DA, Gálvez JC, Elices M, Planas J. Modelling the fracture of concrete under mixed loading. Int J Fract 2000;103:293–310. https://doi.org/10.1023/A:
1007687025575.
[24] Jin Z-H, Sun CT. Cohesive zone modeling of interface fracture in elastic bi-materials. Eng Fract Mech 2005;72(12):1805–17. https://doi.org/10.1016/j.
engfracmech.2004.09.011.
[25] Wang JT. Relating Cohesive Zone Models to Linear Elastic Fracture Mechanics. Natl Aeronaut Sp Adm Langley Res Cent 2010:20.
[26] Rice JR, Sih GC. Plane problems of cracks in dissimilar media. J Appl Mech 1965;32:418–23. https://doi.org/10.1115/1.3625816.
[27] Xiao QiZhi, Karihaloo BL. Asymptotic fields at frictionless and frictional cohesive crack tips in quasibrittle materials. J Mech Mater Struct 2006;1(5):881–910.
[28] Yu H, Wu L, Guo L, Wu H, Du S. An interaction integral method for 3D curved cracks in nonhomogeneous materials with complex interfaces. Int J Solids Struct
2010;47(16):2178–89. https://doi.org/10.1016/j.ijsolstr.2010.04.027.
[29] Yu H, Wu L, Guo L, He Q, Du S. Interaction integral method for the interfacial fracture problems of two nonhomogeneous materials. Mech Mater 2010;42(4):
435–50. https://doi.org/10.1016/j.mechmat.2010.01.001.
[30] Yu H, Wu L, Guo L, Li H, Du S. T-stress evaluations for nonhomogeneous materials using an interaction integral method. Int J Numer Methods Eng 2012;90(11):
1393–413. https://doi.org/10.1002/nme.4263.
[31] Kin Loch AJ, Shaw SJ. The fracture resistance of a toughened epoxy adhesive. J Adhes 1981;12(1):59–77.
[32] Banks-Sills L, Travitzky N, Ashkenazi D, Eliasi R. A methodology for measuring interface fracture properties of composite materials. Int J Fract 1999;99:143–60.
https://doi.org/10.1023/A:1018642200610.
[33] Banks-Sills L, Travitzky N, Ashkenazi D. Interface fracture properties of a bimaterial ceramic composite. Mech Mater 2000;32(12):711–22. https://doi.org/
10.1016/S0167-6636(00)00042-9.
[34] Banks-Sills L. Interface fracture mechanics: theory and experiment. Int J Fract 2015;191(1-2):131–46. https://doi.org/10.1007/s10704-015-9997-1.
[35] Muskhelishvili N. Some basic problems of the mathematical theory of elasticity. Noordhoff, Groningen. 1963;17404(6.2):1.
[36] Liebowitz H, Sih GC. Mathematical theories of brittle fracture 1968; 2 :67–190.
[37] Zavattieri PD, Hector LG, Bower AF. Cohesive zone simulations of crack growth along a rough interface between two elastic-plastic solids. Eng Fract Mech 2008;
75(15):4309–32. https://doi.org/10.1016/j.engfracmech.2007.11.007.
[38] User CAE. ABAQUS 6.14 Analysis User Guide. Dassault Syst 2014.
[39] Shih CF, Moran B, Nakamura T. Energy release rate along a three-dimensional crack front in a thermally stressed body. Int J Fract 1986;30(2):79–102. https://
doi.org/10.1007/BF00034019.
[40] Alam M, Parmigiani JP, Kruzic JJ. An experimental assessment of methods to predict crack deflection at an interface. Eng Fract Mech 2017;181:116–29. https://
doi.org/10.1016/j.engfracmech.2017.05.013.
[41] Carpinteri A, Colombo G. Numerical analysis of catastrophic softening behaviour (snap-back instability). Comput Struct 1989;31(4):607–36. https://doi.org/
10.1016/0045-7949(89)90337-4.
[42] Moës N, Belytschko T. Extended finite element method for cohesive crack growth. Eng Fract Mech 2002;69(7):813–33. https://doi.org/10.1016/S0013-7944
(01)00128-X.
[43] Carpinteri A. Notch sensitivity in fracture testing of aggregative materials. Eng Fract Mech 1982;16(4):467–81. https://doi.org/10.1016/0013-7944(82)90127-
8.

20