Composites Part C: Open Access 12 (2023) 100385

Contents lists available at ScienceDirect

Composites Part C: Open Access

journal homepage: www.sciencedirect.com/journal/composites-part-c-open-access

Simulation of mode II fracture propagation in adhesive joints using a

meshless technique
D.C. Gonçalves a, *, I.J. Sánchez-Arce a, b, L.D.C. Ramalho a, R.D.S.G. Campilho a, c, J. Belinha a, c
a
INEGI, Institute of Mechanical Engineering, Rua Dr. Roberto Frias, 4200-465, Porto, Portugal
b
UTAD, Universidade de Trás-os-Montes e Alto Douro, Quinta de Prados. Vila Real, 5000-801, Portugal
c
School of Engineering, Polytechnic of Porto, ISEP-IPP, Rua Dr. António Bernardino de Almeida, 431, 4200-072, Porto, Portugal

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

Keywords: Adhesive bonding is extensively used by commanding industries such as automotive and aircraft sectors.
Bonded joints Nevertheless, due to the intricate mechanical behaviour of adhesively bonded joints, especially when crack
Meshless method propagation occurs at the adhesive layer, improvement of new numerical techniques to simulate this bonding
Mode II fracture
approach is currently under development. In this work, a recent crack propagation algorithm based on a meshless
Crack propagation
Energy release rate
technique is implemented to analyse mode II fracture propagation in adhesively bonded joints. After the problem
domain is discretized with an independent set of field nodes, a background numerical integration mesh is
constructed. The crack tip advancement is then simulated by iteratively rearranging the field nodes and inte­
gration cells at the crack tip. Radial point interpolation (RPI) functions are constructed using the concept of
influence domains, allowing to obtain sparse and stable global matrices. To assess the suitability of the proposed
method, End-notched flexure (ENF) specimens were manufactured and tested. The numerical simulation justly
reproduces the experimental data provided and, thus, the numerical model can be applied to mode-II shear
loadings and industrial applications for design purposes.

1. Introduction constant number of nodes was employed rather than an iterative com­
plete remeshing procedure around the crack tip. Seminal techniques to
Adhesive bonding is extensively used by commanding industries model crack propagation by iterative remeshing finite element meshes
such as the automotive and aircraft sectors. Nevertheless, due to the can be found in references [6,7]. Nowadays, a wide range of more
intricate mechanical behaviour of adhesively bonded joints, especially complex remeshing procedures is available in recent literature. An
when crack propagation occurs at the adhesive layer, the improvement adaptative remeshing method using p-version FEM was developed by
of new numerical techniques to simulate this bonding approach is Wowk et al. [8] to predict the propagation of cracks with highly irreg­
currently under development. Although research on fracture propaga­ ular crack fronts. Recently, Ramalho et al. [9,10] proposed a local
tion can be traced back to the works of Inglis [1], Griffith [2] or Irwin remeshing procedure suited for finite element simulation of crack
[3], today, a wide range of analytical and numerical techniques are propagation in two-dimensional structures. Additionally, the same
available to predict fracture in engineering structures such as adhesively technique was combined with the meshless radial point interpolation
bonded joints. Analytical approaches are available to model fracture in method (RPIM) for fracture analysis [11,12].
adhesively bonded joints [4], although numerical techniques such as Even though recent FEM variations have been developed for crack
finite element methods (FEM) transcend the limited applicability of propagation, e.g., the extended Voronoi cell finite-element method (X-
analytical and closed-form derivations. VCFEM) [13], meshless methods currently represent a promising alter­
Crack propagation is commonly modelled resorting to advanced native to model fracture propagation problems. Meshless methods can
remeshing techniques of finite element models of the problem domain. be traced back to the primary smoothed particle hydrodynamics (SPH)
Nonetheless, alternative approaches have been developed, e.g., Crusat method [14], originally developed to model astrophysical phenomena.
et al. [5] applied configurational mechanics combined with the FEM Today, a plentiful of meshless techniques are available, either consid­
considering zero thickness interface elements, and a fixed mesh with ering strong or weak forms, or employing approximation or

* Corresponding author.
E-mail address: costa.goncalves.diogo@gmail.com (D.C. Gonçalves).

https://doi.org/10.1016/j.jcomc.2023.100385

Available online 2 August 2023

2666-6820/© 2023 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-
nc-nd/4.0/).
D.C. Gonçalves et al. Composites Part C: Open Access 12 (2023) 100385

Table 1 phase-field modelling has been developed as an efficient technique to

