Computational Materials Science 104 (2015) 98–107
Contents lists available at ScienceDirect
Computational Materials Science
journal homepage: www.elsevier.com/locate/commatsci
XFEM–dislocation dynamics multi-scale modeling of plasticity
and fracture
Amirreza Keyhani a, Mohsen Goudarzi a, Soheil Mohammadi a,⇑, Reza Roumina b
a
b
High Performance Computing Lab (HPC), School of Civil Engineering, College of Engineering, University of Tehran, Tehran, Iran
School of Metallurgical and Materials Engineering, College of Engineering, University of Tehran, Tehran, Iran
a r t i c l e
i n f o
Article history:
Received 2 November 2014
Received in revised form 21 February 2015
Accepted 24 March 2015
Keywords:
XFEM–DD
Extended multi-scale
Plasticity
Cohesive crack
Precipitate
a b s t r a c t
An extended three-dimensional multi-scale framework is presented to model plasticity in precipitate
containing materials in the presence of cracks. The present framework is constructed in two scales. At
the micro scale, plasticity is computed via the line dislocation dynamics (DD) methodology in which
the penetrable and impenetrable precipitates are modeled. At the macro scale, application of the
extended ﬁnite element method (XFEM) allows for accurate analysis of mixed-mode cohesive crack propagation without the need for any remeshing procedure. This is a signiﬁcant improvement especially when
two different scales are involved. While the former simulations have so far been limited to the study of
only mode-I crack propagation to avoid the expensive remeshing procedure, in the present work, the
effects of loading rates and precipitate density on general crack propagations are studied by utilizing
the XFEM–DD multi-scale framework to illustrate the efﬁciency and capability of the proposed approach.
Ó 2015 Elsevier B.V. All rights reserved.
1. Introduction
The fracture behavior of materials is mainly controlled by two
mechanisms, separation of atomic bonds in brittle materials and
plastic deformation in ductile specimens. Dislocations motion is
directly responsible for plastic dissipation near the crack tip region
in a ductile material. The dissipative mechanisms of dislocations
are only rate controlling if there is a possibility of spreading plasticity by dislocations motion around the crack tip before the crack
propagates. Interfacial properties of fracture surfaces can also have
a controlling effect on cracking behavior [1]. In addition, the
microstructural anisotropy near the crack tip region affects the
crack propagation. For instance, random distribution of precipitates or arrangement of dislocations with respect to the crack
direction may create a mixed-mode crack propagation state even
if the crack is loaded purely in mode I. Therefore, a multi-scale
analysis, which is governed by the micro scale dislocations motion
within the ﬁnite element modeling of cracking, is of great importance as it can provide more insight into the role of plasticity
and microstructural features on fracture properties of crystalline
materials.
⇑ Corresponding author. Tel.: +98 21 6111 2258; fax: +98 21 6640 3808.
E-mail addresses: [email protected] (A. Keyhani), [email protected]
(M. Goudarzi), [email protected] (S. Mohammadi), [email protected] (R. Roumina).
http://dx.doi.org/10.1016/j.commatsci.2015.03.032
0927-0256/Ó 2015 Elsevier B.V. All rights reserved.
Dislocation dynamics (DD), basically designed for modeling
evolutions of dislocations [2–4], has been widely used to model
plasticity-related phenomena at the micro scale [5–7]. Based on
the fact that in dislocation dynamics computations, the elastic
ﬁelds of dislocations are usually considered in an inﬁnite homogenous domain, DD has been coupled with the ﬁnite element method
(FEM) or the boundary element method (BEM) in order to treat
boundary conditions properly [8–12]. The combined framework
has been used to analyze plasticity in multi-scale frameworks
[11,13,14] besides its development to study the crack tip plasticity
and its effects on crack propagation [15–24]. For a review of
advances in discrete dislocation modeling, see [25,26].
Several authors have adopted a dislocation-independent cohesive crack model to analyze crack plasticity problems. Cleveringa
et al. [16,17] and Deshpande et al. [18–20] utilized a cohesive crack
model based on the universal binding law of Rose et al. [27] for
analyzing the stationary crack tip plasticity, mode-I crack propagation and fatigue crack growth by dislocation dynamics in a single crystal. More recently, Olarnrithinun et al. [21] and Shishvan
and van der Giessen [23,24] utilized a similar cohesive law to predict the plastic response near the stationary crack tip in anisotropic
crystals and to simulate mode-I crack growth.
Former studies on multi-scale modeling of crack tip plasticity
and crack propagation were mainly two dimensional and limited
to only mode-I fracture due to some computational restrictions.
For instance, a number of studies used symmetry to model half
99
A. Keyhani et al. / Computational Materials Science 104 (2015) 98–107
of a cohesive crack and adopted axial springs to represent cohesive
stresses [15–24]. Such modeling procedures did not allow for simulation of mixed-mode crack propagation because of limitations of
the conventional ﬁnite element method. Therefore, mixed-mode
crack propagations were practically avoided by symmetry
assumptions.
Modeling mixed-mode crack propagation by the ﬁnite element
method requires complicated remeshing algorithms for arbitrary
orientation of crack segments [28,29], which becomes extremely
difﬁcult and expensive for coupled two-scale problems. Among
the available remedies to avoid these restrictions, meshfree methods, the boundary element method and even the isogeometric
analysis may all seem sophisticated methodologies that can
improve the simulations to some degrees [30–33]. However, the
recent trends toward the extended ﬁnite element method
(XFEM) have well proved its efﬁciency and robustness for modeling of arbitrary discontinuities [34].
Since the creation of XFEM [35–38], it has successfully been
applied to cohesive zone modeling. Moës and Belytschko [39] and
Wells and Sluys [40] were the ﬁrst to adopt XFEM to address the
mixed-mode cohesive crack propagation. Zi and Belytschko [41]
proposed special tip elements without branch enrichments for
cohesive cracks. The effect of higher order tip enrichment functions
based on the existing analytical solutions around the tip of the
cohesive cracks was studied by Cox [42], while the use of energy
based criteria for accurate computation of crack propagation angle
was examined numerically by Dumstorff and Meschke [43].
In addition to handling crack propagation, XFEM was adopted to
model dislocations in two and three dimensions as well as in thin
shell carbon nanotubes [44–46]. Coupling XFEM with the atomistic
scale allowed for efﬁcient simulation of dislocations and their
motion [47,48], and XFEM–DD was developed for direct dynamic
dislocation simulations in electromechanical ﬁelds [49]. For a comprehensive review of development of XFEM, see [50,51].
In the present research, an extended multi-scale framework is
presented to model mixed-mode crack propagation in precipitate
contained crystalline materials. At the micro scale, the line dislocation dynamics code DDLab [52,53] is modiﬁed to model dislocation–precipitate interactions [54]. The macro scale cracked
domain is modeled by utilizing XFEM, which allows mixed-mode
crack propagation without the need for any remeshing. The two
adopted methods, XFEM and the modiﬁed dislocation dynamics
are consistent in the way that both are designed to avoid complex
mesh and remeshing procedures. As a result, in the whole process
of multi-scale analysis, the mesh remains independent of the
geometry of cracks and precipitates. This simpliﬁes solution of
large systems with arbitrary distribution of precipitates on a very
simple ﬁnite element mesh, with less DOFs and low computational
costs. The present XFEM–DD framework is an extension of FEM–
DD framework [14], in which dislocation movements are handled
by the ordinary three dimensional dislocation dynamics approach
and the mixed-mode crack propagation is modeled through XFEM.
2. Governing equations
2.1. Modiﬁed line dislocation dynamics formulation
In the line dislocation dynamics analysis, a dislocation curve is
discretized into straight segments; each deﬁned by its two end
nodes. As a result, the overall dislocation movements are determined by the movement of segment nodes [52]. The movement
of each dislocation segment contributes to the total generated
macroscopic plastic strain ep , [52]
ep ¼
N
X
li v i
ðni bi þ bi ni Þ
2V RVE
i¼1
ð1Þ
with summation over the number of segments N. li , v i , ni and bi are
the segment length, magnitude of glide velocity, unit normal vector
to the slip plane, and the Burgers vector of dislocation, respectively.
V RVE is the volume of the ﬁnite element in which the dislocation
glide occurs.
Due to signiﬁcant importance of dislocation–precipitate interaction in computational materials modeling [55], several methods
have been developed to model precipitates in dislocation dynamics
[56–59]. An efﬁcient computational technique to model both penetrable and impenetrable precipitates has recently been presented
by Keyhani et al. [54]. In order to model impenetrable precipitates
by utilizing this technique, it is considered that a dislocation node
which is closer than a speciﬁc distance to a precipitate is locked.
Consequently, the problem of dislocation–precipitate interaction
turns into the Frank–Read (F–R) mechanism, Fig. 1. The critical
stress of two mechanisms must be equal to obtain the same result.
Therefore, it is considered that the dislocation rounds the precipitate with a different diameter than its original size. This modeling
diameter ðDm Þ is deﬁned by
Dm ¼ L þ D 2pLb½lnðD1 =r 0 Þ
1
ð2Þ
where L is the precipitate spacing, D is the precipitate diameter, r 0 is
the core radius of a dislocation and equals to the magnitude
1
of the Burgers vector [52], D1 ¼ ðD1 þ L1 Þ , and b is a constant
representing the anisotropy of dislocation line tension and considered to be 0.81 [55]. A dislocation node is locked when it reaches
a distance of modeling diameter to the precipitate. At each step,
the local shear stress related to the local curvature at the locked
node is compared with the precipitate resistance. If the local shear
stress exceeds the precipitate resistance, it is released.
2.2. XFEM–DD formulation
In order to model crack propagation in the presence of dislocations, the problem can be separated into two parts by applying the
superposition principle, as depicted in Fig. 2 [8–11]. The ﬁrst part is
handled by the dislocation dynamics theory, where dislocations
and their interactions with precipitates in an inﬁnite domain are
computed (Fig. 2b). In the second part, the boundary conditions
are satisﬁed by introducing dislocation image stresses within an
XFEM framework (Fig. 2c). Similar to [16–24], the universal binding law of Rose et al. [27] is adopted to describe the constitutive
response of cohesive crack.
A ﬁnite domain X, surrounded by the boundary C, is deﬁned
with the external tractions t and the prescribed displacement u
on part of C. Boundaries of a cohesive crack are decomposed into
two parts: the cohesive region at the crack tip, denoted by Cc ,
and the crack surfaces which are included in C. The traction t1 is
induced from the dislocations elastic ﬁeld on the boundary C of
the inﬁnite domain (Fig. 2b). The stress ﬁeld resulting from the
image traction t1 on the boundary C of the XFEM framework
along with the dislocations elastic ﬁeld in the inﬁnite domain allow
for the complete account of dislocation–surface and dislocation–
crack interactions.
The governing boundary value equation for an isotropic elastic
body under isothermal conditions and assuming the small deformation regime can be expressed as
rrþb¼0
ð3Þ
where r is the Cauchy stress tensor and b represents the body force,
including the forces due to the presence of dislocations (f B ) and the
movement of dislocations ðf P Þ. Essential and natural boundary conditions of the XFEM framework are
on a part of C
u¼u
ð4Þ
100
A. Keyhani et al. / Computational Materials Science 104 (2015) 98–107
Fig. 1. A dislocation line interacting with an array of precipitates. Nodes positioned closer than a speciﬁc distance to the precipitate are locked. A dislocation segment
between two precipitates with the spacing L acts as an F–R source with length Lf [54].
Fig. 2. Coupled XFEM–DD boundary value problem using the superposition principle, (a) the main problem, (b) the dislocation problem in the inﬁnite domain, and (c) the
continuum scale modeled with XFEM.
r nC ¼ t t1 on C
ð5Þ
r nCc ¼ tc nCc on Cc
ð6Þ
nC and nCc are the unit normal vectors to C and Cc , respectively. The
traction tc deﬁnes the transferring cohesive traction at the cohesive
region at the tip of a crack ðCc Þ, which is extracted from the governing interfacial law in terms of the crack opening displacement. The
ﬁnal weak form of the governing equations can be written as
Z
T
dðLu uÞ rdX X
Z
duT bdX Z
duT ðt t1 ÞdC þ
C
X
Z
T
dsut tc dCc ¼ 0
Cc
ð7Þ
where dðuÞ is the displacement variation, dsut ¼ duþ du and +/
signs correspond to upper and lower edges of the crack, respectively. Lu is deﬁned as
2
@
@x
0
6
Lu ¼ 4 0
3
ð8Þ
@
@x
ð9Þ
I¼1
where Nstd
are the basic ﬁnite element shape functions, ustd
repreI
I
are
the
enrichment
sent the standard nodal point unknowns, uenr
I
DOFs, and HðxÞ is the Heaviside enrichment function,
1 if ðx xcr Þ n P 0
0
ð11Þ
where
f coh ¼
R
R
Cc tc Nu dCc
u dC
f ext ¼ C tN
R
f B ¼ X SD Bu dX
R
f 1 ¼ C t1 Bu dC
R
f p ¼ X Dep Bu dX
u ¼ ½ustd
ð12Þ
uenr Nstd H
Nu ¼ ½N
n
n
X
X
std
enr
uh ðxÞ ¼
Nstd
Nstd
I ðxÞuI þ
I ðxÞHðxÞuI
HðxÞ ¼
rBu dX þ f coh ¼ f ext þ f B þ f 1 þ f P
X
std
The XFEM approximation is now adopted to accurately account for
modeling of arbitrary discontinuity paths,
I¼1
Z
with
@ 7
@y 5
@
@y
ﬁeld. Assuming the inﬁnitesimal strain theory similar to [16], the
ﬁnal discretized form of the XFEM–DD formulation, which accounts
for the full crack–dislocation coupling based on the superposition
strategy of Fig. 2 and the corresponding boundary conditions, is
written based on [14],
otherwise
ð10Þ
where xcr is the nearest point to a point x on the crack surface and n
is the normal vector to crack surface at the point xcr . The second part
of Eq. (9) represents the discontinuous part of the displacements
ð13Þ
Bu ¼ rNu
D is the elastic constitutive tensor, which is assumed to be isotropic,
as widely adopted in dislocation dynamics analysis of plasticity in
single crystals [16–20,22]. For instance, despite the fact that the
elastic response of Cu crystal is anisotropic at the micro scale, several studies have reported low sensitivity of the results with respect
to the elastic anisotropy [23,24,60]. In Eq. (13), u represents the
nodal displacement vector. Superscripts std and enr stand for the
standard and enrichment parts of the solution, respectively, and
SD is the homogenized stress due to the presence of dislocations
in the ﬁnite volume of V RVE ,
SD ¼
1
V RVE
Z
XRVE
rD dX
ð14Þ
101
A. Keyhani et al. / Computational Materials Science 104 (2015) 98–107
1.5
1
1
0.5
[0 1 0] μm
[0 0 1] μm
0.5
0
−0.5
0
−0.5
−1
−1
0
−1
1
−1.5 −1
[1 0 0] μm
−0.5 0
1
0.5
[0 1 0] μm
1.5
−1.5
−1.5
−1
−0.5
0
0.5
1
1.5
[1 0 0] μm
Fig. 3. The initial geometry of the problem: The red lines and the blue surface present the dislocation lines (F–R sources) and the penny crack, respectively. (For interpretation
of the references to colour in this ﬁgure legend, the reader is referred to the web version of this article.)
Fig. 4. (a) The dislocation glide plane, (b) two-dimensional XFEM–DD model, red squares and blue circles represent the heaviside enriched nodes and the domain where f B is
applied. The image force due to the crack free surfaces, f 1 , is applied on the upper crack edge, (c) Typical sub triangulation integration technique, where crosses present the
integration points of triangular sub elements for evaluation of rD . (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version
of this article.)
where rD is the stress ﬁeld induced from the presence of dislocations. As a dislocation line is discretized into straight segments,
rD is evaluated by summation of stress ﬁelds of all segments. To
−5
0
x 10
−0.5
K (N.μm3/2)
−1
−1.5
−2
K Eq. (22.1)
I
−2.5
K XFEM−DD
I
KII Eq. (22.2)
−3
K XFEM−DD
II
−3.5
0
500
1000
1500
2000
2500
3000
3500
r/b
Fig. 5. Comparision of numerical and analytical values of shielding stress intensity
factor for h ¼ p=4.
evaluate the stress ﬁeld generated by a dislocation segment, the
numerical approach used in the DDLab code is adopted. For further
details see [61]. The well-developed sub-triangulation technique,
which has been widely used in XFEM [50,51], is adopted to evaluate
the integral accurately (see typical Fig. 4c). The integrand rD in Eq.
(14) is evaluated on the integration points of the triangular or tetrahedral sub-elements (for two or three dimensional problems,
respectively) of the ﬁnite elements which contain a dislocation.
The sub-triangulation technique ensures that the integration points
are not located on the position of the singular points.
The last three terms of the right hand side of Eq. (11) are related
to the line dislocation dynamics. They are the equivalent force vectors due to boundary conditions and movements of dislocations. f B
is the body force vector due to dislocations long range interaction
in an inﬁnite homogenous domain. Therefore, f 1 is included to
represent the effect of ﬁnite boundary conditions. The last term,
f p , is the body force vector from the plastic strains, which is caused
by the dislocation motions.
Eq. (11) is then linearized using the Newton–Raphson method.
The ﬁnal residual equations for iteration i can then be written as
i
i
i
i
i
Ri ¼ Kui þ f coh f ext f 1 f B f P
ð15Þ
R
where K ¼ X BTu DBu dX. The following incremental algorithm is
iterated in each step of load increment until the solution converges,
Riþ1 ¼ Ri þ J½du
iþ1
¼0
ð16Þ
102
A. Keyhani et al. / Computational Materials Science 104 (2015) 98–107
−5
0
x 10
J¼ Kþ
K (N.μm3/2 )
I
−1
−1.5
2 0
sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ13
2
1
K
KI
I
h ¼ 2 arctan 4 @ þ 8A5
4 K II
K II
−2
−2.5
Eq. (22.1)
Mesh 101x101
Mesh 125x125
Mesh 151x151
−3
0
500
1000
1500
2000
2500
3000
3. Results and discussions
3500
Fig. 6. Mesh sensitivity analysis of the results.
iþ1
¼ J1 Ri
½uiþ1 ¼ ½ui þ ½du
ð17Þ
iþ1
ð18Þ
where J is the Jacobian matrix (derivative of the residual vector),
J¼Kþ
@f coh @f 1 @f B @f P
@u
@u
@u @u
ð19Þ
Variation of forces arising from dislocations depends on dislocation
movements. Nevertheless, due to the fact that dislocation movements in each step are limited to acutely small distances (less than
one percent of element dimension), therefore, without the loss of
accuracy, the equivalent forces due to dislocations are assumed to
be constant in each step. As a result,
Jacobian is further simpliﬁed to
@f 1
@u
B
P
¼ @f
¼ @f
¼ 0 and the
@u
@u
Various simulations are performed to study the crack tip plasticity and its effects on mixed-mode crack propagation by utilizing
the present framework. In Section 3.1, crack shielding due to the
presence of dislocations is compared with the analytical solutions.
The effects of loading ratio and precipitates density on crack propagation are studied in Sections 3.2 and 3.3, respectively. Section 3.4
is dedicated to a special case where the mixed-mode crack propagation occurs while the crack is loaded in mode I. The elastic
behavior of Cu crystal is assumed to be isotropic with the shear
modulus of G ¼ 54:6 GPa, the Poisson’s ratio t ¼ 0:324, and the
magnitude of the Burgers vector is b ¼ 0:255 nm. Except for
Section 3.1, a linear decaying cohesive law with the tensile
strength f t ¼ 0:02G and the fracture energy Gf ¼ 5 N=mm is
adopted. A special case of penny shape crack with isotropic growth
is considered in the present simulations (i.e. the crack remains
penny shape after propagation). The initial diameter of crack is
1 lm and the crack is loaded perpendicular to the crack plane
(mode-I loading). Simulations are continued until the crack diameter reaches to about 2 lm. The initial geometry of the problem, the
penny shape crack and the surrounding F–R sources are presented
in Fig. 3. The red lines present dislocations and the blue surface
δε/δt=4.407e7
δε/δt=6.414e7
1
1
0
0
−1
−1
−1
−1
0
1
−1
0
0
1
1
δε/δt=7.242e7
1
0
0
−1
−1
−1
−1
1
−1
0
−1
0
1
δε/δt=1.051e8
1
0
ð21Þ
where K I and K II are computed by the interaction integral [39].
r/b
½du
ð20Þ
Crack propagation is activated whenever the value of normal stress
at the tip of discontinuity exceeds the tensile strength of materials.
Then, a new segment with a predeﬁned length is added to the existing crack length. The direction of crack propagation is governed by
the maximum hoop stress criterion in terms of the mode I and II
stress intensity factors,
−0.5
−3.5
@f coh
@u
1
0
1
−1
0
Fig. 7. Three-dimensional presentation of F–R sources for different applied strain rates.
1
103
A. Keyhani et al. / Computational Materials Science 104 (2015) 98–107
δε/δt=4.407e7
1.5
1
1
0.5
0.5
0
0
−0.5
−0.5
−1
−1
−1.5
−1.5
−1
−0.5
0
0.5
1
1.5
−1.5
−1.5
−1
δε/δt=7.242e7
1.5
1
0.5
0.5
0
0
−0.5
−0.5
−1
−1
−1
−0.5
0
0.5
1
1.5
−1.5
−1.5
−0.5
0
0.5
1
1.5
1
1.5
δε/δt=1.051e8
1.5
1
−1.5
−1.5
δε/δt=6.414e7
1.5
−1
−0.5
0
0.5
Fig. 8. Geometry of F–R sources in the plane of crack for different applied strain rates.
δε/δt=4.407e7
δε/δt=7.242e7
x 10
−3
δε/δt=6.414e7
−3
x 10
5
5
4
4
3
3
2
2
1
1
x 10
−3
δε/δt=1.051e8
−3
x 10
5
5
4
4
3
3
2
2
1
1
Fig. 9. Contours of the effective plastic strain for different applied strain rates.
shows the initial crack. The initial dislocation density is 0:33 lm2
for Sections 3.2–3.4. In Sections 3.2 and 3.3, the F–R sources are
placed symmetric with respect to the crack plane, whereas in
Section 3.4, the lower F–R sources are shifted about 0:1 lm outward of the crack direction in order to generate a mixed-mode
crack propagation mode due to microstructural non-symmetry.
104
A. Keyhani et al. / Computational Materials Science 104 (2015) 98–107
1
K dis
II ¼
Crack lentgh (μm)
0.95
0.9
0.85
δε/δt=4.407e7
δε/δt=6.414e7
δε/δt=7.242e7
δε/δt=1.051e8
No dislocation
0.8
0.75
0.7
5.2
5.4
5.6
5.8
6
Applied strain εzz
x 10
−3
Fig. 10. The crack length versus the applied strain for different applied strain rates.
3.1. Crack shielding due to the presence of dislocations
app
dis
dis
K I;II ¼ K app
I;II þ K I;II where K I;II and K I;II are the applied stress intensity
factor and the stress intensity factor caused by the presence of
dislocations, respectively. It is called dislocation shielding if the
stress intensity factor induced by presence of dislocations declines
the local stress intensity factor (i.e. K dis < 0). Conversely, it is
denoted dislocation anti-shielding if it intensiﬁes the local stress
intensity factor (i.e. K dis > 0). Analytical solutions for the induced
stress intensity factor due to the presence of dislocations are [62,63]
3Gb
h
pﬃﬃﬃﬃﬃﬃﬃﬃﬃ sin h cos
2
2ð1 tÞ 2pr
ð22:2Þ
where G, t and b are the shear modulus, the Poisson’s ratio and the
magnitude of the Burgers vector, respectively and the subscripts
refer to the crack mode. r and h are the crack tip polar coordinates,
as deﬁned in Fig. 4a.
The two-dimensional XFEM model of crack–dislocation interaction is depicted in Fig 4b. It is governed by the full crack–dislocation coupling of Eq. (11) based on the superposition strategy of
Fig. 2 and its corresponding boundary conditions. The stress ﬁeld
due to the presence of dislocations in an inﬁnite domain is
modeled by f B . While the modeling domain is considered large
enough to simulate an inﬁnite domain condition, f 1 is applied to
treat the boundary conditions due to crack free surfaces. Fig. 5
compares the predicted dislocation-induced stress intensity
factors with the reference analytical solution for h ¼ p=4 and
400b 6 r 6 4000b ’ 1 lm (Fig. 4a), which shows a good agreement. In addition, the mesh sensitivity of the results is evaluated
in Fig. 6 for three different element sizes, which clearly illustrate
that the results are mesh independent.
3.2. The effect of loading rate on mode-I cohesive crack propagation
The elastic ﬁeld arising from the presence of dislocations near the
crack tip may intensify or decline the stress intensity factor. The
local stress intensity factor is calculated by superimposing
K dis
I ¼
Gb
h
h
pﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 cos h cos sin h sin
2
2
2ð1 tÞ 2pr
ð22:1Þ
The isotropic mode-I crack propagation of the penny crack in an
inﬁnite domain for different loading rates is studied by the present
multi-scale framework. As previously mentioned, the dissipative
mechanism can only operate if the fracture is delayed sufﬁciently.
Therefore, by increasing the applied strain rate, the behavior
changes from ductile to brittle. In this simulation, different applied
strain rates in the [0 0 1] direction (perpendicular to the crack
plane) are applied. Geometry of the F–R sources for different
applied strain rates are shown in Fig. 7, which indicates that the
mobility of dislocations is decreased as the applied strain rate is
increased. In order to present the extent of crack propagation more
clearly, the results are shown in the [1 0 0][0 1 0] view in Fig. 8.
The corresponding effective plastic strains for different applied
strain rates are depicted in Fig. 9, which shows that the plastic region
due to dislocation motions becomes smaller by increasing the
Fig. 11. Three-dimensional presentation of F–R sources interacting with the random distribution of precipitates for different precipitate densities.
A. Keyhani et al. / Computational Materials Science 104 (2015) 98–107
105
Fig. 12. Geometry of F–R sources interacting with the random distribution of precipitates for different precipitate densities in [1 0 0][0 1 0] view.
1.1
1.05
Crack lentgh (μm)
1
0.95
0.9
No dislocation
No precipitate
1%
2%
3%
4%
0.85
0.8
0.75
0.7
5.2
5.4
5.6
5.8
Applied strain εzz
6
6.2
−3
x 10
Fig. 13. Variation of the crack length versus the applied strain for different
precipitate densities.
applied strain rate. Therefore, the contribution of dissipation mechanisms in the total fracture energy decreases. As a result, it is expected
that the crack propagates at smaller applied strain values when the
applied strain rate is larger, as clearly shown in Fig. 10.
resistance and diameter. The main purpose of this simulation is
to show anisotropy of the plastic domain near the crack tip in precipitate hardened materials. This anisotropy causes mixed-mode
crack propagation which is the subject of the next simulation. By
increasing the density of precipitates, the dislocation motion
becomes more restricted. The geometry of dislocations and precipitates for different dislocation densities are presented in
Figs. 11 and 12. To improve the presentation of the geometry of
dislocations, precipitates are not shown on the left half of the simulation domain.
Two limiting cases are possible for each applied strain rate: one
is when no dislocation exists and therefore the material behavior is
completely brittle, which means that the plastic dissipation does
not contribute to the total fracture energy. The other limiting case
is when no precipitates exist to resist dislocation motions.
Expectedly, any other cases with a random or arbitrary distribution
of precipitates may be assumed in between these cases. By increasing the density of precipitates, which results in higher resistance
against dislocation motions, the contribution of plastic dissipation
in the total fracture energy decreases signiﬁcantly, and the total
fracture energy is consequently reduced. Variation of the crack
length versus the applied strain for different precipitate densities
are presented in Fig. 13, which shows that the crack propagates
at lower applied strains by increasing precipitates density.
3.3. The effect of precipitate density on mode-I cohesive crack
propagation
3.4. Mixed-mode cohesive crack propagation due to microstructural
non-symmetry
The effect of random distribution of precipitates on the mode-I
crack propagation is studied for different precipitate densities. All
precipitates are considered to be impenetrable with the diameter
of D ¼ 100 nm; however, there is no limitation on the precipitate
Microstructural non-symmetry and anisotropy affect the crack
propagation direction. Such non-symmetry and anisotropy may
arise from different phenomena such as the random distribution
and various properties of F–R sources and precipitates. This ﬁnal
106
A. Keyhani et al. / Computational Materials Science 104 (2015) 98–107
Fig. 14. Three-dimensional propagated crack topologies for two considered cases: (a) mixed-mode crack propagation due to dislocation non-symmetry with respect to the
crack plane, (b) mixed-mode crack propagation due to plastic anisotropy from random distribution of precipitates.
σ (MPa)
σ (MPa)
zz
zz
1500
1500
1000
1000
500
500
0
(a)
0
(b)
Fig. 15. Stress plot ðrzz Þ for two cases: (a) dislocation non-symmetry with respect to the crack plane, (b) plastic anisotropy arising from random distribution of precipitates.
propagated crack topologies for both cases are plotted in Fig. 14.
It should be noted that only the normal axis to the crack plane
(and not dislocation paths) is scaled by a factor of 200.
The stress plots ðrzz Þ, shown in Fig. 15, illustrate that the stress
is not distributed symmetrically due to non-symmetric nature of
F–R sources with respect to the crack plane (Fig. 15a). This clearly
affects the crack propagation direction, as shown in Fig. 16.
Comparing the results in Figs. 15 and 16, shows that in case (a)
−4
15
x 10
Initial crack
Case (a)
Case (b)
Height (μm)
10
5
1.1
0
1.05
1
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
Crack length (μm)
Fig. 16. The initial and propagated crack paths; case (a) for dislocation nonsymmetry with respect to the crack plane and case (b) for random distribution of
precipitates in front of the crack.
Crack lentgh (μm)
−5
Mode−I: No dislocation
Mode−I: Precipitate density = 2%
Mixed−mode: Case (b)
Mixed−mode: Case (a)
Mode−I: Precipitate density = 0%
0.95
0.9
0.85
0.8
simulation is devoted to show the capability of the present framework to handle the mixed-mode crack propagation. Here, the
mixed-mode crack propagation due to microstructural non-symmetry and random distribution of precipitate are studied while
the crack is purely loaded in mode I. In case (a), the lower F–R
sources are shifted about 0:1 lm in the outward crack direction
to disturb the symmetry. In case (b), a random distribution of precipitates with the volume density of 2% is considered. The
0.75
0.7
5.2
5.4
5.6
5.8
6
Applied strain ε
zz
6.2
−3
x 10
Fig. 17. Comparison of the crack length versus the applied strain for mode-I and
mixed-mode crack propagations.
A. Keyhani et al. / Computational Materials Science 104 (2015) 98–107
when the crack tip meets the ﬁrst dislocation, it starts to propagate
upward due to high stress ratio resulted from the dislocation glide.
After a speciﬁc propagation in the direction of ﬁrst dislocation
glide, it is affected by the second dislocation and moves downward.
This procedure repeats until the crack tip passes all dislocations.
Afterwards, the crack propagates in an almost pure mode I, as
shown in Fig. 16.
In Fig. 15b, the stress ﬁeld is not symmetric due to plastic anisotropy arising from the random distribution of precipitates. The
ﬁrst lower F–R source faces high density of precipitates on its glide
plane, which consequently decreases the plastic strain. As a result,
the crack propagates along the ﬁrst higher F–R until it reaches the
second lower and higher F–R sources. After passing all F–R sources,
the crack propagates in mode I. Expectedly, the glide distance of
dislocations in Fig. 15b is less than Fig. 15a due to precipitate resistance against the dislocation glide.
In Fig. 17, the results of crack length versus the applied strain
for two cases of mixed-mode crack propagations are compared
with the previous results, in which the crack was forced to propagate in mode I. Even in mixed-mode crack propagation, the presence of precipitates causes the crack to propagate at lower
applied strain rates. However, after passing all F–R sources, all
cases converge to the same value (Fig. 15).
4. Conclusions
An extended three-dimensional multi-scale framework is presented to model plasticity and fracture in the vicinity of precipitates by merging the line dislocation dynamics analysis and the
extended ﬁnite element method. At the micro-scale, the modiﬁed
line dislocation dynamics is applied to analyze dislocation motions
and dislocation–precipitate interaction. Since the methodology
which is used to model precipitate is independent of the macro
analysis, it is not required to generate a mesh to deﬁne the precipitates geometry. At the macro-scale, XFEM is applied to provide
mixed-mode crack propagation without remeshing. Other available studies intentionally avoided mixed-mode crack propagation
by a symmetric assumption for the crack. The present multi-scale
framework is more efﬁcient since it has to no restrictions for
mixed-mode crack propagation and requires less computational
cost due to decreasing DOFs and avoiding the remeshing procedure. Various simulations are performed to study the crack tip
plasticity and demonstrate the capability of the present framework. It is shown that the crack propagates at lower loadings by
increasing the loading rates or precipitate density since these factors decrease the share of dissipation mechanism in the total fracture energy. In addition, it is shown that the crack may propagate
in mixed-mode even if it is loaded purely in mode I due to
microstructural non-symmetry. Considering elastic anisotropy
and improving some of the limitations of the present approach
can be considered for future studies.
Acknowledgements
The authors would like to acknowledge the technical support of
High Performance Computing Lab., School of Civil Engineering,
University of Tehran. Furthermore, the ﬁnancial support of Iran
National Science Foundation (INSF) is gratefully acknowledged.
References
[1]
[2]
[3]
[4]
[5]
J.R. Rice, J.S. Wang, Mater. Sci. Eng., A 107 (1989) 23–40.
A.N. Gulluoglu, C.S. Hartley, Modell. Simul. Mater. Sci. Eng. 1 (1992) 17.
A.N. Gulluoglu, C.S. Hartley, Modell. Simul. Mater. Sci. Eng. 1 (1993) 383.
L.P. Kubin, G. Canova, Scr. Metall. Mater. 27 (1992) 957–962.
T.A. Khraishi, H.M. Zbib, Comput. Mater. Sci. 24 (2002) 310–322.
107
[6] F. Roters, D. Raabe, G. Gottstein, Comput. Mater. Sci. 7 (1996) 56–62.
[7] R.N. Yellakara, Z. Wang, Comput. Mater. Sci. 75 (2013) 79–85.
[8] J.A. El-Awady, S. Bulent Biner, N.M. Ghoniem, J. Mech. Phys. Solids 56 (2008)
2019–2035.
[9] M.C. Fivel, C.F. Robertson, G.R. Canova, L. Boulanger, Acta Mater. 46 (1998)
6183–6194.
[10] A. Needleman, Acta Mater. 48 (2000) 105–124.
[11] E. Van der Giessen, A. Needleman, Modell. Simul. Mater. Sci. Eng. 3 (1995) 689.
[12] H. Yasin, H.M. Zbib, M.A. Khaleel, Mater. Sci. Eng., A 309–310 (2001) 294–299.
[13] M. Wallin, W.A. Curtin, M. Ristinmaa, A. Needleman, J. Mech. Phys. Solids 56
(2008) 3167–3180.
[14] H.M. Zbib, T. Diaz de la Rubia, Int. J. Plast 18 (2002) 1133–1163.
[15] T.K. Bhandakkar, A.C. Chng, W.A. Curtin, H. Gao, J. Mech. Phys. Solids 58 (2010)
530–541.
[16] H. Cleveringa, E. Van der Giessen, A. Needleman, J. Mech. Phys. Solids 48
(2000) 1133–1157.
[17] H. Cleveringa, E. Van der Giessen, A. Needleman, Mater. Sci. Eng., A 317 (2001)
37–43.
[18] V. Deshpande, A. Needleman, E. Van der Giessen, Acta Mater. 49 (2001) 3189–
3203.
[19] V. Deshpande, A. Needleman, E. Van der Giessen, Acta Mater. 50 (2002) 831–
846.
[20] V. Deshpande, A. Needleman, E. Van der Giessen, Acta Mater. 51 (2003) 1–15.
[21] S. Olarnrithinun, S.S. Chakravarthy, W.A. Curtin, J. Mech. Phys. Solids 61 (2013)
1391–1406.
[22] E. Van der Giessen, V. Deshpande, H. Cleveringa, A. Needleman, J. Mech. Phys.
Solids 49 (2001) 2133–2153.
[23] S.S. Shishvan, E. van der Giessen, Modell. Simul. Mater. Sci. Eng. 21 (2013)
065006.
[24] S.S. Shishvan, E. van der Giessen, Modell. Simul. Mater. Sci. Eng. 21 (2013)
065007.
[25] M.C. Fivel, C.R. Phys. 9 (2008) 427–436.
[26] H.M. Zbib, J. Eng. Mater. Technol. 131 (2009). 041209-041201.
[27] J.H. Rose, J. Ferrante, J.R. Smith, Phys. Rev. Lett. 47 (1981) 675–678.
[28] P. Bocca, A. Carpinteri, S. Valente, Int. J. Solids Struct. 27 (1991) 1139–1153.
[29] A. Carpinteri, G. Colombo, Comput. Struct. 31 (1989) 607–636.
[30] T. Belytschko, D. Organ, C. Gerlach, Comput. Methods Appl. Mech. Eng. 187
(2000) 385–399.
[31] E. Ooi, Z. Yang, Eng. Anal. Bound. Elem. 33 (2009) 915–929.
[32] T. Rabczuk, G. Zi, Comput. Mech. 39 (2007) 743–760.
[33] C.V. Verhoosel, M.A. Scott, R. de Borst, T.J. Hughes, Int. J. Numer. Meth. Eng. 87
(2011) 336–360.
[34] T. Belytschko, R. Gracie, G. Ventura, Modell. Simul. Mater. Sci. Eng. 17 (2009)
043001.
[35] T. Belytschko, T. Black, Int. J. Numer. Meth. Eng. 45 (1999) 601–620.
[36] J.M. Melenk, I. Babuška, Comput. Methods Appl. Mech. Eng. 139 (1996) 289–
314.
[37] N. Moës, J. Dolbow, T. Belytschko, Int. J. Numer. Meth. Eng. 46 (1999) 131–150.
[38] T. Strouboulis, I. Babuška, K. Copps, Comput. Methods Appl. Mech. Eng. 181
(2000) 43–69.
[39] N. Moës, T. Belytschko, Eng. Fract. Mech. 69 (2002) 813–833.
[40] G. Wells, L. Sluys, Int. J. Numer. Meth. Eng. 50 (2001) 2667–2682.
[41] G. Zi, T. Belytschko, Int. J. Numer. Meth. Eng. 57 (2003) 2221–2240.
[42] J.V. Cox, Int. J. Numer. Meth. Eng. 78 (2009) 48–83.
[43] P. Dumstorff, G. Meschke, Int. J. Numer. Anal. Meth. Geomech. 31 (2007) 239–
259.
[44] R. Gracie, J. Oswald, T. Belytschko, J. Mech. Phys. Solids 56 (2008) 200–214.
[45] J. Oswald, R. Gracie, R. Khare, T. Belytschko, Comput. Methods Appl. Mech. Eng.
198 (2009) 1872–1886.
[46] T. Belytschko, R. Gracie, Int. J. Plast. 23 (2007) 1721–1738.
[47] R. Gracie, T. Belytschko, Int. J. Numer. Meth. Eng. 86 (2011) 575–597.
[48] R. Gracie, T. Belytschko, Int. J. Numer. Meth. Eng. 78 (2009) 354–378.
[49] O. Skiba, R. Gracie, S. Potapenko, Modell. Simul. Mater. Sci. Eng. 21 (2013)
035003.
[50] S. Mohammadi, Extended Finite Element Method: For Fracture Analysis of
Structures, John Wiley & Sons, 2008.
[51] S. Mohammadi, XFEM Fracture Analysis of Composites, John Wiley & Sons,
2012.
[52] V.V. Bulatov, W. Cai, Computer Simulations of Dislocations, Oxford University
Press, 2006.
[53] W. Cai, J. Deng, K. Kang, A short course on DDLab and ParaDiS (2005).
[54] A. Keyhani, R. Roumina, S. Mohammadi, Submitted for publication, 2015.
[55] A.J. Ardell, Metall. Trans. A 16 (1985) 2131–2165.
[56] C.S. Shin, M.C. Fivel, M. Verdier, K.H. Oh, Phil. Mag. 83 (2003) 3691–3704.
[57] A. Takahashi, N.M. Ghoniem, J. Mech. Phys. Solids 56 (2008) 1534–1553.
[58] Y. Xiang, D.J. Srolovitz, L.T. Cheng, Acta Mater. 52 (2004) 1745–1760.
[59] K. Yashiro, F. Kurose, Y. Nakashima, K. Kubo, Y. Tomita, H.M. Zbib, Int. J. Plast.
22 (2006) 713–723.
[60] S.S. Shishvan, S. Mohammadi, M. Rahimian, E. Van der Giessen, Int. J. Solids
Struct. 48 (2011) 374–387.
[61] W. Cai, A. Arsenlis, C.R. Weinberger, V.V. Bulatov, J. Mech. Phys. Solids 54
(2006) 561–587.
[62] S.M. Ohr, J. Phys. Chem. Solids 48 (1987) 1007–1014.
[63] J.R. Rice, R. Thomson, Phil. Mag. 29 (1974) 73–97.