ARTICLE IN PRESS
Engineering Fracture Mechanics xxx (2007) xxxxxx
[Link]/locate/engfracmech
Three-dimensional crack initiation, propagation, branching
and junction in non-linear materials by an extended
meshfree method without asymptotic enrichment
Stephane Bordas
a,1
, Timon Rabczuk
b,2
, Goangseup Zi
c,*
Civil Engineering Department, University of Glasgow, Rankine Building, Glasgow G12 8LT, UK
Department of Mechanical Engineering, University of Canterbury, Christchurch, New Zealand
Department of Civil and Environmental Engineering, Korea University, 5 Ga 1, An-Am Dong, Sung-Buk Gu,
Seoul 136-701, Republic of Korea
b
Received 18 December 2006; received in revised form 4 May 2007; accepted 10 May 2007
Abstract
This paper presents a three-dimensional, extrinsically enriched meshfree method for initiation, branching, growth and
coalescence of an arbitrary number of cracks in non-linear solids including large deformations, for statics and dynamics.
The novelty of the methodology is that only an extrinsic discontinuous enrichment and no near-tip enrichment is required.
Instead, a Lagrange multiplier eld is added along the crack front to close the crack. This decreases the computational cost
and removes diculties involved with a branch enrichment. The results are compared to experimental data, and other simulations from the literature to show the robustness and accuracy of the method.
2007 Elsevier Ltd. All rights reserved.
Keywords: Extended element free Galerkin method (XEFG); Discontinuous enrichment; Cohesive cracks; Lagrange multipliers; Dynamic
fracture; Large deformations; High velocity impact
1. Introduction
To model growing cracks in numerical methods such as nite element and meshfree methods, numerous
techniques have been developed and those may be classied into two categories: the strong discontinuity
approach and the weak discontinuity approach. In the former, cracks are modeled as discontinuities in the displacement eld. In the latter, cracks are smeared within a certain region in which the deformation is localized.3
*
Corresponding author. Tel.: +82 2 3290 3324.
E-mail addresses: bordas@[Link] (S. Bordas), [Link]@[Link] (T. Rabczuk), g-zi@[Link] (G. Zi).
1
Tel: +44(0)141 330 4075.
2
Tel: +64 3 364 8836.
3
In this approach, the so-called localization limiter must be introduced to keep the consistency of the energy release rate during the
crack growth.
0013-7944/$ - see front matter 2007 Elsevier Ltd. All rights reserved.
doi:10.1016/[Link].2007.05.010
Please cite this article in press as: Bordas S et al., Three-dimensional crack initiation, propagation, branching ..., Eng
Fract Mech (2007), doi:10.1016/[Link].2007.05.010
ARTICLE IN PRESS
2
S. Bordas et al. / Engineering Fracture Mechanics xxx (2007) xxxxxx
Recently, the strong discontinuity approach is of more interest than the weak discontinuity approach because
one can study the behavior of solids near the sharp crack front with better resolution using the strong discontinuity approach.
Among the many methods of the strong discontinuity approach, the extended nite element method
(XFEM) [1] is one of the most versatile and accurate. This method has been successfully applied to static problems in two and three dimensions, (see e.g. [27]) and to dynamic problems [810] in two dimensions and three
dimensions including contact along the crack faces, and small-scale plasticity [11]. The extended nite element
method is now utilized in industrial settings to assess damage tolerance of complex structures [12,13] and open
source C++ libraries are available, such as [14,15]. Note that the stress eld in the XFEM is not smooth
because the method is an extended version of the standard nite element method as it is implied from its name.
Recently, a posteriori error indicators for the XFEM were devised by Bordas and Duot [1618].
Meshfree methods benet from a higher order of continuity which naturally smoothes the stress eld in the
crack tip region. Fracture simulation work in two and three dimensions using meshfree methods is extensive
[1923], in which discontinuities are treated by the visibility criterion or some modications of it. Other novel
approaches which were able to treat kinked and curved cracks were proposed by Ventura et al. [24]. They
enrich the moving least squares (MLS) basis functions around the crack tip and improve signicantly the convergence behavior.4
Rabczuk et al. proposed the extended element free Galerkin (XEFG) method for cohesive crack initiation,
growth and junction in two and three-dimensional statics and dynamics, but the closure of the crack along the
front is ensured through near-tip enrichment which vanishes along this front [2527] . As noted in [26],
the polar coordinate system along a crack front is not well-dened at the kinks. While the elds spanning
the space of linear elastic fracture mechanics solutions are known for small strains, this is the case neither
in large strain nor for non-linear materials in general. This makes the selection of near-tip enrichment elds
dicult in these circumstances, and calls for an alternate method. In this work, we do not use front enrichment, but only discontinuous enrichment; the crack is closed using a constraint eld enforced by Lagrange
multipliers. Similar accuracy compared to [26] is obtained on the examples treated. Furthermore, we add
and test a crack branching algorithm.
The paper is arranged as follows: in the next section, we will describe the element-free Galerkin method.
Section 3 describes the discontinuous enrichment employed and Section 4 the Lagrange multiplier method
devised to close the crack fronts. The crack tracking algorithm follows, in Section 5 before the presentation
of the governing equations in strong and weak form, in Section 6. The discretized version of the weak form
is then presented in Section 7. We then show, in Section 8, several static and dynamic crack initiation/propagation problems and compare the results to experimental data or/and other numerical results in the literature, and close the work with some conclusions in Section 9.
2. Element-free Galerkin (EFG) approximation
The construction of the MLS approximation employed to describe the displacement eld in the elementfree Galerkin method is standard [19,28,29], and only recalled briey for completeness. In the whole paper,
we denote by X0 the body of interest in the reference conguration, and X X ; Y ; Z 2 X0 R3 a point in
the reference conguration. Let W be the set of all particles in X0. Particles (also called nodes, or points)
are denoted by upper-case indices. Let I be a particle, its approximation radius (smoothing length) is hI,
and its position in X0 is noted XI. For a point X 2 X0, the distance between X and particle I is
rI(X) = kX XIk, where k k is the Euclidian norm of R3 . Each particle is associated with a weight (or kernel)
function W : R ! R. In this paper, we used the cubic B-spline. Continuity in meshfree methods is governed by
the continuity of the kernel function W, in our case, the MLS shape functions derived below are C2-contiIk
nuous, as W. Dening the dimensionless distance rI X rIhX
kXX
, the weight function employed in this
hI
I
work writes
4
In these papers, the MLS basis is enriched extrinsically by adding special functions to the approximation. This is to be opposed to
intrinsic enrichment at the level of the MLS basis.
Please cite this article in press as: Bordas S et al., Three-dimensional crack initiation, propagation, branching ..., Eng
Fract Mech (2007), doi:10.1016/[Link].2007.05.010
ARTICLE IN PRESS
S. Bordas et al. / Engineering Fracture Mechanics xxx (2007) xxxxxx
8I 2 W;
8 X 2 X0
8
2
2
3
>
< 3 4rI 4rI
W rI X 43 4rI 4r2I 43 r3I
>
:
0
for rI 6 12
for 12 < rI 6 12
for rI > 1
The next required ingredient for the MLS approximation is the basis vector denoted by p = [1, X, Y, Z], which
ensures that the MLS approximation can reproduce linear elds in R3 . Denote by IIX the set of NIX particles
whose domain of inuence contain a given point X 2 X0. The MLS shape function UI associated with particle
I evaluated at point X is obtained by the matrix/vector product
1
UI X pT X AX; XI DX; XI
|{z}
|{z}
|{z} |{z}
N IX 1
N IX 4
AX; XI
44
41
pXI p XI W rI X
I2IIX
DX; XI
pXI W rI X
I2IIX
With these denitions, the standard part (C2-continuous and linear-complete) of the MLS approximation
(denoted by a superscript s) of the displacement eld us at positive times t can be written as
X
8X 2 X0 ; 8t > 0; us X; t
UI XusI t
5
I2W
where the unknown coecients usI are functions of time t.
In addition to the fact that the order of continuity can be increased quite easily, meshfree methods have
advantages over nite elements because of their smoothness and nonlocal interpolation character. Better stress
distributions around the crack tip are obtained, leading to a non-oscillatory crack path.
3. Discontinuous enrichment of the displacement eld
To model the discontinuity due to cracks, the displacement approximation is written as the sum of a standard (us) part as in Eq. (5) and an enriched part (ue), discontinuous through the faces of the cracks
8X 2 X0 ; 8t > 0;
uX; t us X; t ue X; t
Dene E as the set of all the cracks in the domain and uea the enriched part (discontinuous part) of the displacement approximation due to crack a. The enriched part of the displacement approximation due to all cracks in
E is the sum of the discontinuous contribution associated with each crack:
X
8X 2 X0 ; 8t > 0; ue X; t
uea X; t
7
a2E
Let us now derive an expression for uea . Let Wa be the set of particles whose domain of inuence is cut by crack
a, and WaI be the enrichment function associated with particle I and crack a, discontinuous through this crack
and dened in detail below. Dene by aaI the additional degrees of freedom for the enrichment WaI . The discontinuous part of the displacement approximation due to crack a 2 E writes
X
8X 2 X0 ; 8t > 0; uea X; t
UI XWaI XaaI t
8
I2Wa
Note that the shape functions UI in Eq. (8) need not be the same as UI in Eq. (5) [3032], it will however be the
case in the numerical examples to follow.
We now dene the enrichment functions WaI . If a domain of inuence is cut by a crack, it is enriched with
the sign function given by
1 for x > 0
signx
9
1 for x < 0
Please cite this article in press as: Bordas S et al., Three-dimensional crack initiation, propagation, branching ..., Eng
Fract Mech (2007), doi:10.1016/[Link].2007.05.010
ARTICLE IN PRESS
4
S. Bordas et al. / Engineering Fracture Mechanics xxx (2007) xxxxxx
f (X ) < 0
f (X ) < 0
f (X ) > 0
f (X ) > 0
Fig. 1. The enriched displacement eld in one dimension enriched (a) by just the sign of the signed distance function f and (b) by the
shifted sign of the signed distance function as in Eq. (11).
Let n be the crack normal and Cac represent the surface of crack a. The choice of the orientation of the crack
normal is completely arbitrary as long as it is consistent throughout the entire computation.
Let Xa be a point on the surface Cac of crack a 2 E and X a point in X0. The sign of quantity n (X Xa)
denes on which side of Cac X is located. The distance from X to Cac is minXC 2Cac kX XC k. The signed distance
f a(X) from X to crack a then writes
8X 2 X0 ; 8Xa 2 Cac ;
f a X signn X Xa min kX XC k
XC 2Cac
10
The enrichment for crack a and particle I writes
8a 2 E; 8I 2 W; 8X 2 X0 ;
WaI X signf a X signf a XI
11
The enriched displacement eld is shown in Fig. 1 for one-dimensional case. Note that the enrichment function was shifted by its nodal values in Eq. (11). The use of the shifted enrichment function makes the enriched
region narrower. The magnitude of the displacement jump is scaled by aaI t in Eq. (8).
4. Enforcing crack closure along their front
Let us now see how the cracks may be closed along their fronts without recourse to near-front enrichment
in meshfree methods. When cracks are not discretized, i.e. when nodes are not present along the crack faces, as
is the case in enriched FEM and meshfree methods, a strategy must be fashioned to close each crack along its
front. Crack closure along the front is a natural outcome of near-front enrichment vanishing along the front,
as in [26] for instance. In the context of the extended nite element method, a methodology was devised to
close the crack without near-front enrichment [4,33] as long as the crack tip is located on element edges. In
C0-FEM-based methods, the shape function associated with each node is only coupled with those of the nodes
contained in its support. Imagine a crack tip in a two-dimensional nite element mesh; see Fig. 2a, where the
crack tip is positioned on the edge connecting nodes A and B. Because the crack must close at its tip, i.e. the
crack opening displacement must be zero, nodes A and B should not be enriched.5 Fig. 2b shows the case of
meshless methods. The domain of inuence for a particle in meshless methods overlaps heavily with that of
other particles. The idea presented in [4,33] would be applicable if and only if the crack tips were located
on the boundary of domains of inuence, but never inside, i.e. if all domains of inuence were either not
cut by the crack, or completely cut by the crack, but never partially cut.
Recently, Zi et al. [34] developed a simple method to close the crack tip without near-tip enrichment. The
method is similar to the technique developed for XFEM by Zi and Belytscko [35]; the domains of inuence
5
If they were enriched with a function discontinuous through the crack interior, the displacement approximation would be
discontinuous along line segment [AB], and the opening of the crack would be non-zero.
Please cite this article in press as: Bordas S et al., Three-dimensional crack initiation, propagation, branching ..., Eng
Fract Mech (2007), doi:10.1016/[Link].2007.05.010
ARTICLE IN PRESS
S. Bordas et al. / Engineering Fracture Mechanics xxx (2007) xxxxxx
Fig. 2. The enrichment for the crack tip by using the step function in (a) the nite element method and (b) meshless methods; solids are
enriched nodes and circles unenriched nodes. (a) Finite element method, (b) meshless method.
containing a crack tip are suitably adjusted so that the crack tip is placed on their edge. However, their
approach is limited for two-dimensional problems only and, as shown in Fig. 3, supports must be articially
modied so as to be entirely split by the crack, and that the crack tip always lie on a support boundary, but
never in its interior. Extending this to three dimensions is not straightforward. The crack front is composed of
a series of segments that would need to be made conform with the boundaries of the spherical supports.
We propose the construction of a Lagrange multiplier eld to close the cracks along their fronts. If only the
sign function enrichment of Eq. (11) is used, the discontinuity Cc extends beyond the crack front, as depicted
Fig. 3. Method proposed by Zi et al. [34] to keep the crack tip on the boundaries of domains of inuence. This methodology is dicult to
scale up to three dimensions where the front is a broken line and the supports are spheres. (a) Standard domains of inuence: the heavy
lines indicate supports which are partially cut by the crack, and would therefore lead to the articial extension of the crack represented in
Fig. 4. (b) Domains of inuence are modied so that the tip is always on their boundary and that there exists no support that is partially
cut by the crack.
Please cite this article in press as: Bordas S et al., Three-dimensional crack initiation, propagation, branching ..., Eng
Fract Mech (2007), doi:10.1016/[Link].2007.05.010
ARTICLE IN PRESS
6
S. Bordas et al. / Engineering Fracture Mechanics xxx (2007) xxxxxx
Fig. 4. When a domain of inuence is cut by a crack, the particle is enriched with the Heaviside (jump) function. If no special care is taken,
the crack articially extends throughout the domain of inuence of the enriched particle, which leads to an inaccurate crack representation.
To avoid this, we add a Lagrange multiplier eld along the crack front, to close it.
in Fig. 4. To correctly model the crack, the discontinuity on the front Cc should vanish. Because the condition
should be satised along a line in 3D space, the Lagrange multiplier eld must be discretized. To avoid introducing additional nodes for the discretization of the Lagrange multipliers, we use the same shape functions as
those for the domain partially cut by the crack. The detailed formulation is given later.
5. Description of cracks
5.1. Geometric description
The cracks are described similarly as in [26]. The crack surfaces are non-planar and represented by the
union of planar segments obtained by slicing the tetrahedral background mesh by the failure planes obtained
through material stability analysis (described in Section 5.2). These crack segments are either triangles or
quadrangles depending on how the tetrahedron is split.
5.2. Initiation and propagation of cracks
We employ the loss of hyperbolicity (in dynamics) and loss of ellipticity (in quasi-statics) criterion for crack
initiation and propagation, as proposed by Belytscko et al. 8. Therefore, a crack is initiated or propagated if
the minimum eigenvalue of the acoustic tensor Q is smaller or equal to zero:
min eigQ 6 0
with Q n A n
12
where n is the normal to the crack surface, A = Ct + r I, r is the stress tensor and Ct is the fourth-order
tangential modulus.
Note that we allow crack branching in a given background cell, as in Fig. 5. The loss of ellipticity/hyperbolicity is checked at the Gau points ahead of the crack fronts; if two signicantly dierent normal directions
(1 and 2 in Fig. 5) are found, the crack is assumed to branch.
5.3. Measure of the crack opening and sliding
Since the standard part, us of the displacement eld is continuous everywhere, the jump in the displacement
eld is governed only by the enriched part ue and is given by
X X
suXt 2
UaI XaaI
13
a2E I2W a X
Please cite this article in press as: Bordas S et al., Three-dimensional crack initiation, propagation, branching ..., Eng
Fract Mech (2007), doi:10.1016/[Link].2007.05.010
ARTICLE IN PRESS
S. Bordas et al. / Engineering Fracture Mechanics xxx (2007) xxxxxx
n1
existing crack
n2
Fig. 5. Crack branching in a given background cell. If more than one Gau point is found to loose stability ahead of an existing crack
front, and the normals obtained by the eigenvalue analysis are signicantly dierent, the crack branches. Cross symbols represent failed
Gau points ahead of the crack front.
Fig. 6. The change of the signed distance function when one crack joins with another, in which the signed distance function of cracks in (a)
and (b) is changed to that in (c) and (d), respectively; the hatched plane represents the side with positive sign.
The normal part dn, i.e. the crack opening and the tangential part su(X)bs, the crack sliding is given by
dn n suXt
dt ksuXt ndn k
14
15
5.4. Tracking the crack path
The algorithmic procedure to track the crack path is detailed in Rabczuk et al. [26]. We extend the algorithm such that it can handle crack branching within the same background cell. Moreover, we allow multiple
cracks to overlap during initiation (at the same time step.6) Especially for problems with excessive cracking,
this facilitates the implementation and also decreases computational cost.
Let us now consider branching and joining cracks as shown in Fig. 6. The discontinuities are located using
the signed-distance function that is called f in the following. Let W1b be the set of nodes whose domain of
6
In contrast to our method in [26] where we did not allow overlapping of simultaneously initiated cracks, see Fig. 7.
Please cite this article in press as: Bordas S et al., Three-dimensional crack initiation, propagation, branching ..., Eng
Fract Mech (2007), doi:10.1016/[Link].2007.05.010
ARTICLE IN PRESS
8
S. Bordas et al. / Engineering Fracture Mechanics xxx (2007) xxxxxx
Fig. 7. Allowing overlap during simultaneous crack initiation. Left: two simultaneously initiating cracks are cut at their intersection, as
in [26]. Right: the two cracks are allowed to overlap, this is the formulation adopted in this paper.
inuence is cut by the discontinuity f1(X) = 0 and W2b the corresponding set for f 2(X) = 0. W3b W1b \ W2b .
By using the signed distance functions of the pre-existing and approaching crack, the signed distance function
of the approaching crack is modied. Consider Fig. 6. Three dierent subdomains have to be considered [5]:
(f 1 < 0, f 2 < 0), (f 1 > 0, f 2 > 0), (f 1 > 0, f 2 < 0) as in Fig. 6b or (f 1 > 0, f 2 < 0), (f 1 > 0, f 2 > 0), (f 1 < 0, f 2 < 0) as
in Fig. 6d. The signed distance function of crack 1 of a point X is then obtained by:
(
f01 X; if f02 X1 f02 X > 0
1
f X
16
f02 X; if f02 X1 f02 X < 0
in which f01 X, f02 X represent the signed distance functions of cracks 1 and 2 without consideration of the
junction, respectively, and X(1) is any point on crack 1. The approximation for a branching crack reads then
nc
X
X
X
a
a
UI XuI
UI X H fI X aI
17
uX
I2WX
a1 I2Wb X
where nc is the number of cracks.
6. Governing equations
6.1. The momentum equation and the boundary conditions
Dening Cc0 as the union of all the crack surfaces, and denoting r0 the gradient operator in the reference
conguration,the strong form of the momentum equation in the total Lagrangian description is given by
.0
18
u r0 P .0 b in X0 n Cc0
with boundary conditions:
uX; t
uX; t on Cu0
n0 PX; t tX; t on Ct0
n0 P n0 P tc0 on Cc0
tc0 tc0 sut on Cc0
19
20
21
22
u is the acceleration, P denotes the nominal stress tensor, b designates the
where .0 is the initial mass density,
body force,
u and t are the prescribed displacement and traction, respectively, n0 is the outward normal to the
domain and Cu0 [ Ct0 [ Cc0 C0 , Cu0 \ Ct0 [ Ct0 \ Cc0 [ Cc0 \ Cu0 ;. Moreover, we assume that the stresses
P are bounded on the crack surface Cc0 . Since the stresses are not well dened along the crack, the crack surface
Cc0 is excluded from the domain X0 which is considered as an open set.
6.2. Cohesive cracks
We exclusively used initially rigid cohesive models as shown in Fig. 8 for one-dimensional stressdisplacement laws. When a potential for the cohesive crack is dened, the unidirectional stressdisplacement relation
Please cite this article in press as: Bordas S et al., Three-dimensional crack initiation, propagation, branching ..., Eng
Fract Mech (2007), doi:10.1016/[Link].2007.05.010
ARTICLE IN PRESS
S. Bordas et al. / Engineering Fracture Mechanics xxx (2007) xxxxxx
Fig. 8. The types of the cohesive laws frequently used in practice; (a) linear (or triangular), (b) bilinear and (c) exponential cohesive laws.
can be extended to general mixed mode problems, as in [3638]. Note that if only one-dimensional cohesive
models for mode-I fracture are employed, traction continuity is in general violated. We only employed cohesive models that fulll the traction continuity condition, i.e. the cohesive traction tc = n r at crack initiation is
compatible with the stress state at crack initiation. We note that the fulllment of traction continuity is a complementary problem and cannot simply be enforced by Lagrange multipliers since the traction cannot exceed
the tensile strength of the material.
Diculties would occur at crack initiation due to the innite slope of the tangent, especially when unloading occurs at an early stage. However, because Gau points for the integration of the traction are never placed
at the crack front, it is not a problem in practice.
7. Discretized equations
7.1. Discretization
The weak form of the momentum equation is given by
dW dW int dW kin dW ext dW coh dW L
23
in which dWint, dWkin, dWext, dWcoh are the parts composing the virtual work of the internal stress, the inertia
force, the external traction and the cohesive traction, respectively; dWL is introduced to close the crack at its
crack front. The four parts of the virtual work write7
Z
T
dW int
r0 du : P dX0
24
X0 nCc0
dW kin
dW ext
Z
Z
X0 nCc0
X0 nCc0
dW coh
.0 du
u dX0
.0 du b dX0
25
Z
du t0 dC0
26
Ct0
sdut s dC0
27
K sut dC0
28
Cc0
dW L d
Cc;ext
0
in which dWL is the general variation with constraint and K is the Lagrange multiplier vector. As the Lagrange
multiplier is dened for Cc;ext
and discretized using the same shape functions as the enrichment, the discretized
0
Lagrange multiplier is given by
7
Recall that Cc0 is the union of all crack surfaces.
Please cite this article in press as: Bordas S et al., Three-dimensional crack initiation, propagation, branching ..., Eng
Fract Mech (2007), doi:10.1016/[Link].2007.05.010
ARTICLE IN PRESS
10
S. Bordas et al. / Engineering Fracture Mechanics xxx (2007) xxxxxx
K 2Uk
29
where k are the coecients of the discretization of the Lagrange multiplier K. Substituting the continuous and
discontinuous displacement elds us and ue in Eqs. (5) and (7), and the crack opening displacement sub in
Eq. (13) into the weak form, we obtain.8
Z
X
T
T
dW int
duI
$0 UI X : P dX
X0 nCc0
I2W
X X
daaK T
a2E K2Wa
dW kin
duTI
XZ
duTI
X X Z
X X
X X
I2W
dW ext
duTI
dW coh 2
X X
dkaI T
a2E I2Wa
X0 nCcb
0
.0 UI X b dX
daTK
a2E K2Wa
dW L 4
T
.0 UK XWaK X UI X dXuI
X0 nCc0
X X
X0 nCca
0
b2E M2Wb
I2W
.0 UI X UK XWaK X dXaa lK
X X Z
daaK T
a2E K2Wa
X0 nCca
0
XZ
daaK T
a2E K2Wa
X X
30
.0 UI X UJ X dXuJ
a2E K2Wa
I2W
$0 UK XWaK X : P dX
T
X0 nCc0
J 2W
I2W
X0 nCca
0
duTI
T
.0 UK XWaK X UM XWbM X dXaLM
Z
UI X t0 dC
31
32
Ct0
UK X tc dC
33
Cca
0
X X Z
b2E J 2Wb
a2E I2Wa
Cca
I2W
kaI T
Cc;ext
0
UaI X UbJ X dC0 abJ
X X Z
b2E J 2Wb
Cc;ext
0
UaI X UbJ X dC0 dabJ
34
Using the fundamental lemma of calculus, we obtain the discretized equations
M
q f ext f coh f int f L
Ga 0
35
36
where M is the consistent mass matrix, q is the vector of generalized parameters, fextR, fint, fcoh are the discrete
external, internal and cohesive force vectors, respectively; fL = kT G in which G 4 Cc; ext UT U dC0 is the force
0
term due to the constraint to close the crack at its front. The expressions for M, q, fext, fcoh and fcoh are given
by
" T
#
Z
T
U0 U0 U0 Ue
M
.0
dX
37
X0 nCc0
Ue T U0 Ue T Ue
Z
Z
0T
B P dX
Be T P dX
38
f int
X0 nCc0
f ext
X0 nCc0
X0 nCc0
.0 U b dX
UTt0 dC
39
Ct0
ca
c
We denote by Cca
0 the surface of crack a in the reference conguration. The union of the C0 for a 2 E forms set C0 .
Please cite this article in press as: Bordas S et al., Three-dimensional crack initiation, propagation, branching ..., Eng
Fract Mech (2007), doi:10.1016/[Link].2007.05.010
ARTICLE IN PRESS
S. Bordas et al. / Engineering Fracture Mechanics xxx (2007) xxxxxx
11
Fig. 9. Numerical integration.
f coh 2
UT tc0 dC
40
Cc0
u
a
s T
T
u uI 8 I 2 W and a aaK
8K 2 Wa ; 8a 2 E
a
0
e
U UI 8I 2 W and U Ua WK 8 K 2 Wa ; 8a 2 E
0
B $0 U and B $0 U
41
42
43
44
7.2. Numerical integration
To integrate the discrete Eqs. (37)(40), tetrahedral background integration cells are constructed. If no
crack crosses an integration cell, standard Gauss quadrature is sucient. Instead of using a sub-tetrahedration
procedure, we use a procedure described in Rabczuk and Areias [27] and modify the quadrature weights
according to Fig. 9. Since in meshfree methods, a larger number of Gauss points is needed for accurate integration, this procedure is preferable since the mapping of quantities and state variables that occurs when a
crack enters a background cell and that generates additional inaccuracies9 can be omitted. The modied quadrature weights are easily obtained by a local Vorono procedure10 on background element level followed
by the determination of the crack intersections and the volumes V+ and V as shown in Fig. 9.
Note that since there is no near-tip enrichment, no special care is required to integrate the non-polynomial,
or high-order polynomial elds that were used in [26], for instance.
8. Numerical examples
8.1. Pull-out test
Consider a pull out test of reinforced concrete as shown in Fig. 10. This example was studied previously by
Gasser and Holzapfel [39] and Areias and Belytschko [40] by the PUFEM and XFEM, respectively. We also
employed symmetry conditions and modeled only one quarter of the specimen. A vertical displacement
boundary condition is applied to pull the reinforcement bar out of the concrete specimen as illustrated in
Fig. 10. We adopted the same constitutive and cohesive model as [39]. At crack initiation, the cohesive traction
are computed from the bulk, i.e. tc = n r to ensure traction continuity. The material parameters are
9
10
Due to the mapping.
That is already implemented in most meshfree codes.
Please cite this article in press as: Bordas S et al., Three-dimensional crack initiation, propagation, branching ..., Eng
Fract Mech (2007), doi:10.1016/[Link].2007.05.010
ARTICLE IN PRESS
12
S. Bordas et al. / Engineering Fracture Mechanics xxx (2007) xxxxxx
168
40
42
168
1060
1200
600
[mm]
1060
1200
plane view
side view
Fig. 10. Test-setup of the pull-out test.
j = 16,670 MPa and m = 12,500 MPa. For the cohesive model, we use, according to [39], t0 = 3 MPa,
a = 11.32 mm1, b = 0.674 and a = 1.
We tested two discretizations starting with 15,000 and 45,000 nodes and rened adaptively [41]. The crack
pattern at dierent load stages is shown in Fig. 11 for the ne discretization and dierent view points. Note
that the surface of the crack is rippling, and is captured nicely. The loaddeection curves for two dierent
renement are shown in Fig. 12 and is similar to the results presented in [26,39,40].
8.2. Chalk under torsion
In this section, we analyze numerically the problem of a circular-cylindrical chalk bar under torsion. The
test setup is shown in Fig. 13. A non-uniform traction boundary condition is applied and a small pertubation
Fig. 11. Crack pattern of the pull-out test.
Please cite this article in press as: Bordas S et al., Three-dimensional crack initiation, propagation, branching ..., Eng
Fract Mech (2007), doi:10.1016/[Link].2007.05.010
ARTICLE IN PRESS
S. Bordas et al. / Engineering Fracture Mechanics xxx (2007) xxxxxx
13
Fig. 12. Loaddeection curves of the pull-out test for dierent renement.
Fig. 13. Set-up of the chalk torsion test.
is introduced on the lateral surface on the midplane normal to the cylinders director. The traction is applied in
the tangential direction to the surface of the chalk and its magnitude t is given by
# p2
t tmax sin
45
2
in which tmax is the maximum traction and # is the angle as shown in Fig. 13. The material properties used are
E = 2 GPa, m = 0.18 and Gf = 50 N/m. The linear cohesive law is adopted for simplicity. tmax was 200 kPa.
If the numerical chalk specimen were perfect, the failure would happen at an arbitrary position in the
longitudinal direction of the chalk. The experimental failure surface of the chalk is compared to the numerical
failure surface in Fig. 14. We note that the numerical results are relatively sensitive to boundary conditions
and to the dimensions of the model specimen; more specically, the ratio between the length and the cross
section is inuential. If the length is large compared to the diameter of the chalk bar, and the tractions applied
uniformly, then the elastic strain energy stored is large, resulting in a relatively straight fracture surface. For
non-uniform boundary conditions, as illustrated in Fig. 14, we obtain a helical fracture surface, comparable to
the experimental failure surface.
8.3. Crack branching
In this section, we examine the performance of the method for a crack branching problem. Consider a rectangular prenotched specimen as shown in Fig. 15. The length of the rectangle is 0.1 m, the width 0.04 m and
depth 0.004 m. Initially, there is a horizontal crack from the left edge to the center of the plate over the entire
thickness. A tensile traction of 1 MPa is applied on the top and bottom edges.
We used Lemaitres damage law [42], loss of hyperbolicity and an exponential decaying cohesive law
(Fig. 8c). The material constants are . = 2450 kg/m3, E = 32 MPa, m = 0.2 and A = 1.0, B = 7300 and
D0 8:5 105 for the Lemaitre model. Two-dimensional computations of this problem was previously
reported in [4345]. A three-dimensional computation was carried out in [46]. Experimental data of this problem is available; see [4749].
The crack pattern is shown in Fig. 16 at dierent time steps. Two crack branches appear, similar to the
previous results in [45,46]. The time history of the crack speed is shown in Fig. 17. The crack starts to
propagate at about 0.012 ms and the crack speed increases almost up to the maximum theoretical value of
Please cite this article in press as: Bordas S et al., Three-dimensional crack initiation, propagation, branching ..., Eng
Fract Mech (2007), doi:10.1016/[Link].2007.05.010
ARTICLE IN PRESS
14
S. Bordas et al. / Engineering Fracture Mechanics xxx (2007) xxxxxx
Fig. 14. Chalk bar under torsion in which the simulation result was depicted using the integration cells.
the Rayleigh wave speed. Then the crack branches and the crack speed quickly decreases again. A similar
result was reported by Belytschko et al. [50]. Following the upper path of the crack11, the crack speed soon
11
The same observations apply for the lower crack path.
Please cite this article in press as: Bordas S et al., Three-dimensional crack initiation, propagation, branching ..., Eng
Fract Mech (2007), doi:10.1016/[Link].2007.05.010
ARTICLE IN PRESS
S. Bordas et al. / Engineering Fracture Mechanics xxx (2007) xxxxxx
15
Fig. 15. A plate with an edge crack loaded by a uniform traction on the top and bottom edges.
Fig. 16. Cracking branching at dierent time steps: (a) t = 30 ms; (b) t = 45 ms; (c) t = 60 ms.
Fig. 17. Crack speed time history for the crack branching problem in which the crack branches as the speed reached the peaks.
Please cite this article in press as: Bordas S et al., Three-dimensional crack initiation, propagation, branching ..., Eng
Fract Mech (2007), doi:10.1016/[Link].2007.05.010
ARTICLE IN PRESS
16
S. Bordas et al. / Engineering Fracture Mechanics xxx (2007) xxxxxx
accelerates and at a certain speed, the crack branches again. Directly after crack branching, the crack speed
decreases monotonically until the crack almost hits the right hand end of the specimen. We note that the crack
speed at the second time of bifurcation is smaller compared to the one when the crack bifurcates rst. The
relation between crack speed and crack branching is still an unsolved problem.
8.4. Taylor bar impact
To test the method for multiple cracks with crack junction, we consider a Taylor bar impact. There are
experimental results available; see Teng et al. [51]. The Taylor bar has a diameter of 6 mm and length
30 mm. We consider an impact velocities of 600 m/s. The failure mechanism is petalling. We have performed
similar computations in [26] but for a much smaller impact velocity where relatively moderate strains occur.
The material is a 2024-T351 aluminium alloy. We use the Johnson Cook model [52] as in the previous section
with Youngs modulus E = 74 GPa, Poisson ratio m = 0.3, density 2700 kg/m3, a reference strain rate of
3.33 104, A = 352 MPa, B = 440 MPa, C = 0.0083, n = 0.42, m = 1, cv = 875 J/kg C, Tr = 296 K, Tm =
775 K and b = 1. We tested two dierent discretizations, with approximately 7000 particles and 22,000
particles. The nal deformation of the Taylor bar is shown in Fig. 18 for both discretizations. As can be seen,
Fig. 18. Final crack pattern of the Taylor bar problem: (a) 22,000 particles; (b) 7000 particles.
Please cite this article in press as: Bordas S et al., Three-dimensional crack initiation, propagation, branching ..., Eng
Fract Mech (2007), doi:10.1016/[Link].2007.05.010
ARTICLE IN PRESS
S. Bordas et al. / Engineering Fracture Mechanics xxx (2007) xxxxxx
17
multiple cracking occurs including crack junctions. A similar failure mode was observed in Teng et al. [51].
There are basically four major cracks that cause the petalling. Our failure mode is slightly too ductile, which
is likely related to diculties in the discontinuous bifurcation analysis since the Johnson Cook material does
not lose stability easily. The results look almost identical for both discretizations, and agree well with the
experimental data in [51].
9. Conclusions
This paper presented a three-dimensional adaptive meshfree method for fracture in statics and dynamics.
The initiation, growth, coalescence and branching of an arbitrary number of cracks is handled simply and
eectively.
The discontinuities are introduced through extrinsic enrichment of the moving least squares basis, but no
near-front enrichment is required. To close the cracks, a Lagrange multiplier eld is instead added along the
front of the cracks. Geometrically, the cracks are non-planar surfaces composed of triangular and quadrangular planar sections obtained by cutting the tetrahedral background integration cells by planes whose normals are provided by a material stability analysis.
The results show accurate simulations of large deformation failure problems including fragmentation,
where the exibility of the meshfree method coupled with the ecient crack interaction procedure is most
clear. The simulations agree well with experimental results available in the literature.
Acknowledgement
G. Zi would like to thank the Agency for Defense Development (ADD) for grant ADD-06-05-06.
References
[1] Belytschko T, Black T. Elastic crack growth in nite elements with minimal remeshing. Int J Numer Meth Engng 1999;45(5):60120.
[2] Moes N, Dolbow J, Belytschko T. A nite element method for crack growth without remeshing. Int J Numer Meth Engng
1999;46(1):13350.
[3] Moes N, Gravouil A, Belytschko T. Non-planar 3-D crack growth by the extended nite element method and level sets, part I:
Mechanical model. Int J Numer Meth Engng 2002;53(11):254968.
[4] Zi G, Belytschko T. New crack-tip elements for xfem and applications to cohesive cracks. Int J Numer Meth Engng
2003;57(15):222140.
[5] Zi G, Song J-H, Budyn E, Lee S-H, Belytschko T. A method for growing multiple cracks without remeshing and its application to
fatigue crack growth. Model Simulat Mater Sci Engng 2004;12(1):90115.
[6] Gravouil A, Moes N, Belytschko T. Non-planar 3D crack growth by the extended nite element and level sets part II: level set
update. Int J Numer Meth Engng 2002;53:256986.
[7] Areias PMA, Belytschko T. Non-linear analysis of shells with arbitrary evolving cracks using XFEM. Int J Numer Meth Engng
2005;62:384415.
[8] Belytschko T, Chen H, Xu J, Zi G. Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment.
Int J Numer Meth Engng 2003;58(12):1873905.
[9] Zi G, Chen H, Xu J, Belytschko T. The extended nite element method for dynamic fractures. Shock Vibrat 2005;12(1):923.
[10] Song J-H, Areias PMA, Belytschko T. A method for dynamic crack and shear band propagation with phantom nodes. Int J Numer
Meth Engng 2006;67:86893.
[11] Rethore J, Gravouil A, Combescure A. An energy-conserving scheme for dynamic crack growth using the eXtended fnite element
method. Int J Numer Meth Engng 2005;63:63159.
[12] Bordas S, Moran B. Enriched nite elements and level sets for damage tolerance assessment of complex structures. Engng Fract Mech
2006;73:1176201.
[13] Bordas S, Conley JG, Moran B, Gray J, Nichols E. A simulation-based design paradigm for complex cast components. Engng
Comput 2007;23(1):2537.
[14] Bordas S, Nguyen VP, Dunant C, Nguyen-Dang H, Guidoum A. An extended nite element library. Int J Numer Meth Engng, in
press, doi:10.1002/nme.1966.
[15] Dunant C, Nguyen P, Belgasmia M, Bordas S, Guidoum A, Nguyen-Dang H. Architecture trade-os of including a mesher in an
object-oriented extended nite element code. Eur J Comput Mech 2007;(16):23758.
[16] Bordas S, Duot M, Le P. A simple a posteriori error estimator for the extended nite element method, Commun Numer Meth
Engng, in press. doi:10.1002/cnm.1001.
Please cite this article in press as: Bordas S et al., Three-dimensional crack initiation, propagation, branching ..., Eng
Fract Mech (2007), doi:10.1016/[Link].2007.05.010
ARTICLE IN PRESS
18
S. Bordas et al. / Engineering Fracture Mechanics xxx (2007) xxxxxx
[17] Bordas S, Duot M. Derivative recovery and a posteriori error estimation in extended nite element methods, Computer Meth Appl
Mech Engng, in press. doi:10.1016/[Link].2007.03.011.
[18] Duot M, Bordas S. An extended global recovery procedure for a posteriori error estimation in extended nite element methods, Int J
Numer Meth Engng, in press.
[19] Belytschko T, Tabbara M. Dynamic fracture using element-free Galerkin methods. Int J Numer Meth Engng 1996;39(6):92338.
[20] Lu Y, Belytschko T, Tabbara M. Element-free Galerkin method for wave-propagation and dynamic fracture. Comput Meth Appl
Mech Engng 1995;126(12):13153.
[21] Belytschko T, Lu Y. Element-free Galerkin methods for static and dynamic fracture. Int J Solids Struct 1995;32:254770.
[22] Belytschko T, Lu Y, Gu L. Crack propagation by element-free Galerkin methods. Engng Fract Mech 1995;51(2):295315.
[23] Krysl P, Belytschko T. The element free Galerkin method for dynamic propagation of arbitrary 3-D cracks. Int J Numer Meth Engng
1999;44(6):767800.
[24] Ventura G, Xu J, Belytschko T. A vector level set method and new discontinuity approximations for crack growth by EFG. Int J
Numer Meth Engng 2002;54(6):92344.
[25] Rabczuk T, Zi G. A meshfree method based on the local partition of unity for cohesive cracks. Comput Mech 2007;39(6):74360.
[26] Rabczuk T, Bordas S, Zi G. A three-dimensional meshfree method for continuous multiple-crack initiation, propagation and junction
in statics and dynamics, Comput Mech, 2006, in online, doi:10.1007/s00466-006-0122-1.
[27] Rabczuk T, Areias P. A meshfree thin shell for arbitrary evolving cracks based on an external enrichment. Comput Model Engng Sci
2006;16(2):11530.
[28] Belytschko T, Lu Y, Gu L. Element-free Galerkin methods. Int J Numer Meth Engng 1994;37:22956.
[29] Fleming M, Chu YA, Moran B, Belytschko T. Enriched element-free Galerkin methods for crack tip elds. Int J Numer Meth Engng
1997;40:1483504.
[30] Strouboulis T, Copps K, Babuska I. The generalized nite element method: an example of its implementation and illustration of its
performance. Int J Numer Meth Engng 2000;47(8):140117.
[31] Chessa J, Wang H, Belytschko T. On the construction of blending elements for local partition of unity enriched nite elements. Int J
Numer Meth Engng 2003;57(7):101538.
[32] Stazi F, Budyn E, Chessa J, Belytschko T. XFEM for fracture mechanics with quadratic elements. Comput Mech 2003;31:3848.
[33] Dolbow J, Devan A. Enrichment of enhanced assumed strain approximation for representing strong discontinuities: Addressing
volumetric incompressibility and the discontinuous path test. Int J Numer Meth Engng 2003;59:4767.
[34] Zi G, Rabczuk T, Wall W. Extended meshfree methods without the branch enrichment for cohesive cracks. Comput Mech
2007;40(2):36782.
[35] Zi G, Belytschko T. New crack-tip elements for XFEM and applications to cohesive cracks. Int J Numer Meth Engng
2003;57:222140.
[36] Ortiz M, Pandol A. Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis. Int J
Numer Meth Engng 1999;44:126782.
[37] Bazant ZP, Caner FC. Microplane model M5 with kinematic and static constraints for concrete fracture and anelasticity I: Theory.
ASCE J Engng Mech 2005;131(1):3140.
[38] Camacho GT, Ortiz M. Computational modeling of impact damage in brittle materials. Int J Solids Struct 1996;33:2899938.
[39] Gasser T, Holzapfel G. Modeling 3D crack propagation in unreinforced concrete using pufem. Comput Meth Appl Mech Engng
2005;194:285996.
[40] Areias P, Belytschko T. Analysis of three-dimensional crack initiation and propagation using the extended nite element method. Int
J Numer Meth Engng 2005;63:76088.
[41] Rabczuk T, Belytschko T. Adaptivity for structured meshfree particle methods in 2D and 3D. Int J Numer Meth Engng
2005;63(11):155982.
[42] Lemaitre J, Chaboche JL. Mechanics of solid materials. Cambridge: Cambridge University Press; 1990.
[43] Xu X-P, Needleman A. Numerical simulations of fast crack growth in brittle solids. J Mech Phys Solids 1994;42:1397434.
[44] Falk M, Needleman A, Rice J. A critical evaluation of cohesive zone models of dynamic fracture. J Phys IV 2001;11(PR5):4350.
[45] Rabczuk T, Belytschko T. Cracking particles: a simplied meshfree method for arbitrary evolving cracks. Int J Numer Meth Engng
2004;61(13):231643.
[46] Rabczuk T, Belytschko T. A three dimensional large deformation meshfree method for arbitrary evolving cracks. Comput Meth Appl
Mech Engng 2007;196(2930):277799.
[47] Ravi-Chandar K. Dynamic fracture of nominally brittle materials. Int J Fract 1998;90(12):83102.
[48] Sharon E, Fineberg J. Microbranching instability and the dynamic fracture of brittle materials. Phys Rev B 1996;54(10):712839.
[49] Fineberg J, Sharon E, Cohen G. Crack front waves in dynamic fracture. Int J Fract 2003;121(12):5569.
[50] Belytschko T, Chen H, Xu J, Zi G. Dynamic crack propagation based on loss of hyperbolicity with a new discontinuous enrichment.
Int J Numer Meth Engng 2003;58(12):1873905.
[51] Teng X, Wierzbicki T, Hiermaier S, Rohr I. Numerical prediction of fracture in the Taylor test. Int J Solids Struct 2005;42:292948.
[52] Johnson G, Cook W. A constitutive model and data for metals subjected to large strains, high strain rates, and high temperatures.
In: Proceedings of the 7th international symposium on ballistics; 1983.
Please cite this article in press as: Bordas S et al., Three-dimensional crack initiation, propagation, branching ..., Eng
Fract Mech (2007), doi:10.1016/[Link].2007.05.010