Adherend and adhesive materials properties [29–32] (a manufacturer’s data). model complex fracture problems, and meshless methods such as the
AW 6082-T651 Araldite® AV138 RPIM have been incorporated in such formulations [27], showing
interesting capabilities. Morever, meshless techniques such as the
Young’s modulus E [GPa] 70.07 ± 0.83 4.89 ± 0.81
Poisson’s ratio ν [− ] 0.33a 0.35a localized method of fundamental solutions (LMFS) [28], a meshless
Tensile yield strength σy [MPa] 261.67 ± 7.65 36.49 ± 2.47 collocation method, have also been developed to analyse crack propa­
Tensile failure strength σf [MPa] 324.00 ± 0.16 39.45 ± 3.18 gation with bimaterial interfaces.
Tensile failure strain εf [%] 21.70 ± 4.24 1.21 ± 0.10 As previously discussed, it can be concluded that meshless methods’
Shear modulus G [GPa] 26.9 1.56 ± 0.01
Shear yield strength τy [MPa] n.a 25.1 ± 0.33
interest to fracture propagation has been increasing recently. However,
Shear failure strength τf [MPa] n.a. 30.2 ± 0.40 their application to non-benchmark practical applications such as
Shear failure strain γf [%] n.a. 7.8 ± 0.7 adhesively bonded joints is still scarce, and deep studies are required to
develop suitable techniques to analyse fracture propagation in
demanding structures. In this work, it is intended to extend a fracture
interpolation shape functions. Popular approximating meshless methods
propagation algorithm, based on the RPIM and energy release rate cri­
include the diffuse element method (DEM) [15], the element-free
terion, to mode II fracture propagation in adhesively bonded joints, and
Galerkin method (EFGM) [16] or the reproducing kernel particle
assess its performance against experimental data.
method (RKPM) [17]. Although these techniques and their variations
are well established in the community, meshless methods implementing
2. Experimental work
interpolating shape functions present the delta Kronecker property
advantage, allowing the direct imposition of boundary conditions in the
2.1. Geometry and materials
global matrix system. The RPIM [18] used in this work is derived from
the addition of radial basis functions (RBF) to the original point inter­
The ENF adhesive joints assembled in this work were manufactured
polation method (PIM) [19] formulation. In this manner, singularity
using aluminium adherends, bonded with a brittle and strong epoxy
problems arising originally are circumvented, rendering the RPIM one of
adhesive material, the Araldite® AV138 (Araldite® from Huntsman,
the most accurate meshless interpolating techniques.
Basel, Switzerland). The aluminium used was the AA6082 T651 alloy,
In recent years, meshless methods have found their way into the
characterized by its high strength and ductility. Table 1 presents the
fracture propagation field due to their advantages over traditional
mechanical properties of the adherend and adhesive materials, previ­
techniques. Recently, an improved SPH method [20] was extended to
ously characterized by Campilho et al. [29]. The ENF specimen geom­
crack propagation in brittle fracture, showing numerical results in
etry is shown in Fig. 1, being the specimens’ targeted dimensions the
agreement with experimental data. In reference [21] the background
following: adherend thickness ts = 3mm, adhesive thickness ta = 0.2mm,
decomposition method (BDM) is considered to numerically integrate the
total length LT = 230mm, width w = 25mm, and a the initial crack
domain integrals of the RPIM and an enriched EFGM formulation,
length. The boundary conditions were imposed experimentally as shown
developed for crack tip problems. The benchmark analyses and the high
by Fig. 1 and Fig. 4. The specimens are placed on two rollers and an
accuracy results further demonstrated the potential of meshless methods
imposed displacement δ is applied at midspan.
in linear elastic fracture mechanics. Analysis of cracks in plates was also
investigated using the EFGM in the work of Peng et al. [22]. Zhuang
et al. [23] proposed an RPIM technique for linear fracture mechanics 2.2. Fabrication and test setup
using the Williams expansion as trial functions near the crack tip,
addressed as the sub-region radial point interpolation method The complete fabrication of the specimens was performed at room
(MS-RPIM). In comparison with other meshless techniques, it was found temperature. Firstly, a surface preparation procedure was carried out.
that the RPIM requires a less dense nodal discretization around the crack The aluminium adherends were grit blasted with corundum sand and
tip, while smooth and accurate stress fields are obtained. The coupled wiped with acetone. The curing procedure was achieved in a steel
extended meshfree–smoothed meshfree method (CXSMM) [24] was mould, guaranteeing the adherends’ alignment up to complete curing of
developed to predict crack growth in two-dimensional solids. It pos­ the adhesive. The desired ta within the bonded adhesive layer was
sesses both the assets of extended EFGM and smoothed meshless achieved by inserting calibrated steel spacers between the adherends.
methods allowing to accurately model linear elastic fracture. An RPIM Also, a demoulding agent was applied to posteriorly ease its removal.
approach was also used to analyse linear fracture in viscoelastic mate­ The pre-crack (a0) was induced at the adhesive layer using a 0.1mm
rials [25], demonstrating the efficiency of the method to model complex razor blade between two 0.05mm spacers. After applying the adhesive,
wave propagation problems with multiple crack tips in linear visco­ the specimens were closed carefully to circumvent air entrapment and
elasticity environments. Moreover, Hamidpour et al. [26] likewise avoid air bubbles within the adhesive layer. Finally, curing was achieved
considered the RPIM to model fracture in viscoelasticity, and numerical at room temperature. Fabricated specimens during collage are shown in
predictions close to experimental tests were also obtained. Recently, Fig. 2.
Before testing, the steel spacers were removed and the specimens

Fig. 1. ENF test geometry and boundary conditions.

2
D.C. Gonçalves et al. Composites Part C: Open Access 12 (2023) 100385

