Computational Mechanics manuscript No.
(will be inserted by the editor)
Finite Element Modelling of Internal and Multiple Localized
Cracks
Savvas Saloustros · Luca Pelà · Miguel Cervera · Pere Roca
Received: date / Accepted: date
Abstract Tracking algorithms constitute an efficient
numerical technique for modelling fracture in quasi-
brittle materials. They succeed in representing localized
cracks in the numerical model without mesh-induced di-
rectional bias. Currently available tracking algorithms
have an important limitation: cracking originates ei-
ther from the boundary of the discretized domain or
from predefined “crack-root” elements and then prop-
agates along one orientation. This paper aims to cir-
cumvent this drawback by proposing a novel tracking
algorithm that can simulate cracking starting at any
point of the mesh and propagating along one or two
orientations. This enhancement allows the simulation of
structural case-studies experiencing multiple cracking.
The proposed approach is validated through th e simu-
lation of a benchmark example and an experimentally
tested structural frame under in-plane loading. Mesh-
bias independency of the numerical solution, computa-
tional cost and predicted collapse mechanisms with and
without the tracking algorithm are discussed.
Savvas Saloustros · Luca Pelà · Miguel Cervera · Pere Roca
Department of Civil and Environmental Engineering
Technical University of Catalonia, UPC-BarcelonaTech
Jordi Girona 1-3, 08034 Barcelona, Spain.
International Center for Numerical Methods in Engineering
(CIMNE)
Gran Capità, S/N, 08034 Barcelona, Spain.
E-mail: savvas.saloustros@upc.edu
E-mail: luca.pela@upc.edu
E-mail: miguel.cervera@upc.edu
E-mail: pere.roca.fabregat@upc.edu
Keywords Continuum Damage Mechanics · Crack-
Tracking · Damage Localization · Quasi-brittle
materials · Shear/Flexural/Tensile Cracks
1 Introduction
Almost half a century after the pioneering works of Ngo
& Scordelis [
1] and Rashid [2], the numerical simulation
of cracking in quasi-brittle materials is still a challeng-
ing task. Although a wide range of novel formulations
has been proposed in the field of computational fail-
ure mechanics, up to date there is no such thing as a
panacea for the realistic and efficient numerical analy-
sis of failure in quasi-brittle materials. As a fact, the
analyst has to consider at least three aspects before
choosing the adequate numerical tool: the realistic nu-
merical modelling of cracking, the accurate simulation
of the material behaviour and the bearable computa-
tional cost.
Realistic modelling of cracking implies that the loca-
tion and direction of cracks are not spuriously affected
by the used mesh topology. This situation is commonly
referred as mesh-induced directional bias. The depen-
dency of the numerical simulation to the used mesh is
a common limitation of both smeared and discrete crack
approaches [
3]. This defect has triggered the research,
first, on the possible limitations in the numerical dis-
crete problem and, second, on the adequate remedies to
overcome them. The proposed solutions depend on the
perception of the origin of the numerical problem. Some
approaches intervene in the mathematical formulation
of the continuum problem, such as gradient-enhanced
[
4, 5], non-local [6, 7] or micropolar [8, 9] mo dels. Other
formulations aim to overcome the error produced by the
spatial discretization of the continuum domain, such as
2 Savvas Saloustros et al.
mixed formulations [
10], crack-tracking techniques with
[
11, 12, 13, 14] or without [15, 16, 17] enrichment of the
used finite elements or nodes. Finally, another approach
reported in the literature bases on the modification of
the material model [
18]. For a review on the issue the
interested reader is referred to [
19, 20].
Of equal importance to the numerically objective
modelling of cracking, is the accurate simulation of the
material behaviour. The description of the mechani-
cal behaviour in quasi-brittle materials necessitates a
proper failure criterion and a constitutive law with soft-
ening behaviour controlled by the fracture energy. To
account for this behaviour various models have been ap-
plied for the simulation of quasi-brittle materials based
on plasticity [
21, 22, 23], continuum damage models
[
24, 25] or a combination of both [26, 27, 28]. The pre-
vious isotropic formulations have been complemented
by anisotropic models [
29, 30, 31]. In the last years,
an increasing effort has been made to develop numer-
ical methods that consider the interaction of compo-
nents in composite materials (e.g. concrete, masonry).
Such methodologies can b e based on computational ho-
mogenization [
32, 33, 34, 35, 36, 37, 38], or on micro-
modelling techniques, also known as direct numerical
simulations, including all the information about the ma-
terial’s micro-structure [
39, 40].
The micro- and multi-scale modelling approaches,
even though very sophisticated and characterized by a
remarkable level of accuracy, are still hardly applica-
ble to the study of large-scale structural problems ex-
periencing multiple cracking. This is due to the signifi-
cant computational cost required by these models. As a
consequence, the Continuum Finite Element Modelling
(also known as macro-modelling) based on the classi-
cal smeared crack-approach is still on the foreground
of applications in large-scale concrete [
41] and masonry
structures [
42, 43, 44].
The aim of this work is to apply the know-how of
objective mesh-localization approaches to the analysis
of structures made of quasi-brittle materials that exp e-
rience multiple cracking. The adopted numerical tech-
nique is the classical smeared crack approach enhanced
with a local crack-tracking algorithm [
16]. This choice
is justified by the following reasoning. Firstly, tracking
algorithms provide numerical solutions that are free of
mesh-induced directional bias [
45, 16, 17, 46, 47, 48].
Secondly, the particular tracking algorithm can be ef-
fectively combined with constitutive models suitable
for quasi-brittle materials, like anisotropic [
31, 49] or
isotropic ones [
26]. Additionally, an important feature
of the algorithm is the possibility to define a minimum
distance between the discrete cracks. These exclusion
zones surrounding a crack render it possible to con-
sider indirectly the size of the material’s mesoscale (e.g.
masonry units), when the composite is simplified as a
homogenized continuum in the numerical simulation.
Finally, the overhead in the computational cost is lim-
ited, allowing the efficient simulation of multiple crack-
ing problems in large-scale structures.
Up to date approaches using tracking algorithms
have focused mostly on cracks nucleating from the bound-
ary of the domain or from pre-existent flaws and prop-
agating in a single orientation (see for instance [
46, 47,
50, 51, 52, 53]). This confines the application of tracking
algorithms to a very narrow family of structural prob-
lems, in which the origin of cracks is either known prior
to the analysis, or it is imposed by including material
or geometrical perturbations within the analy sed nu-
merical domain. Additionally, structural problems with
cracking starting from an arbitrary point within the
analysed domain and propagating in two opposite ori-
entations have not been addressed.
The main novelty of the present work resides in
the enhancement of a crack-tracking algorithm with
a technique that allows the initiation and propagation
of cracks at any element of the discretized domain. In
particular, the proposed methodology extends the use
of local crack-tracking algorithms to the simulation of
cracking starting from the interior of the discretized
domain and propagating along two opposite orienta-
tions. This novel contribution to the original methodol-
ogy [
16] makes possible the efficient simulation of local-
ized cracks under shear, flexure and traction, without
any a priori definition of their location by the analyst.
Another original aspect of the proposed method is the
identification of the potential crack path prior to the oc-
currence of the damage, differently from the approach
recently proposed in [54].
The paper is structured in the following way. Section
2 presents the underlying continuum damage model.
Section
3 outlines the methodology for the simulation
of formation and propagation of internal cracking with
the use of a crack-tracking algorithm. The proposed
numerical strategy is validated with the simulation of
a single-crack problem with internal fracture (Brazilian
splitting test) and then with the analysis of a complex
multi-crack problem (frame wall with one-door open-
ing tested against horizontal loading), see Section
4.
Finally, conclusions are summarized in Section
5.
2 Continuum Damage model
In this work, a constitutive model based on the con-
tinuum damage mechanics presented in [
25] is used for
the simulation of fracture. The model benefits from the
use of two separate damage scalar indexes to distinguish
Finite Element Modelling of Internal and Multiple Localized Cracks 3
between tensile and compressive damage. This is an im-
portant characteristic for the simulation of degradation
in materials such as concrete and masonry which have
quite different capacity under tension and compression.
Other local models based on the same concept have
been proposed by Lee & Fenves [
27], Comi & Perego
[
55], Wu et al. [28], Pelà et al. [56], Voyiadjis et al. [57],
Mazars et al. [
58]. He et al. [59] and Pereira et al. [60]
have recently proposed nonlocal counterparts of some
of the these models.
The constitutive model establishes on the concept of
the effective stress tensor ¯σ, which basing on the strain
equivalence principle [
61] is related to the strain tensor
ǫ according to the following equation
¯σ = C : ǫ (1)
where C stands for the i sotropic linear-elastic consti-
tutive tensor and (:) the tensor product contracted on
two indices. Aiming to model a different response under
tension and compression, the effective stress is split into
a tensile ¯σ
+
and a compressive part ¯σ
−
. Such tensors
are respectively defined as
¯σ
+
=
3
X
j=1
h¯σ
j
ip
j
⊗ p
j
(2)
¯σ
−
= ¯σ − ¯σ
+
(3)
where ¯σ
j
stands for the j-th principal stress value from
tensor ¯σ and p
j
is the unit vector of the respective prin-
cipal direction j. The symbol h•i denotes the Macaulay
brackets (hxi = x, if x ≥ 0 ,hxi = 0, if x < 0).
Following that, two internal damage variables are
introduced d
+
and d
−
, with the first denoting the ten-
sile damage and the second the compressive one. The
constitutive relation can be expressed then as
σ = (1 −d
+
)¯σ
+
+ (1 − d
−
)¯σ
−
(4)
Loading, unloading and reloading conditions are dis-
tinguished with the use of two scalar positive quantities,
one for tension τ
+
and a second for compression τ
−
,
termed as the equivalent stresses. Their values depend
on the stress tensor and the assumed failure criteria. In
this work, the failure criterion presented in [
26] is used
so that the values of the equivalent stresses for tension
and compression are
τ
+
= H [¯σ
max
]
1
1 − α
h
p
3
¯
J
2
+ α
¯
I
1
+ βh¯σ
max
i
i
f
+
f
−
(5)
τ
−
= H [−¯σ
min
]
1
1 − α
q
3
¯
J
2
−
+ α
¯
I
1
−
+ βh¯σ
max
i
(6)
α =
f
−
b
/f
−
− 1
2
f
−
b
/f
−
− 1
(7)
β = (1 − α)
f
−
f
+
− (1 − α) (8)
In the above equations f
−
b
and f
−
are the biaxial
and uniaxial compressive strengths respectively,
¯
I
1
de-
notes the first invariant of the effective stress tensor,
¯
J
2
the second invariant of the effective deviatoric stress
tensor, while ¯σ
max
and ¯σ
min
stand for the maximum
and minimum principal stress respectively. The failure
surface for the case of the plane stress is shown in Fig-
ure
1. H [ •] denotes the Heaviside step function. Tensile
damage is activated for stress states within the first, sec-
ond and fourth quadrants of the principal stress space,
see Figure
1, while compressive damage for stress states
only within the third quadrant.
σ
1
= σ
3
σ
1
/f
+
σ
3
/f
+
σ
2
= 0
Fig. 1: Adopted damage surface under plane-stress con-
ditions [
26].
The damage criteria are defined then as [
62]
Φ
±
r
±
, τ
±
= τ
±
− r
±
≤ 0 (9)
where r
±
are internal stress-like variables that represent
the current damage thresholds and the respective ex-
pansion of the damage surface. As a consequence, their
initial values are equal to the uniaxial stress under ten-
sion r
+
0
= f
+
and compression r
−
0
= f
−
and thereafter
vary according to
r
±
= max
r
±
0
, max
i∈(0,n)
τ
±
i
(10)
The evolution of the internal damage variables d
±
is defined as [63]
d
±
= 1 −
r
±
0
r
±
exp
(
2H
±
d
r
±
0
− r
±
r
±
0
)
(11)
4 Savvas Saloustros et al.
In the above Equation (
11), H
±
d
≥ 0 stands for the
discrete softening parameter, included to ensure mesh-
size objective results considering the compressive and
tensile fracture energy of the material G
±
f
and the char-
acteristic finite element length. In particular, the spe-
cific dissipated energy (i.e. dissipated energy per unit
of volume) in tension or compression D
±
is regularized
considering the characteristic crack width l
dis
related
to the area (for two-dimensional elements) or volume
(for three-dimensional elements) of each finite element
in the crack band according to the equation [
64]
D
±
l
dis
= G
±
f
(12)
The specific dissipated energy for a damage model
with exponential softening is [
65]
D
±
=
1 +
1
H
±
d
(f
±
)
2
2E
(13)
From Equations (
12) and (13) the previously intro-
duced softening parameter is defined as
H
±
d
=
l
dis
l
mat
− l
dis
(14)
For the current work, two-dimensional linear trian-
gular elements are used with l
dis
=
√
2A, where A de-
notes the area of each element. This selection can be
refined according to the work of Oliver [
66]. The mate-
rial characteristic length l
mat
and the discrete softening
parameter H
mat
depend only on the material properties
according to [
16, 65]
l
mat
=
1
H
mat
=
2EG
±
f
(f
±
)
2
(15)
It is noted that regularized stress versus strain mo d-
els, as the present one, can be shown to be equivalent to
traction versus displacement jumps models such those
used in fracture mechanics [
67], X-FEM and E-FEM
formulations [
68].
3 Modelling of cracking with a local tracking
technique
As discussed in the Introduction, most of the currently
available crack-tracking algorithms focus on the simu-
lation of cracks starting from the boundary of the dis-
cretized domain and propagating towards a single ori-
entation. This drawback limits the application of such
algorithms, which cannot be applied to the analysis of
structures experiencing internal cracking such as shear
cracks. To overcome this limitation, in this work the
local-crack tracking algorithm presented by Cervera et
al. [
16] is adequately enhanced for the simulation of
cracking initiating from internal elements of the mesh
and propagating along two opposite orientations. This
section presents the main features of the algorithm fo-
cusing on its novel contributions. The method is here
applied to constant strain three node elements but it
can be extended to other types of two-dimensional fi-
nite elements (see [17]).
The local crack-tracking algorithm constitutes an
enhancement of the classical smeared crack approach.
The algorithm is called at the beginning of each load
increment of the numerical analysis prior to the eval-
uation of the stresses. Its purpose is to identify and
“label” the elements pertaining to a crack path for the
current increment. For these elements, the evaluation
of the stresses is computed according to the nonlinear
constitutive law defined in Section
2. Contrariwise, the
elements falling out of the crack path will keep their
linear elastic stress-strain relation.
The first procedure carried-out by the crack-tracking
algorithm is the identification of new cracks. A new
crack starts at an element according to the tensile dam-
age criterion defined by Equations (
5) and (9). The
above check is performed throughout the whole dis-
cretized domain and not only at the elements lying over
the boundary as in existing crack-tracking algorithms.
This is necessary to identify and allow the initiation
of internal cracking. The elements satisfying the failure
criterion of Equation (
9) (Φ
+
= 0) are defined as crack
root elements.
The control of the damage dispersion over a small
part of the discretized domain, and thus the simulation
of separate and individual cracks, is possible with the
use of an exclusion radius criterion. This criterion, in-
troduced in [
16], defines as a crack root element the one
with the highest value of the tensile equivalent stress
τ
+
within a radius r
excl
(see Figure
2a). The value of
r
excl
is defined a priori, and may be according to the
mesoscale geometry of the heterogeneous material.
Following the above procedure, the coordinates of
the crack origin are defined and stored. These depend
on the location of the crack root element within the dis-
cretized domain. For the family of internal crack root
elements, the crack origin is defined as the centroid of
the triangular element, i.e. the intersecting point of the
medians (Figure
2a). The same holds for corner ele-
ments, whereas for the remaining boundary elements
the midp oint of the side lying on the mesh boundary is
selected.
The second part of the algorithm, after the defi-
nition of the crack root and crack tip elements, is the
identification of the following elements pertaining to the
crack path. The procedure is different for the two crack
root elements, i.e. boundary and internal ones. For each
Finite Element Modelling of Internal and Multiple Localized Cracks 5
Fig. 2: Simulation of internal cracking with the crack-tracking technique: (a) internal crack root element with the
two opposite vectors of the crack propagation, (b) labelling of the next potential elements towards the first side of
the crack, (c) labelling of the next potential elements towards the second side of the crack and use of the maximum
curvature criterion.
boundary crack root element, a vector is drawn, starting
from the crack origin location, with a direction perpen-
dicular to the one defined by the maximum principal
stress. The intersection of this vector with the neigh-
bouring element defines the exit point and the next po-
tential element of the crack. Similarly, starting from
this point the following next potential elements of the
crack are recognised. The same procedure is followed for
identifying the propagation path of consolidated cracks
from the crack tip elements. In this case the crack ori-
gin point is the exit point of the crack at the previous
cracked element.
Contrarily to cracking starting from the boundary,
other cracks, such as shear ones, initiate from the in-
terior and propagate along two opposite orientations.
To account for this damage typology, the algorithm is
enhanced with a different procedure. Starting from the
crack origin point of the internal crack root element,
two vectors (
−→
V
e,1
) and (
−→
V
e,2
) are defined, having a
direction p erpendicular to that of the maximum prin-
cipal stress but opposite orientations (Figure
2a). Fol-
lowing this, the identification of the potential cracking
elements within the current increment takes part in two
steps. First, the elements pertaining to the path defined
by the orientation of vector (
−→
V
e,1
) are identified start-
ing from the crack origin point and following the same
process as described above for the boundary crack root
elements (Figure
2b). Upon concluding the labelling to-
wards that side of the crack, the elements lying at the
opposite face can be recognised starting again from the
crack origin point of the internal crack root, but using
the orientation of vector (
−→
V
e,2
) (Figure
2c). Figure 3
presents the main steps of the labelling in case of inter-
nal cracks.
The described procedure for the definition of crack
propagation stops on three alternative conditions: (i)
when the next potential element belongs to a different
crack, i.e. when the crack joins another one (meeting
criterion in Figure
3); (ii) wh en a crack reaches the
mesh boundary (boundary criterion in Figure
3); (iii)
when the stress-state of a potential element is lower
than a pre-defined threshold (stress threshold criterion
in Figure
3). This threshold can be conveniently de-
fined in terms of the failure criterion and experience has
demonstrated that labelling can be completed when the
inequality
τ
+
f
+
< 0.75 holds [
16].
The selection of the elements of a crack depending
on the local values of the stresses justifies the “local”
nature of the presented crack-tracking algorithm. Even
if this choice is very convenient in terms of computa-
tional efficiency and for cases of multiple cracks com-
paring to the global crack-tracking algorithms [
16], it
can meet some difficulties under bending stress states,
when the local calculation of the principal stress di-