Fig. 2. Fabricated specimens during cure.

Fig. 3. Assembled specimen with measuring scale.

Fig. 4. Mechanical test setup.

3
D.C. Gonçalves et al. Composites Part C: Open Access 12 (2023) 100385

Fig. 5. Basic influence domain procedure (a) and influence domain of a node near the crack (b).

were sprayed with a brittle white paint allowing clear identification of

3.1.3. CBBM
the crack tip path. Additionally, a measuring scale was glued to the
In opposition to the previously described methods to measure GII, a
adherends, as shown in Fig. 3, to assist the crack length measurement.
new technique independent of a measurement during propagation is
The mechanical tests were performed using a Shimadzu AG-X 100
hereby presented, the Compliance Based Calibration Method (CBBM)
(Shimadzu, Kyoto, Japan) electro-mechanical testing machine, equip­
[33]. Moreover, the CBBM techniques takes into account a Fracture
ped with a 100kN load cell and using a 0.5mm/s test velocity. In order to
Process Zone (FPZ) arising at the crack tip. GII is evaluated using,
capture the evolution of a during fracture propagation, images of the
specimens at approximately 100mm were captured every five seconds 9P2 a2eq
(Fig. 4), allowing to posteriorly synchronize the a evolution with the GII = (6)
16Ew2 ts3
imposed value of δ.
being aeq an equivalent crack length given by,
3. Methods [ ( ) ]13
Cc 3 2 Cc
aeq = a0 + − 1 LT3 (7)
3.1. Mode II data reduction methods C0 C 3 C0 C

3.1.1. DBT where LT is the specimen’s total length, Cc the corrected compliance, and
The Direct Beam Theory (DBT) method is the most straightforward C0C the corrected compliance for a = a0, obtained as,
technique to obtain the mode II fracture energy (GII) and is obtained as, 3L
Cc = C − (8)
9P2 a2 10GBh
GII = (1)
16Ew2 ts3
with
being the P the applied load and E the elasticity modulus. 3a3 + 2L3 3L
C= + . (9)
8Ewts3 10Gwts
3.1.2. CBT
Although the DBT method provides a direct estimation of the energy Using an apparent elasticity modulus, Ef, allows to take into account
release rate, more complex techniques, such as the Corrected Beam the changeability of the material properties between different specimens
Theory (CBT), have been developed. The general expression is derived and several effects such as stress concentration near the crack tip,
by the DBT formula, but a corrected crack length affected by the factor 3a30 + 2L3T
ΔII is considered. The energy release rate is then given by Ef = . (10)
8wts3 C0c

GII =
9P2 (a + ΔII )2
(2) Full details on the CBBM development and formulation can be con­
16Ew2 ts3 sulted in the seminal work of Morais et al. [33].
The correction factor ΔII is related to ΔI for mode I,
ΔII = 0.42⋅ΔI (3) 3.2. Meshless method

which is given by, The meshless RPIM is used to numerically analyse the cracked
√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
[ ̅ structure during fracture propagation. In opposition to the FEM, where
( Γ )2 ]
E the numerical model is discretized by rigid elements naturally con­
ΔI = 0.42⋅h 3− 2 (4)
11G 1+Γ necting the field nodes, in the RPIM, the numerical model is discretized
with a set of field nodes with no predetermined connectivity. None­
being theless, a background integration mesh is necessary to obtain the inte­
E gration points required to numerically integrate the RPI shape functions.
Γ = 1.18 (5) In this work, triangular integration cells are used considering the inte­
G
gration scheme proposed in reference [34]. First, each triangular cell is

4
D.C. Gonçalves et al. Composites Part C: Open Access 12 (2023) 100385

divided into three independent quadrilateral subcells, allowing Gauss where Φ(xI) is the shape function vector and Ψ(xI) is a residual vector
integration scheme to be performed at each subcell. Then, influence with no physical meaning that can be completely negligible. The partial
domains may be determined for every integration point, as shown by derivatives of Φ(xI) are easily computed [37]. The RPI shape functions
Fig. 5a, allowing nodal connectivity to be imposed by the overlap of possess the Kronecker delta property, which allows simple direct
influence domains of distinct integration points, e.g., xI,xJ. Influence imposition of the boundary conditions [37].
domains represent the set of nodes contributing to the shape functions
and local stiffness matrix construction of a given integration point. As 3.2.2. Discrete system of equations
Fig. 5a illustrates, influence nodes are obtained by radial searching the n Considering a solid with domain Ω ∈ R and boundary Γ, the linear-
closest neighbouring nodes. Studies on the optimal value of n can be static equilibrium equations are defined by:
found in early RPIM literature. However, in this work, no predetermined
∇σ + b = 0, (16)
value must be considered most advantageous as it may vary significantly
with the practical application and accuracy/computational cost where σ is the Cauchy stress tensor and b the body force per unit volume
compromise. In this work, n = 9 nodes inside the influence domains are vector. The boundary conditions are given by:
considered, since it is within the lower but acceptable range recom­
mended by the literature, and allows satisfactory computational cost. σn = t onnaturalboundary Γt
(17)
In the presence of a crack, the concept of influence domains to u = u onessentialboundaryΓu
impose nodal connectivity near the crack tip presents a drawback, since
such influence domains must be restricted to account for the crack tip where t is the traction force on Γt and u the prescribed displacement on
discontinuity. In the presence of a crack, the influence domain near one Γu. Considering Eq. (16), the Galerkin weak form can be derived and
edge of the crack should not contain influence nodes belonging to the presented as:
opposite side of the crack, as the numerical model would ignore the ∫ ∫ ∫
presence of the discontinuity. Hence, the influence domain of xI is L = δεT σdΩ − δuT bdΩ − δuT tdΓt = 0. (18)
restricted to field nodes of the same crack side domain, as shown by Ω Ω Γt

Fig. 5b. Assuming the two-dimensional formulation, the displacement vector

of a given field node is defined by two degrees of freedom: u(xI) = {u(xI),
3.2.1. RPI shape functions v(xI)}T. Using the shape function approximation defined in Eq. (15) for a
With the integration points and influence domains established, the given integration point, the displacement vector can be defined by:
RPI shape functions may be calculated. The interpolation u(xI), defined ⎧ ⎫
within the n nodes inside the influence domain of xI, is constructed by ⎪
⎪
⎪ ⎪
u1 ⎪
⎪
{ } [ ]⎪
⎨ v1 ⎪
adding a RBF to a polynomial basis function (PBF): u(xI ) φ1 0 ⋯ φn 0
⎬
{ } u(xI ) = = ⋮ (19)
∑n ∑m
{ } a(xI ) v(xI ) 0 φ1 ⋯ 0 φn ⎪ ⎪ ⎪ ⎪
⎪
⎪ un ⎪
⎪
u(xI ) = Ri (xI )⋅ai (xI ) + pj (xI )⋅bj (xI ) = RT (xI ), pT (xI ) H(xI ) ⎩ ⎭
b(xI ) vn
i=1 j=1
us
(11)
The deformation vector ε(xI) can be obtained through derivation of
Different PBF can be used depending on the number of monomials m. the displacement field:
In this work, a constant PBF is considered: p(xI) = {1}. Although several ⎡ ⎤
RBF are available in the literature, this work used the Multi-Quadrics ∂
0
RBF (MQ-RBF), originally proposed by Hardy [35]. The MQ-RBF is ⎧
⎢
⎫ ⎢ ∂x
⎥
⎥
defined as: ⎨ εxx (xI ) ⎬ ⎢
⎢
⎥{
∂ ⎥ u(xI )
}
⎢
ε(xI ) = εyy (xI ) = ⎢ 0 ⎥ (20)
( )p ⎩ γ (x ) ⎭ ⎢ ∂y ⎥
⎥ v(xI )
R(dIi ) = dIi2 + c2 , (12) xy I ⎢ ⎥
⎣∂ ∂⎦
being dIi the distance dIi between xI and the nodes in the influence ∂y ∂x
domain, c and p constant shape parameters that may be optimized. Early L

works show that c = 0.0001 and p = 0.9999 provide accurate results in

which, after combining Eqs. (19) and (20), can be written as:
solid mechanics applications [34]. The coefficients ai(xI) and bj(xI) of ⎡ ⎤
the RBF and PBF are obtained by constructing a system of equations that ∂φ1 ∂φn ⎧ ⎫
0 … 0
forces the interpolation function u(xI) in Eq. (11) to pass through the n ⎧
⎢
⎫ ⎢ ∂x ∂x ⎥⎪ u1 ⎪
⎥⎪
⎪ ⎪
⎢ ⎥⎪ ⎪ ⎪
field nodes in the influence domain of xI: ⎨ εxx (xI ) ⎬ ⎢ ∂φ1 ∂φn ⎥⎨ v1 ⎬
ε(xI ) = εyy (xI ) = ⎢ 0 … 0 ⎥ ⋮ (21)
{ } ⎩ γ (x ) ⎭ ⎢ ⎢ ∂y ∂y ⎥
⎥⎪
⎪ un ⎪
⎪
a(xI ) ⎥⎪
⎪ ⎪
us = [ R P ] . (13) xy I ⎢
⎣ ∂φ1 ∂φ1 ∂φn ∂φn ⎦⎩ vn ⎭
⎪
b(xI ) …
∂y ∂x ∂y ∂x us
To obtain a system with n + m equations and n + m unknowns, an B(xI )
extra set of m equations is added to the system [36]: { }T
{ } [ ]{ } The stress field σ = σ xx σ yy τxy can be defined by the material
us R P a(xI )
= T (14) constitutive relation σ = cε = cBus, being c the material constitutive
0 P 0 b(xI )
M
matrix. The Galerkin weak form can then be written as:
∫ ∫ ∫ ΩI
Coefficients a(xI) and b(xI) can now be calculated through inversion δB(xI )T cB(xI )us dΩ − δH(xI )T bdΩ − δH(xI )T tdΓt = 0 (22)
of the moment matrix M and substituted in Eq. (11), which can now be Ω Ω Γt
written as: Ki f bi f ti

{ } { }
{ }
u(xI ) = RT (xI ), pT (xI ) M− 1
us u
= {Φ(xI ), Ψ(xI )} s = Φ(xI )us , and assembled in a global system of algebraic equations in the form:
0 0
(15) δuT [Ku − Fb − Ft ] = 0 ⇔ Ku = f b + f t . (23)
The boundary conditions can be enforced directly in global stiffness

5
D.C. Gonçalves et al. Composites Part C: Open Access 12 (2023) 100385

Fig. 6. Algorithm flowchart.

matrix since the RPI shape functions possess the Kronecker delta material properties, RPIM parameters, and fracture parameters. Then,
property. the RPIM analysis takes place. Considering the current crack propaga­
tion remeshing increment, the integration points, influence domains and
shape functions are updated. The discrete system of equations is then
3.3. Crack growth algorithm assembled, and the structure response is obtained by solving Ku = f. The
crack initiation load at the current step is then estimated using the
In this section, the fracture propagation algorithm is briefly critical energy release rate criterion. The energy release rate at the
described. The local remeshing algorithm allowing crack propagation current step is obtained using the data reduction schemes presented in
simulation is based on the works of Ramalho et al. [11,12], but applied section 7 and compared with its critical value, (GIIC)
to the RPIM. Using this technique, the integration cells and integration
GII = GIIC (24)
points are incrementally rearranged to account for the crack tip propa­
gation. Fig. 6 presents a flow diagram of the implemented fracture
being GIIC derived experimentally for the respective adhesive. When the
propagation algorithm. Initially, the numerical model is built and pa­
propagation criterion is satisfied, the crack tip is advanced using the
rameters defined: nodal discretization, visibilities for influence domains,

Fig. 7. Experimental P − δ curves.

6
D.C. Gonçalves et al. Composites Part C: Open Access 12 (2023) 100385

Table 2 results are discussed and compared with the experimental data,
Specimen’s initial crack length. permitting to assess the suitability of the proposed numerical method­
Initial crack length, a [mm] ology to predict the mode II fracture propagation in adhesive joints.
Specimen 1 65.90
Specimen 2 51.65 4.1. Experimental data
Specimen 3 50.60
Specimen 4 −
Specimen 5 66.10 Fig. 7 shows the P δ curves of all the specimen’s tested. The curves
Specimen 6 64.15 with higher failure loads (specimens 2 and 3) achieve higher peaks than
Specimen 7 67.80 remaining samples due to the lower initial crack length (Table 2), which
Specimen 8 65.15
increases overall structure stiffness. The test of specimen 4 was consid­
ered invalid as it deviates from the remaining data or expected ENF
local remeshing algorithm, fully detailed in references [11,12], and a adhesive joint behaviour.
new RPIM analysis is performed. Fig. 8 presents the resistance curves obtained using the DBT (Fig. 8a),
CBT (Fig. 8b) and CBBM (Fig. 8c) theories presented in Section 3.1. Data
4. Results from specimens 6, 7 and 8 was found to be the most reliable and
coherent, and are therefore hereby presented. Analysis of the resistance
The obtained results are presented in this section. First, an experi­ curves of Fig. 8 allows to obtain GIIC, defined by the stable value of GII
mental evaluation of the ENF tests is carried out, allowing to fully after crack initiation. The retrieved values of GIIC are shown in Fig. 8 by
characterise the ENF adhesive joints and obtain the experimental frac­ the respective horizontal lines and are also presented in Table 3 for each
ture properties. Posteriorly, the correspondent numerical simulation specimen and data reduction scheme.
In order to proceed to the numerical simulation of the ENF adhesive

Fig. 8. Experimental resistance curves for DBT (a), CBT (b) and CBBM (c) theories.

7
D.C. Gonçalves et al. Composites Part C: Open Access 12 (2023) 100385

Table 3 adhesive joint model and its nodal discretization are presented in Fig. 9.
Critical energy release rate GIIC[N/mm]. The specimen dimensions are those of Fig. 1. The value of a considered
Theory DBT CBT CBBM in the numerical model was obtained by averaging the a0 measured
Specimen experimentally. Therefore, a0 = 66 mm was considered, which corre­
6 0.4606 0.4696 0.5010 sponds to the average of the initial crack length of specimens 1, 5, 6, 7
7 0.4603 0.4695 0.5101 and 8, presenting the basis for comparison purposes. In this manner, a
8 0.4780 0.4887 0.5269 unique numerical model can be simulated and the results consistently
Mean 0.4663 0.4756 0.5127 compared to the experimental data. The nodal discretization model
Standard Deviation (S.D.) 0.0083 0.0086 0.0108
presented in Fig. 9 is composed by 10,513 field nodes and 20,612
integration cells for numerical integration of the RPIM shape functions.
tests, a critical GIIC must be defined, which was obtained by averaging Notice that the nodal discretization is significantly denser at the crack
the values of each method (Table 3). The numerical work was carried tip region and gets sparser at the specimen’s end, allowing reliable
considering the DBT, CBT and CBBM methods distinctively and the discretization at the concern regions and reducing computational
following values were considered within the numerical algorithm, GDBT expense. The deformed configuration of the ENF model is also shown in
IIC
= 0.4663N/mm, GCBT CBBM Fig. 9, where the deformed in shear arising at the crack tip and speci­
IIC = 0.4756N/mm and GIIC = 0.5127N /mm.
mens’ end is visible.
The numerical approximation of the load-displacement curve is
4.2. Numerical work presented in Fig. 10. The numerical model shows a slightly lower stiff­
ness, nonetheless the failure loads are accurately predicted using the
The numerical results are presented in the following section. The

Fig. 9. Model nodal discretization.

Fig. 10. Load-displacement predictions.

8
D.C. Gonçalves et al. Composites Part C: Open Access 12 (2023) 100385

Table 4 with slightly higher failure load prediction. It can also be concluded by
Failure loads Pmax[N]. Fig. 10 that the numerical model predicts the ENF joint mechanical
Experimental DBT CBT CBBM behaviour during crack propagation. The compliance of the adhesive
joint during the ENF test is presented in Fig. 11 related to the crosshead
Pmax[N] 445.3 474.4 473.1 457.0
Error[%] – 6.52 6.23 2.62 displacement imposed. The numerical approximation is again accurate
and precisely reproduces the experimental data.
Using the numerical load and displacement data, the correspondent
presented meshless technique. Table 4 presents the numerical failure numerical resistance curves were obtained using the data reduction
loads and the respective error regarding the experimental value, ob­ schemes considered. Fig. 12 shows the numerical prediction of the
tained by averaging the failure loads of the tested specimens. The DBT resistance curves for the DBT, CBT and CBBM theories. It can be noticed
and CBT present similar errors of 6.52% and 6.23%, respectively, being that the GIIC values from the numerical curves correspond to the values
the most accurate method the CBBM, with a 2.62% deviation. Moreover, introduced in numerical algorithm, and deriving out of Table 3, enabling
after crack initiation the CBBM presents a closer approximation to the the positive validation of the numerical approach.
experimental data. The DBT and CBT schemes show overlapped results In spite of the advantages of the meshless technique implemented in

Fig. 11. Specimen’s compliance.

Fig. 12. Numerical resistance curves.

9
D.C. Gonçalves et al. Composites Part C: Open Access 12 (2023) 100385

Fig. 13. Computational cost of the fracture analysis: RPIM vs FEM.

this work, the computational cost of the analysis represents the major experimental data further validated the numerical model. Hereby, the
drawback of the current method, as Fig. 13 demonstrates the compari­ RPIM is demonstrated to be an accurate computational tool to model
son of a FEM and meshless RPIM analysis for the equivalent fracture mode II fracture propagation in adhesively bonded joints. Despite the
propagation algorithm. The analyses were performed in an Intel i9 7900 higher computational cost compared with the FEM for an equal nodal
× 4 GHz with 64Gb RAM. Despite the higher computational cost of the density, the inherent meshless assets present interesting capabilities of
RPIM analysis in comparison with straightforward finite element tech­ the RPIM, allowing combination and development of accurate and
nique, such cost can be reduced by lowering the nodal density of the efficient fracture propagation algorithms.
model, without severely compromising the numerical accuracy of the
meshless method, which is known to present high rates of convergence, Declaration of Competing Interest
as opposed to finite element techniques.
The authors declare that they have no known competing financial
5. Conclusions interests or personal relationships that could have appeared to influence
the work reported in this paper.
This work presented the application of a crack propagation algo­
rithm based on the RPIM to the mode II fracture propagation in adhe­ Data availability
sively bonded joints. The crack tip propagation was computationally
achieved by iteratively remeshing the field nodes and integration cells/ The data that has been used is confidential.
points at the crack tip region. The RPI shape functions were calculated at
the newer increments, and then the mechanical response of the structure
was obtained. The crack initiation was predicted by a GIIC criterion, Acknowledgements
calculated using distinct fracture energy theories.
Foremost, an experimental assessment of ENF adhesive joint tests This work has been funded by the Ministério da Ciência, Tecnologia e
was executed. Experimental results present characteristic ENF test Ensino Superior through the Fundação para a Ciência e a Tecnologia
behaviour, demonstrated good reproducibility and provided a well- (Portugal), under project funding ‘2022.13841.BD’, ‘MIT-EXPL/ISF/
founded basis to compare the obtained algorithm results and validate 0084/2017′, ‘POCI-01–0145-FEDER-028351′, and ‘SFRH/BD/147628/
the numerical model. Considering the DBT, CBT, and CBBM data 2019′. Additionally, the authors acknowledge the funding provided by
reduction schemes, resistance curves were obtained, and it was the Associated Laboratory for Energy, Transports and Aeronautics
demonstrated that the brittle adhesive AV138 presents GIIC values of (LAETA), under project ‘UIDB/50022/2020′.
0.4663N/mm, 0.4756N/mm and 0.5127N/mm, respectively.
Regarding the numerical application, the local remeshing algorithm References
permitted obtaining refined and well-conditioned nodal discretizations
around the crack tip, stably propagate the crack tip, and therefore in­ [1] C.E. Inglis, Stresses in a plate due to the presence of cracks and sharp corners, Proc.
Inst. Nav. Archit. 60 (1913) 219–241.
crease the RPIM analysis accuracy at each crack increment before
[2] A.A. Griffits, The phenomena of rupture and flow in solids, Phil. Trans. Roy. Soc.
specimen’s failure. The energy criterion applied was proved to be suit­ London A 221 (1920) 163–197.
able as the maximum failure loads are accurately predicted by the nu­ [3] G.R. Irwin, Analysis of stresses and strains near the end of a crack traversing a
merical simulation with errors of 6.52%, 6.23% and 2.62%, regarding plate, J. Appl. Mech. 24 (1957) 361–364.
[4] M. Brandtner-Hafner, Structural safety evaluation of adhesive bonds: a fracture
the DBT, DBT and CBBM methods, respectively. Moreover, the experi­ analytical approach, Eng. Fail. Anal. 123 (2021), 105289, https://doi.org/
mental P − δ behaviour is also closely approximated by the presented 10.1016/j.engfailanal.2021.105289.
methods, the CBBM being the technique that provides the most accurate [5] L. Crusat, I. Carol, D. Garolera, Application of configurational mechanics to crack
propagation in quasi-brittle materials, Eng. Fract. Mech. 241 (2020), 107349,
results. The obtained resistance curves computed against the https://doi.org/10.1016/j.engfracmech.2020.107349.

10
D.C. Gonçalves et al. Composites Part C: Open Access 12 (2023) 100385

[6] P.O. Bouchard, F. Bay, Y. Chastel, I. Tovena, Crack propagation modelling using an [22] L.X. Peng, et al., Simulation of a crack in stiffened plates via a meshless formulation
advanced remeshing technique, Comput. Methods Appl. Mech. Eng. 189 (3) (2000) and FSDT, Int. J. Mech. Sci. 131–132 (2017) 880–893, https://doi.org/10.1016/j.
723–742, https://doi.org/10.1016/S0045-7825(99)00324-2. ijmecsci.2017.07.063.
[7] P.O. Bouchard, F. Bay, Y. Chastel, Numerical modelling of crack propagation: [23] X. Zhuang, Y. Cai, C. Augarde, A meshless sub-region radial point interpolation
automatic remeshing and comparison of different criteria, Comput. Methods Appl. method for accurate calculation of crack tip fields, Theor. Appl. Fract. Mech. 69
Mech. Eng. 192 (35–36) (2003) 3887–3908, https://doi.org/10.1016/S0045-7825 (2014) 118–125, https://doi.org/10.1016/j.tafmec.2013.12.003.
(03)00391-8. [24] W. Ma, G. Liu, W. Wang, A coupled extended meshfree–smoothed meshfree
[8] D. Wowk, L. Alousis, P.R. Underhill, An adaptive remeshing technique for method for crack growth simulation, Theor. Appl. Fract. Mech. 107 (2020),
predicting the growth of irregular crack fronts using p-version finite element 102572, https://doi.org/10.1016/j.tafmec.2020.102572.
analysis, Eng. Fract. Mech. 207 (2019) 36–47, https://doi.org/10.1016/j. [25] M. Hamidpour, M.R. Nami, A. Khosravifard, An effective crack identification
engfracmech.2018.12.002. method in viscoelastic media using an inverse meshfree method, Int. J. Mech. Sci.
[9] L.D.C. Ramalho, J. Belinha, R.D.S.G. Campilho, A New Crack Propagation 212 (2021), 106834, https://doi.org/10.1016/j.ijmecsci.2021.106834.
Algorithm Combined with the Finite Element Method, J. Mech. 36 (4) (2020) [26] M. Hamidpour, M.R. Nami, A. Khosravifard, M. Lévesque, Modeling fracture in
405–422, https://doi.org/10.1017/jmech.2020.1. viscoelastic materials using a modified incremental meshfree RPIM and DIC
[10] L.D.C. Ramalho, J. Belinha, R.D.S.G. Campilho, A novel robust remeshing finite technique, Eur. J. Mech. A/Solids 92 (2022), https://doi.org/10.1016/j.
element technique for fracture propagation, Int. J. Comput. Methods 18 (2) (2021), euromechsol.2021.104456.
2050040, https://doi.org/10.1142/S0219876220500401. [27] L. Novelli, L. Gori, R.L. da Silva Pitangueira, Phase-field modelling of brittle
[11] L.D.C. Ramalho, J. Belinha, R.D.S.G. Campilho, Fracture propagation using the fracture with Smoothed Radial Point Interpolation Methods, Eng. Anal. Bound.
radial point interpolation method, Fatigue Fract. Eng. Mater. Struct. (2019) 1–15, Elem. 138 (2022) 219–234, https://doi.org/10.1016/j.enganabound.2022.01.011.
https://doi.org/10.1111/ffe.13046. [28] X. Wang, Y. Gu, M.V. Golub, Analysis of bimaterial interface cracks using the
[12] L.D.C. Ramalho, J. Belinha, R.D.S.G. Campilho, The numerical simulation of crack localized method of fundamental solutions, Results Appl. Math. 13 (2022),
propagation using radial point interpolation meshless methods, Eng. Anal. Bound. 100231, https://doi.org/10.1016/j.rinam.2021.100231.
Elem. 109 (2019) 187–198, https://doi.org/10.1016/j.enganabound.2019.10.001. [29] R.D.S.G. Campilho, M.D. Banea, A.M.G. Pinto, L.F.M. Da Silva, A.M.P. De Jesus,
[13] H. Li, R. Guo, H. Cheng, Extended Voronoi cell finite element method for multiple Strength prediction of single- and double-lap joints by standard and extended finite
crack propagation in brittle materials, Theor. Appl. Fract. Mech. 109 (2020), element modelling, Int. J. Adhes. Adhes. 31 (5) (2011) 363–372, https://doi.org/
102741, https://doi.org/10.1016/j.tafmec.2020.102741. 10.1016/j.ijadhadh.2010.09.008.
[14] R.A. Gingold, J.J. Monaghan, Smoothed particle hydrodynamics: theory and [30] J.A.B.P. Neto, R.D.S.G. Campilho, L.F.M. Da Silva, Parametric study of adhesive
application to non-spherical stars, Mon. Not. R. Astron. Soc. 181 (1977) 375–389. joints with composites, Int. J. Adhes. Adhes. 37 (2012) 96–101, https://doi.org/
[15] B. Nayroles, G. Touzot, P. Villon, Generalizing the finite element method: diffuse 10.1016/j.ijadhadh.2012.01.019.
approximation and diffuse elements, Comput. Mech. 10 (5) (1992) 307–318, [31] R.D.S.G. Campilho, D.C. Moura, D.J.S. Gonçalves, J.F.M.G. Da Silva, M.D. Banea, L.
https://doi.org/10.1007/BF00364252. F.M. Da Silva, Fracture toughness determination of adhesive and co-cured joints in
[16] T. Belytschko, Y.Y. Lu, L. Gu, Element-free Galerkin methods, Int. J. Numer. natural fibre composites, Compos. Part B 50 (2013) 120–126, https://doi.org/
Methods Eng. 37 (2) (1994) 229–256, https://doi.org/10.1002/nme.1620370205. 10.1016/j.compositesb.2013.01.025.
[17] L.Wing Kam, S. Jun, Z.Yi Fei, Reproducing kernel particle methods, Int. J. Numer. [32] T. Faneco, R. Campilho, F. Silva, R. Lopes, Strength and fracture characterization of
Methods Fluids 20 (8–9) (1995) 1081–1106. a novel polyurethane adhesive for the automotive industry, J. Test. Eval. 45 (2)
[18] J.G. Wang, G.R. Liu, A point interpolation meshless method based on radial basis (Mar. 2017) 398–407, https://doi.org/10.1520/JTE20150335.
functions, Int. J. Numer. Methods Eng. 54 (11) (2002) 1623–1648, https://doi.org/ [33] M.F.S.F. de Moura, A.B. de Morais, Equivalent crack based analyses of ENF and ELS
10.1002/nme.489. tests, Eng. Fract. Mech. 75 (9) (2008) 2584–2596, https://doi.org/10.1016/j.
[19] G.R. Liu, Y.T. Gu, A point interpolation method for two-dimensional solids, Int. J. engfracmech.2007.03.005.
Numer. Methods Eng. 50 (4) (2001) 937–951, https://doi.org/10.1002/1097-0207 [34] J. Belinha, Meshless Methods in Biomechanics - Bone Tissue Remodelling Analysis,
(20010210)50:4<937::AID− NME62>3.0.CO;2-X. 1st ed., Springer International Publishing Switzerland, Porto, 2014, p. 2014.
[20] D. Mu, A. Wen, D. Zhu, A. Tang, Z. Nie, Z. Wang, An improved smoothed particle [35] R.L. Hardy, Theory and applications of the multiquadric-biharmonic method,
hydrodynamics method for simulating crack propagation and coalescence in brittle Comput. Math. with Appl. 19 (8–9) (1990) 163–208, https://doi.org/10.1016/
fracture of rock materials, Theor. Appl. Fract. Mech. 119 (April) (2022), 103355, 0898-1221(90)90272-L.
https://doi.org/10.1016/j.tafmec.2022.103355. [36] M.A. Golberg, C.S. Chen, H. Bowman, Some recent results and proposals for the use
[21] A. Khosravifard, M.R. Hematiyan, T.Q. Bui, T.V. Do, Accurate and efficient analysis of radial basis functions in the BEM, Eng. Anal. Bound. Elem. 23 (4) (1999)
of stationary and propagating crack problems by meshless methods, Theor. Appl. 285–296, https://doi.org/10.1016/s0955-7997(98)00087-3.
Fract. Mech. 87 (2017) 21–34, https://doi.org/10.1016/j.tafmec.2016.10.004. [37] J. Belinha, Meshless Methods in Biomechanics - Bone Tissue Remodelling Analysis,
1st ed., Springer, Cham, Heidelberg, New York, Dordrecht, London, 2014.

11