1-s2 0-S0013794412002809-main

Prévia do material em texto

Engineering Fracture Mechanics 95 (2012) 13–28
Contents lists available at SciVerse ScienceDirect
Engineering Fracture Mechanics
journal homepage: www.elsevier .com/locate /engfracmech
Simulation of crack propagation using a gradient-enriched ductile
damage model based on dilatational strain
T. Linse a,⇑, G. Hütter b, M. Kuna b
a Technische Universität Dresden, Institute for Solid Mechanics, 01062 Dresden, Germany
b Technische Universität Bergakademie Freiberg, Institute of Mechanics and Fluid Dynamics, Lampadiusstrasse 4, 09596 Freiberg, Germany
a r t i c l e i n f o a b s t r a c t
Keywords:
Damage mechanics
Ductile damage
Non-local damage
Strain-localization
Fracture mechanics
0013-7944/$ - see front matter � 2012 Elsevier Ltd
http://dx.doi.org/10.1016/j.engfracmech.2012.07.004
⇑ Corresponding author.
E-mail address: thomas.linse@tu-dresden.de (T.
An implicit gradient formulation based on the micro-dilatational approach is applied to the
GURSON–TVERGAARD–NEEDLEMAN ductile damage model. The stability of the model is proven
analytically for the considered class of materials. The damage model is implemented in
an implicit finite-element code in an UPDATED LAGRANGIAN formulation incorporating an objec-
tive time integration scheme and consistent tangent moduli for quadratic convergence. The
model is applied to simulate a fracture mechanical test of a typical pressure vessel steel.
The investigations are focused on the evolution of ductile damage and stress state at the
crack tip. The results show that the model captures the different stages of crack initiation
and propagation realistically. A condition for mesh-convergence in terms of the internal
length parameter is suggested.
� 2012 Elsevier Ltd. All rights reserved.
1. Introduction
When standard local continuum models are used to describe ductile damage of typical engineering metals [1,2], a high
sensitivity of the results with respect to the spatial discretization length is observed. Thus, the usual numerical implemen-
tation by means of the finite element method leads to a pathological mesh sensitivity for crack problems. For this class of
damage models, material softening due to the evolution of damage localizes in a zone of a size that is determined by the
element size [3,4]. The finite element size becomes an additional model parameter, which influences all other model param-
eters. As the mesh is refined, the localization zone narrows and finally vanishes. Hence, the dissipated energy converges to
zero for infinitesimal small elements, which contradicts the physical observation that both localization width and energy
dissipated by fracture are material properties. As a consequence, it is generally not possible to use the same damage param-
eters in dissimilar FE simulation with different element sizes. In addition, problems arise if the ductile–brittle transition re-
gion is to be modeled which is an important engineering task. In this region large plastic deformations and hence ductile
damage can precede final failure by cleavage. Therefore, it is necessary to account for both ductile and brittle damage. Dif-
ferent approaches are used to evaluate the probability of cleavage fracture at a certain load state such as the models by BERE-
MIN [5,6] or RITCHIE–KNOTT–RICE [7]. These models require very fine discretizations in order to resolve the high stress gradients
at crack tips accurately, which is in conflict with the pathological mesh sensitivity of local ductile damage models.
The effect of localization and the resulting mesh sensitivity of local strain-softening damage models have been addressed
using different physically motivated and phenomenological approaches, for critical comparisons see [8–10]. A common idea
of different approaches is to consider the influence of a finite region of surrounding material at a given material point. This
results in a formulation that somehow introduces a characteristic length scale into the model. So-called non-local
. All rights reserved.
Linse).
http://dx.doi.org/10.1016/j.engfracmech.2012.07.004
mailto:thomas.linse@tu-dresden.de
http://dx.doi.org/10.1016/j.engfracmech.2012.07.004
http://www.sciencedirect.com/science/journal/00137944
http://www.elsevier.com/locate/engfracmech
Nomenclature
A void nucleation pre-factor
a surface in current configuration
b relative dilatational plastic flow
c non-local length parameter
ep, eq volumetric and equivalent plastic strain
�ep non-local volumetric plastic strain
�e equivalent plastic strain of the matrix material
en void nucleation constant
E YOUNG’s modulus
f, f⁄ actual and effective void volume fraction
f0 initial void volume fraction
U yield function
fc void volume fraction for coalescence of voids
ff void volume fraction at final failure
fn void nucleation constant
fu effective void volume fraction at final failure
C region boundary
~g coupled non-local tangent stiffness
G shear modulus
hm hardening modulus of the matrix material
H denominator of plastic compliance
h hardening modulus
e EULER-ALMANSI strain tensor
DeE incremental logarithmic strain tensor
g vectorial weight function
F, Fr absolute and relative deformation gradient
f vector of nodal forces
G coupled non-local tangent stiffness tensor
I, I(4) second and forth order unit tensor
K element stiffness matrix
L velocity gradient tensor
N deviatoric flow direction
N shape function matrix
K bulk modulus
_kpl plastic multiplier
lc characteristic intrinsic length
l relative hydrostatic sensitivity
m POISSON’s ratio
X region
p hydrostatic pressure
q VON MISES equivalent stress
qi GTN parameters
rDII second principal value of stress deviator
rm yield stress of the matrix material
sn void nucleation constant
t time
v volume in the current configuration
w scalar weight function
A principal symbol
A geometrical stiffness
B shape function derivatives matrix
Cel elasticity tensor
del, dpl elastic and plastic rate of deformation tensor
d EULERian rate of deformation tensor
D; eD local and non-local tangent stiffness tensor
�ep nodal non-local strain vector
n normal vector in the current configuration
Rr relative rotation tensor
r CAUCHY stress tensor
14 T. Linse et al. / Engineering Fracture Mechanics 95 (2012) 13–28
_rJ JAUMANN stress rate
{r} elemental stress components vector
t̂ traction vector
ft̂g elemental traction vector
Ur relative right stretch tensor
u nodal displacement vector
v velocity vector
w spin tensor
T. Linse et al. / Engineering Fracture Mechanics 95 (2012) 13–28 15
modifications that can be applied to strain-softening damage models may be divided into three main types: formulations of
the integral type [11,12], explicit gradient formulations [13] and implicit gradient models [14,15]. While integral and impli-
cit gradient models are referred to as strongly non-local, explicit gradient formulations can be considered weakly non-local
[10]. Both explicit and implicit gradient methods can be correlated (under certain circumstances) to the integral expression
as shown in [14]. However, explicit gradient formulations result in some physical and numerical problems [9,16]. FOREST [17]
showed that implicit non-local formulations fit in a formal generalization of ERINGEN’s theory of micro-morphic media [18]
whereas ERINGEN’s micro-mechanical considerations leading to the theory can hardly be transferred to this formalism.
A number of non-local modifications have been proposed for ductile damage models, such as those by GURSON [19–22] and
ROUSSELIER [23], or classical isochoric [15,24] and crystal plasticity [25], respectively. There, spatial averaging is applied to
internal variables of the model, e.g. the void volume fraction or the equivalent plastic strain. These models were successfully
applied to simulate localization phenomena as well as crack initiation and propagation [e.g.[21,23,25–27]]. The focus was
mostly laid on the correct simulation of the global behavior of the investigated fracture mechanical specimen such as the
load displacement curves.
BARGELLINI ET AL. [28]proposed a higher-order damage model within the framework of a constrained micro-dilatational con-
tinuum. In this theory, the volumetric strain becomes an additional degree of freedom. A basic ansatz for the damage func-
tion was incorporated.
In the present study, a non-local modification of the micro-mechanically well-founded GURSON-model is employed based
on the micro-dilatational approach. It is shown analytically that the boundary-value problem arising from the application of
the model is well-posed. That means, convergence of finite element implementation with respect to the element size is en-
sured. The model is applied to ductile fracture in a compact tension specimen, whereby emphasis is put on the predicted
evolution of the fields near the crack tip.
2. Gradient-enriched ductile damage model
2.1. Ductile damage model
The continuum damage model established by GURSON and modified by TVERGAARD and NEEDLEMAN (GTN, [1,29,30]) provides
the basis for the constitutive equations applied in the present work. It was developed to describe the growth and coalescence
of initially present or nucleating voids in isotropic ductile materials.
The scalar porosity, i.e. the void volume fraction, f is used as a measure for ductile damage. The yield condition
U ¼ q
rm
� �2
þ 2q1f � cosh �
3
2
q2
p
rm
� �
� ð1þ ðq1f �Þ
2Þ 6 0; _kpl P 0; _kplU ¼ 0 ð1Þ
incorporates both hydrostatic stress p and the deviatoric part of the stress tensor rD according to the decomposition
p ¼ �1
3
r : I; q ¼
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
3
2
rD : rD
r
; N ¼ 3
2
rD
q
: ð2Þ
The associate flow rule for the plastic rate of deformation dpl can be written as
dpl ¼ _kpl @U
@r
¼ 1
3
_epI þ _eqN ð3Þ
with two scalar variables
_ep ¼ � _kpl
@U
@p
; _eq ¼ _kpl
@U
@q
ð4Þ
describing the dilatational and deviatoric plastic flow for convenient use in the next section.
16 T. Linse et al. / Engineering Fracture Mechanics 95 (2012) 13–28
It is assumed that elastic strains are small compared to plastic strains and the rate of deformation tensor can be decom-
posed additively into elastic and plastic part d = del + dpl. Using a hypoelastic–plastic formulation, the objective Jaumann
stress rate _rJ is determined by the rate of elastic deformation and the constant elasticity tensor:
_rJ ¼ Cel : del ¼ Cel : ðd� dplÞ ð51Þ
Cel ¼ 2GIð4Þ þ K � 2
3
G
� �
II: ð52Þ
In the above equation, I(4) and I are the fourth and second order unit tensors, K and G the bulk and shear modulus,
respectively.
The yield stress of the matrix material rm in (1) may be a function of the equivalent plastic strain e whose evolution equa-
tion reads
_�e ¼ �p
_ep þ q _eq
ð1� f Þrm
: ð6Þ
The porosity evolves due to growth of existing voids and the nucleation of additional voids at non-metallic inclusions
_f ¼ _f G þ _f N: ð7Þ
The existence of voids allows for volumetric plastic deformation since the matrix material is incompressible, thus void
growth is related to the dilatational part of plastic flow
_f G ¼ ð1� f Þ _ep; ð8Þ
whereas nucleation of voids is driven by the equivalent plastic matrix strain until the total nucleated void volume fraction
equals the volume fraction fn of nuclei initially in the matrix, as proposed by CHU and NEEDLEMAN [31]
_f N ¼ Að�eÞ _�e ð91Þ
Að�eÞ ¼ fn
sn
ffiffiffiffiffiffiffi
2p
p exp �1
2
�e� en
sn
� �2" #
: ð92Þ
Usually, void coalescence starts in the microstructure at a volume fraction of voids fc < 1. To incorporate this effect, an effec-
tive void volume fraction f⁄ was proposed by TVERGAARD and NEEDLEMAN [30]
f � ¼ f �ðf Þ ¼
f f 6 fc
fc þ ðf � fcÞK fc < f 6 ff
fu ff < f
8><>:
K ¼ fu�fcff�fc ; f u ¼
1
q1
:
ð10Þ
Thereby, coalescence of voids is modeled in a simplified form using an accelerated damage growth for f > fc. The total loss of
stress carrying capacity, i.e. material failure, occurs when ff is reached.
2.2. Non-local modification
As long as an internal length is incorporated and the local behavior is recovered for the elastic case, a number of variables
may be formulated in a non-local way for the considered class of softening materials. The mechanics of generalized continua
(for a recapitulation see e.g. [32]) offers one opportunity to take the characteristic length of a material into consideration.
Hereby, additional degrees of freedom are utilized that are independent from the translational degrees of freedom.
Motivated by the approach of generalized continua, the GTN continuum damage model is modified by replacing the dila-
tational part of the plastic strain rate _ep by its non-local spatial average _�ep in the evolution equation for the growth of existing
voids
_f nlG ¼ ð1� f Þ _�ep ð11Þ
_f ¼ _f nlG þ _f N: ð12Þ
All other equations of the GTN model remain untouched. For the considered case of elastic strains being small compared to
plastic strains, this approach is denoted micro-dilatational [32], i.e. the change of the micro-volume becomes an additional
degree of freedom. The local behavior is sustained for the elastic case.
Implicit gradient models are numerically advantageous compared to other non-local methods [14,15]. Therefore, the
additional field variable �ep is found by solving the partial differential equation of the HELMHOLTZ type
�ep � cr2�ep : I ¼ ep: ð13Þ
T. Linse et al. / Engineering Fracture Mechanics 95 (2012) 13–28 17
In the field Eq. (13), the local dilatational part of the plastic strain ep acts as a source term, the second order tensorr2�ep is the
second gradient of �ep with respect to the current configuration, i.e. r2():I is the LAPLACE operator. The parameter
ffiffiffi
c
p
with a
dimension of length can be regarded as an internal length that controls the influence of the surrounding material (c > 0).
To complete the implicit gradient formulation, the homogeneous NEUMANN boundary condition is chosen
r�ep � n ¼ 0: ð14Þ
Although rather mathematically motivated, this choice ensures that the overall plastic volume change of the continuum is
preserved
Z
v
�ep dv ¼
Z
v
ep dv ð15Þ
and �ep ¼ ep holds for homogeneous deformations. See [9] for a discussion of the effects of different boundary conditions on
the behavior of implicit gradient models.
In previous work [22,23] averaging was applied to the rates of f or f⁄. Note that choosing the non-local variable �ep based on
the micro-dilatational approach is numerically advantageous since ep can adopt arbitrary large values and is not bound to
certain limits, while f and thus f⁄ are. Furthermore, a model based on averaging f⁄ can show locking effects when f approaches
ff, as discussed for the one-dimensional case in [33].
The Eqs. (13) and (14) for the non-local field variable �ep are coupled with the standard elastoplastic boundary value prob-
lem governed by the static equilibrium in the absence of volume forces
r � r ¼ 0: ð16Þ
2.3. Stability analysis
2.3.1. General non-local material
Within a classical theory of plasticity with associated flow a material instability coincides with a loss of ellipticity of the
underlying system of partial differential equations so that bifurcation modes of infinitesimal width, i.e. jumps in the gradi-
ents of the essential field variables, become possible, see e.g. [4]. These modes of infinitely small wave length are responsible
for the observed mesh sensitivity. Thus, the regularizing effect of a general class of non-local strain-softening damage models
including the proposed gradient-enriched GTN model can be checked by analyzing the type of the underlying system of par-
tial differential equations. If the latter remains elliptic, a material instability leads to a bifurcation into modes of finite wave
length.
Following the approach of [34], the conditions for continuing equilibrium are obtained by differentiating (16) and the
HELMHOLTZ Eq. (13) with respect to time. For this purpose, the differential operators are transformed to the reference config-
uration by the deformation gradient before the differentiation with respect to time is carried out. Subsequently,changing
back to the current configuration yields
r � ð _rJ þ A : LÞ ¼ 0 ð171Þ
cr2 _�ep : I � 2cr2�ep : d ¼ _�ep � _ep: ð172Þ
Thereby, the geometrical stiffness A (see e.g. [34]) and correspondingly the second term in (172) enter through the differen-
tiation of the deformation gradient in the product rule. The term L =rv denotes the velocity gradient.
A general material law reads
_rJ ¼ D : dþ G _�ep ð181Þ
_ep ¼ eD : dþ ~g _�ep ð182Þ
with stiffness tensors D, G, eD and ~g. Combining Eqs. (17) and (182) leads to
r � ½ðDþ AÞ : Lþ G _�ep� ¼ 0 ð191Þ
cr2 _�ep : I þ ðeD � 2cr2�epÞ : dþ ð~g � 1Þ _�ep ¼ 0: ð192Þ
The principal symbol of this system of partial differential equations for the velocities v and the rate of non-local strains _�ep is
AðnÞ ¼
n � ðDþ AÞ � n 0
0 cn � n
� �
: ð20Þ
Ellipticity is lost whenever AðnÞ ¼ 0 for any n – 0, i.e. when
c n � n det½n � ðDþ AÞ � n� ¼ 0: ð21Þ
18 T. Linse et al. / Engineering Fracture Mechanics 95 (2012) 13–28
Since c > 0 the pre-factor of the determinant is positive and thus condition (21) restricts the tangent stiffness D but not the
coupling coefficients G, eD or ~g. Whether the boundary value problem is elliptic and therefore well posed is only determined
by that part of the constitutive description that is not affected by the non-local variable. The same conclusion was drawn in
[34] for an explicit gradient-enriched formulation. REUSCH [35] investigated the bifurcation behavior of an implicit gradient-
enriched damage model. As mentioned above, a loss of ellipticity coincides with the possibility of bifurcation modes of infin-
itesimal wave length. In this limit case, the bifurcation condition of this author corresponds to the condition (21) for a loss of
ellipticity.
2.3.2. Non-local GTN-model
In order to obtain a representation according to the general material law (182) for the proposed non-local GTN-model,
firstly the plastic rate of deformation is computed by inserting the flow rule (3) into the consistency condition _U ¼ 0 and
eliminating the plastic multiplier _kpl as
dpl ¼ @U
@r
@U
@r : _r
J
H
þ @U
@r
@U
@f ð1� f Þ _�ep
H
ð221Þ
with
H ¼ �
@U
@f Aþ @U@rm hm
� �
r : @U
@r
ð1� f Þrm
: ð222Þ
The term hm ¼ drm=d�e describes the hardening tangent modulus of the matrix material. Inserting dpl into the elastic law (52)
allows to identify the generalized tangent stiffness tensors. As shown, the stiffness with respect to the non-local strain rate _�ep
does not affect the possible loss of ellipticity and so only the first part of (221) needs to be investigated. The localization for a
general dilatational plastic material has been studied by RUDNICKI and RICE [36]. They assumed a constitutive behavior
dpl ¼ 1
h
1ffiffiffi
3
p N þ b
3
I
� �
1ffiffiffi
3
p N þ l
3
I
� �
: _rJ: ð23Þ
The symbol h denotes a hardening modulus and the values l and b measure the pressure sensitivity of the yield condition
and the dilatancy of plastic flow under pure shearing. For this material RUDNICKI and RICE derived that no localization occurs if
h
E
> � NIIffiffiffi
3
p þ l
3
� �2
ð24Þ
in the considered case of associative plastic flow (l = b) within a first order expansion of the rotation terms. Therein, E is
YOUNG’s modulus and NII is the second principal value of the deviatoric flow direction N.
As performed in [37] for the classical GURSON-model and for an integral non-local formulation in [20], inserting the yield
condition (1) and comparing the first part of (221) with the constitutive description (23) of RUDNICKI and RICE allows to identify
the parameters h and l. Using these values condition (24) can be rearranged to
rm
3Eð1� f Þ
hm
rm
� A df
�
df
q1ch� q21f �
1þ ðq1f �Þ
2 � q1f �ð2ch� RshÞ
" #
> � ð3N � q1q2f
�shÞ2
4½1þ ðq1f �Þ
2 � q1f �ð2ch� RshÞ�
; ð25Þ
where the following abbreviations have been introduced for convenience
R ¼ �3
2
q2
p
rm
ch ¼ cosh R sh ¼ sinh R N ¼ rDII=rm: ð26Þ
Thus, for a hardening matrix material (hm > 0) and without nucleation (A = 0) the material is always stable, since the term on
the left side is positive whereas the one on the right side is non-positive, as seen in [20]. The denominator on the right side of
(25) is always positive.
If nucleation is considered (A > 0), typically the parameters are such that voids nucleate in the initial and thus sub-critical
state of plastic deformation (f⁄ < fc). In this region void volume fraction is still small (f� 1) and even f � exp (R)� 1 holds
[37]. Under these circumstances (25) reduces to
rm
3E
hm
rm
� Aq1ch
� �
> �1
4
ð3N � q1q2f �shÞ
2
: ð27Þ
In this stage, the matrix material hardens considerably such that the hardening modulus typically lies in the range hm/
rm J 10. In the considered case, nucleation occurs for certain plastic deformation leading to a nucleation pre-factor
A [ 0.1. Even for pure (stably dilatational) crack tip blunting the relative hydrostatic stress in front of the crack tip lies below
R 6 4 [38] such that ch [ 25. Consequently, condition (27) is fulfilled and the material behavior for the non-local GTN-mod-
el is well-defined during all stages of ductile deformation for the material under consideration.
If a higher nucleating porosity is envisaged for the initial range so that the left side of (25) could become negative, the
easiest solution is to incorporate the nucleating porosity fn into the initial void volume fraction f0. If considerable void nucle-
T. Linse et al. / Engineering Fracture Mechanics 95 (2012) 13–28 19
ation and growth occur simultaneously, a non-local modification of the evolution Eq. (92) in an analogous manner as per-
formed for the void growth term could be applied. The further partial differential equation according to (13) for the equiv-
alent plastic strain formulated for this purpose offers the opportunity to introduce a second intrinsic length scale for the
nucleation.
2.4. Numerical implementation
The non-local modification of the GTN model is implemented into the commercial FEM software ABAQUS via a two-dimen-
sional isoparametric eight-node user element in an UPDATED LAGRANGIAN formulation [39] for plane strain.
With the help of the vectorial and scalar test functions g and w, the governing Eqs. (16) and (13) are written in weak form
to give the required expressions for implementation in the finite element method
Z
v
r � rg dv ¼
Z
a
t̂ � g da ð281Þ
c
Z
v
r�ep � rw dv ¼
Z
v
ðep � �epÞw dv; ð282Þ
where t̂ is the external traction vector. The finite element discretization is done using quadratic interpolation functions for
the displacement degrees of freedom u. The non-local volumetric part of the plastic strain �ep is interpolated using linear
functions, thus each of the corner nodes of the finite element is enriched by one additional degree of freedom contained
in the vector �ep. After consistent linearization [39], the system of algebraic equations for the increments of the degrees of
freedom (vectors Du and D�ep) resulting from (28) is rewritten into matrix form for each element as
Kuu Ku�ep
K�epu K�ep�ep
" #
Du
Dep
� �
¼ f
u
fep
" #
; ð29Þ
which is assembled using the entries
Kuu ¼
Z
X
BuT½Dre�Bu
� �
dXþ
Z
X
KGdX ð301Þ
Ku�ep ¼
Z
X
BuT½Dr�ep �N
�ep dX ð302Þ
K�epu ¼ �
Z
X
N�ep ½Depe�B
u dX ð303Þ
K�ep�ep ¼
Z
X
N�ep
Tð1� Dep�ep ÞN
�ep þ cB�ep TB�ep
� �
d X; ð304Þ
where ½Dre�; ½Dr�ep � and ½Depe� being matrices of dimension 4 � 4, 1 � 4 and 4 � 1, respectively. They contain the algorithmic
derivatives of stress with respect to strain (Dre) and with respect to the non-local volumetric plastic strain ðDr�ep Þ. The deriv-
atives of the local volumetric plastic strain with respect to strain and with respect to the non-local micro-dilatational plastic
strain are represented by Depe and Dep�ep . The corresponding force vectors are given as
fu ¼
Z
C
NuTft̂g dC�
Z
X
BuTfrg dX ð311Þ
feq ¼
Z
X
N�ep
Tðep � N�ep �epÞ þ cB�ep
T
B�ep�epdX: ð312Þ
In the above equations, Nu and N�ep are matrices containing interpolation functions for each degree of freedom, the corre-
sponding derivativesare found in Bu and B�ep . The terms ft̂g and {r} are external load vector and stress tensor in VOIGT nota-
tion. The geometric stiffness matrix KG in (301) accounts for the non-linear strain–displacement relations, for a derivation see
e.g. [40]. A similar term arising in (303) is neglected. Standard GAUSSian integration is applied.
2.4.1. Incrementally objective integration algorithm
Since the proposed damage model is formulated using a hypoelastic–plastic formulation, the numerical integration of the
constitutive equations in rate form must be incrementally objective [41]. This can be achieved by applying rotation neutral-
ization, allowing the use of standard algorithms for the time integration. Here, a rotationally neutralized algorithm proposed
by WEBER [42] is applied to integrate (51).
For the motion between t and t+ (t 6 t+) an orthogonal rotation tensor is defined by the initial value problem
_Rr ¼ w � Rr ð321Þ
20 T. Linse et al. / Engineering Fracture Mechanics 95 (2012) 13–28
RrðtÞ ¼ I; ð322Þ
where w is the spin tensor. The relative rotation tensor is obtained by polar decomposition of the relative deformation
gradient
FrðtþÞ ¼ FðtþÞ � F�1ðtÞ ¼ Rr � Ur: ð33Þ
Resulting from (32), the rate of the rotated stress
~r ¼ RTr � r � Rr ð34Þ
contains only material time derivatives
_~r ¼ RTr � _rJ � Rr: ð35Þ
As shown in [42], the time integration algorithm
~rðtþÞ ¼ RTr � rðtþÞ � Rr ¼ rðtÞ þ DrðDeEÞ ð36Þ
is incrementally objective if the deformation path between t and t+ is invariant under superposed rigid body motions, which
is achieved by defining a strain measure that is a function of the relative right stretch tensor only
DeE ¼ ln Ur: ð37Þ
The time integration (36) in the rotated (material) configuration is performed using a backward EULER method based on the
scheme for the integration of pressure-dependent plasticity described in [43]. In the following, rate forms are replaced by
increments. If stress and internal variables are given from the last converged state at t, a trial stress is calculated using
the relative strain increment. Furthermore, void growth caused by the non-local volumetric plastic strain increment is
considered
~rtrial ¼ rðtÞ þ Cel : DeE ð381Þ
f trial ¼ f ðtÞ þ ð1� f ðtÞÞD�ep: ð382Þ
If the new stress state causes plastic deformation according to (1), the following non-linear system of equations is solved
iteratively using a standard NEWTON scheme, thereby choosing Dep and Deq as primary unknowns
0 ¼ Uðp; q; �e; f Þ ð391Þ
0 ¼ Dep
@U
@q
þ Deq
@U
@p
ð392Þ
D�e ¼ �Deppþ Deqqð1� f Þrmð�eÞ
ð393Þ
Df ¼ Að�eÞD�e ð394Þ
p ¼ ptrial þ KDep ð395Þ
q ¼ qtrial � 3GDeq ð396Þ
f ¼ f trial þ Df ð397Þ
�e ¼ �eðtÞ þ D�e: ð398Þ
Note that for the iterative solution procedure a well-suited initial value Dep – 0 can be found in accordance with the limit of
R as discussed in section 2.3.2. After convergence, stress and local volumetric plastic strain are updated
~rðtþÞ ¼ �pI þ 2
3
qN ð401Þ
epðtþÞ ¼ epðtÞ þ Dep: ð402Þ
Finally, the updated stress needs to be rotated back into the current configuration
rðtþÞ ¼ Rr � ~rðtþÞ � RTr : ð41Þ
T. Linse et al. / Engineering Fracture Mechanics 95 (2012) 13–28 21
2.4.2. Consistent tangent moduli
Following the procedure proposed by [44], the algorithmic tangent moduli in the rotated configuration are found
explicitly without the need of expensive matrix inversions (for an outline of the derivation of Da; Ga; eDa; ~ga see Section A)
d~r ¼ Da : dDeE þ Ga dD�ep ð421Þ
dDep ¼ eDa : dDeE þ ~gadD�ep: ð422Þ
For a quadratic convergence of the global NEWTON method, the consistent tangent moduli needed in (30) can be obtained by
rotating the algorithmic moduli (42) back into the global configuration
½Dre� ¼ ½Rr � Rr � Da � RTr � R
T
r � ð431Þ
½Dr�ep � ¼ Rr � G
a � RTr
	 
ð432Þ
½Depe� ¼ Rr � eDa � RTrh i ð433Þ
Dep�ep ¼ ~ga: ð434Þ
The matrices ½Dre�; ½Dr�ep � and ½Depe� are obtained from the tensors in brackets on the right hand side using VOIGT notation.
3. Application and results
The implemented non-local modification of the GTN model is applied to simulate a fracture mechanical test. The simu-
lations are carried out to investigate the influence of the internal length
ffiffiffi
c
p
on the fracture toughness for small amounts of
crack growth.
The geometry of a compact tension specimen (CT, W = 50 mm, A0 = 0.57 W, H = 1.2 W) was chosen for the simulations
(Fig. 1). The CT specimen is modeled using meshes of different element sizes ranging from 0.1 mm � 0.1 mm for the coarsest
mesh to 0.00625 mm � 0.00625 mm for the finest one in the vicinity of the crack tip and along the ligament. Due to sym-
metry, only one half of the specimen is modeled, i.e. displacements perpendicular to the plane of symmetry are set to zero
along the whole ligament including the crack tip (Fig. 1). The plane strain formulation applied in the simulations is justifiable
for side-grooved specimens, the results are referred to an effective thickness BN = 20 mm.
The material parameters are chosen for a typical German pressure vessel steel (22NiMoCr3-7, A508 Cl.3), since in future
works its failure behavior is to be simulated in the brittle-ductile transition region. YOUNG’s modulus and POISSON’s ratio are
E = 210 GPa and m = 0.3. The yield curve of the material at room temperature was obtained by means of the small-punch-test
[45,46] and is given via the parameters r1 = 494 MPa, r2 = 200 MPa, r3 = 900 MPa, k = 18.4, eL = 0.01 of a modified VOCE
function.
Fig. 1. FE-Model of the side-grooved compact tension (CT) specimen.
Fig. 2. Load–displacement-curves obtained for the local model (dotted lines) and for the gradient-enriched model (solid lines,
ffiffiffi
c
p
¼ 0:06 mm). The symbols
mark the displacement corresponding to a crack growth Da � 0.2 mm.
22 T. Linse et al. / Engineering Fracture Mechanics 95 (2012) 13–28
rmð�eÞ ¼
r1 0 6 �e 6 eL
r1 þ r2ð�e� eLÞ
þðr3 � r2 � r1Þð1� e�kð�e�eLÞÞ eL 6 �e 6 1
8><>: ð44Þ
The parameters of the ductile damage model are en = 0.15, fn = 0.015, sn = 0.05, fc = 0.05, ff = 0.2, q1 = 1.5, q2 = 1.0,
f0 = 1.78 � 10�4 (initial void volume fraction). The internal length parameter
ffiffiffi
c
p
has been varied.
Using the gradient-enriched damage model, load–displacement-curves calculated for different meshes converge to a un-
ique curve. The load corresponding to a certain crack growth does not decrease if the mesh is refined as observed for the local
model (Fig. 2). The evolution of damage is distributed among a finite region around the crack tip that is not limited by the
element size but can be controlled using the internal length parameter (Fig. 3).
The distribution of f perpendicular to and along the ligament is plotted in detail in Fig. 4. To show the evolution of damage
with increasing load, the diagram contains curves for a crack growth Da = 0.05 mm and Da � 3:5
ffiffiffi
c
p
. The latter corresponds
to the first occurrence of material failure. It should be mentioned that the amount of crack growth Da is defined by a con-
tinuous length measure as
Da ¼
Z W�A0
0
f � � f0
fu � f0
dx ð45Þ
with reference of x to the undeformed configuration.
The results show that the maximum value of the void volume fraction f does not stay at the original crack tip but moves to
a position in front of it (Fig. 4, right). Contrary to the local model, failure due to void growth does not occur first directly at the
crack tip but at a distance of approximately 2
ffiffiffi
c
p
along the ligament.
The mechanism is as follows. A plastic zone forms for arbitrary small loadings. In the initial stage this zone is small com-
pared to the intrinsic length
ffiffiffi
c
p
so that the non-local plastic strain �ep remains zero. Corresponding to the evolution Eq. (11)
no void growth occurs then but there is stable plastic deformation. As large deformations are involved, the crack tip blunts.
The GTN yield condition (1) incorporates the MISES and the hydrostatic stress so that the plastic deformation has deviatoric as
well as dilatational components. For a bluntingcrack tip the maximum hydrostatic stresses and thus the dilatational plastic
flow occur in front of the crack tip, see e.g. [47,38]. So if the size of the highly strained region of the plastic zone becomes
comparable to the intrinsic length, this local dilatational plastic flow induces an increase of the non-local strain together
with the corresponding void growth in front of the initial crack tip.
The gradient-enriched model is able to reproduce complete material failure, as shown in Fig. 5: damage continues to grow
in a region around the first occurrence of failure in front of the crack tip, where f exceeds ff and stress reduces to zero. How-
ever, some numerical problems regarding the convergence of the overall NEWTON scheme arise for very fine meshes (see
Fig. 3. Distribution of damage in the vicinity of the crack tip for a crack growth Da � 0.2 mm (Top to bottom: be = 0.1 mm, be = 0.05 mm, be = 0.025 mm,
0.0125 mm, 0.00625 mm. Left: local model, right: gradient-enriched,
ffiffiffi
c
p
¼ 0:06 mm. Highlighted region: f > fc).
T. Linse et al. / Engineering Fracture Mechanics 95 (2012) 13–28 23
Fig. 4. Distribution of the void volume fraction (normalized using void volume fraction at final failure) perpendicular to (left) and along (right) the ligament
(element edge length be = 0.00625 mm). The distances are normalized using the characteristic length lc, which is the internal parameter
ffiffiffi
c
p
for the gradient-
enriched model and the element edge length be for the local model.
Fig. 5. Distribution of stress (top, normalized using the initial yield stress r1) and void volume fraction (bottom, normalized using void volume fraction at
final failure) along the ligament for two different meshes (gradient-enriched model,
ffiffiffi
c
p
¼ 0:06 mm). The distance is normalized using the internal
parameter ðlc ¼
ffiffiffi
c
p
Þ. Dashed lines: evolution for
ffiffiffi
c
p
< be with an element-by-element-failure starting at the first element in front of the crack tip. Solid
lines: evolution for
ffiffiffi
c
p
> be. Each line represents the damage distribution for an increase of Da of 0.05 mm.
24 T. Linse et al. / Engineering Fracture Mechanics 95 (2012) 13–28
0
0
50
100
150
200
250
300
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4
Crack growth Δa [mm]
J-
in
te
gr
al
J
[N
/m
m
]
be=0.1 mm
be=0.05 mm
be=0.025 mm
be=0.0125 mm
be=0.00625 mm
Fig. 6. R-curves calculated for different meshes (be – element edge length; solid lines:
ffiffiffi
c
p
¼ 0:1 mm, dashed lines:
ffiffiffi
c
p
¼ 0:05 mm, dotted lines: local
model). The results for a given
ffiffiffi
c
p
> 0 are mesh-independent if be 6
ffiffiffi
c
p
=4 for the gradient-enriched model.
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4
Crack growth Δa [mm]
0
50
100
150
200
250
300
350
400
J-
in
te
gr
al
J
[N
/m
m
]
Fig. 7. R-curves calculated for different length parameters
ffiffiffi
c
p
(element edge length be = 0.00625 mm). Solid lines: be 6
ffiffiffi
c
p
=4, dashed lines be >
ffiffiffi
c
p
=4, dotted
line: local model. Typical fracture toughness values obtained from different specimen geometries for a crack growth Da � 0.2 mm for this material [23] are
plotted for comparison.
T. Linse et al. / Engineering Fracture Mechanics 95 (2012) 13–28 25
26 T. Linse et al. / Engineering Fracture Mechanics 95 (2012) 13–28
Fig. 2: calculations for curves ending at uLL < 1.2 mm were stopped due to convergence problems). These numerical problems
are caused by two effects. First, elements adjoining the crack tip that did not fail so far are distorting excessively due to the
blunting mechanism described above. This problem gets worse for finer meshes. Numerical problems with a non-local dam-
age model for small element sizes were observed in [23], too. Second, the increment of the local volumetric plastic strain ep,
acting as a source term in (13), is currently set to zero for integration points where complete failure occurred. Hence, the
local term in (312) remains constant throughout subsequent increments of an analysis. Nevertheless, the non-local differen-
tial Eq. (13) respectively its weak formulation (282) is not fulfilled identically but still has to be enforced although the ele-
ments can undergo arbitrary shear deformations. This problem is inherent to all implicit gradient enriched formulations of
the HELMHOLTZ-type (13) independent of the choice of the variable to be treated non-locally. This point has partly been ad-
dressed in [48,49] but to the authors’ knowledge no generally accepted solution exists until now.
Note that using the non-local model it might be expected that evolution of damage is spatially distributed for an internal
length of approximately the element edge length be or even smaller, especially in the case that f stays below fc. However, ifffiffiffi
c
p
is smaller than be and f grows larger than fc, it is very likely that strain localizes in one element, since the maximum value
of damage appears in a distance of approximately 2
ffiffiffi
c
p
in front of the crack tip. For this case an element-by-element-failure
starting at the crack tip will be observed, as seen in [23] and demonstrated for
ffiffiffi
c
p
smaller than be in Fig. 5.
Fig. 6 shows crack growth resistance curves (R-curves) calculated for length parameters
ffiffiffi
c
p
¼ 0:05 mm and 0.1 mm. For
comparison, Fig. 6 contains curves calculated using the same material properties but the local damage model. The influence
of discretization for the local model is obvious: since damage localizes in a zone that is restricted to the element edge length,
the curves do not converge to a unique solution as the mesh is refined. Using the gradient-enriched model and the internal
length
ffiffiffi
c
p
¼ 0:1 mm, results for element edge lengths be = 0.025 mm, 0.0125 mm and 0.00625 mm converge to a unique
solution, in the case of
ffiffiffi
c
p
¼ 0:05 mm equal curves are obtained for be = 0.0125 mm and 0.00625 mm. Hence, a necessary
resolution of the mesh to achieve convergence is proposed:
be 6
ffiffiffi
c
p
4
; ð46Þ
i.e. considering
ffiffiffi
c
p
being a material constant, Eq. (46) defines the largest element edge length to be used.
Note that in contrast to the local approach the employed non-local model predicts R-curves with a distinct plateau
whose level correlates with the dissipated energy up to the first occurrence of complete material failure for a crack
growth Da � 3:5
ffiffiffi
c
p
, which can be seen as macroscopic point of crack initiation. This realistic behavior is attributed to
the transition from the initial stage of stable blunting to the evolution of a damage zone as described above.However,
the aforementioned numerical problems arise here. As Fig. 6 indicates, this plateau in the R-curves is followed by a tear-
ing stage.
The influence of the internal length parameter on the calculated R-curves is demonstrated in Fig. 7. A higher crack growth
resistance is observed as
ffiffiffi
c
p
increases. Although plotted for completeness, it must be emphasized that condition (46) does
not apply to curves for
ffiffiffi
c
p
6 0:02 mm (dashed lines). For comparison, typical fracture toughness values obtained from dif-
ferent specimen geometries for a crack growth Da � 0.2 mm for this material [23] are plotted.
4. Summary
In this paper, a gradient-enriched extension of the GURSON–TVERGAARD–NEEDLEMAN ductile damage model was presented
corresponding to the micro-dilational framework. In the model the evolution of damage is controlled by introducing an
internal length that represents a characteristic material parameter independent of the mesh size. The stability of the model
was shown analytically for material softening induced by void growth. The numerical implementation was realized as user-
defined subroutines in the FEM software ABAQUS by means of incrementally objective hypoelastic large deformation
algorithms.
The proposed damage model was successfully used to simulate a fracture mechanical CT-test. The ability to regularize theboundary value problem could be verified by converging solutions for mesh refinement. The results showed that in contrast
to local damage approaches the employed non-local model predicts an initial stage of stable crack tip blunting followed by a
distinct point of crack initiation which is in accordance with experimental experience. Investigating the global behavior of
the specimen as well as the fields directly at the crack tip, a condition for the mesh size was derived from the numerical
simulations. However, it was pointed out that the initial blunting limits the possible mesh resolution also from below if
no measures are taken. Thus, to get an adequate solution of the boundary-value problem, the range of the element size at
the crack tip is relatively narrow. Possibly, the problem of element distortion during initial blunting can be handled with
special (collapsed) crack tip elements which is work in progress.
Acknowledgements
The financial support of the Projects 1501298 and 1501343 by the German Federal Ministry of Economics and Technology
(BMWi) and the Project KU 929/14-1 by the Deutsche Forschungsgemeinschaft (DFG) is gratefully acknowledged. The
authors thank the reviewers for their valuable hints helping to improve the manuscript considerably.
T. Linse et al. / Engineering Fracture Mechanics 95 (2012) 13–28 27
Appendix A. Derivation of the consistent tangent moduli
Using the notation from [44] with H1 = f, H2 ¼ �e, h1 and h2 from the right hand sides of Eqs. (397), (398), the expressions
Da; Ga; eDa; ~ga are obtained by applying the procedure described in [44]. An equation similar to (5.26) of [44] is derived from
the linearization of Eqs. (391), (392) as
A11 dDep þ A12 dDeq þ E1 dD�ep ¼ ½B11I þ B12N� : dr ðA:11Þ
A21 dDep þ A22 dDeq þ E2 dD�ep ¼ ½B21I þ B22N� : dr; ðA:12Þ
where the additional constants
E1 ¼ wij
@hj
@D�ep
Dep
@2U
@q@Hi
þ Deq
@2U
@p@Hi
 !" #
E2 ¼
@U
@Hi
wij
@hj
@D�ep
ðA:2Þ
are introduced since the evolution of damage depends on D�ep (382). Eq. (A.12) is rearranged to give
dDep ¼ ðC11I þ C12NÞ : Z : dDeE þ F1 dD�ep ðA:31Þ
dDeq ¼ ðC21I þ C22NÞ : Z : dDeE þ F2 dD�ep: ðA:32Þ
Substituting (A.32) into (5.24) of [44] gives (422), where the expressions D
a; Ga; eDa; ~ga are formulated as
Da ¼ d0Ið4Þ þ d1II þ d2NN þ d3IN þ d4NI eDa ¼ ~d0I þ ~d1N
Ga ¼ g0I þ g1N ~ga ¼ F1
ðA:4Þ
with the new constants
~d0 ¼ 3KC11 g0 ¼ �KF1
~d1 ¼ 2GC12 g1 ¼ �2GF2
F1 ¼ 1D ½E1ðA22 þ 3GB22Þ � E2ðA12 þ 3GB12Þ� F2 ¼ 1D ½E2ðA11 þ 3KB11Þ � E1ðA21 þ 3KB21Þ�:
ðA:5Þ
All other terms appearing in Eqs. (A.1)–(A.5) (di,Aij,Bij,Cij,wij,D,Z) are identical to those given in [44].
References
[1] Gurson AL. Continuum theory of ductile rupture by void nucleation and growth: part I – yield criteria and flow rules for porous ductile media. J Engng
Mater Technol, Trans ASME 1977;99(1):2–15.
[2] Rousselier G. Ductile fracture models and their potential in local approach of fracture. Nucl Engng Des 1987;105(1):97–111.
[3] Pijaudier-Cabot G, Baant ZP, Tabbara M. Comparison of various models for strain-softening. Engng Comput 1988;5(2):141–50.
[4] de Borst R, Sluys L, Mühlhaus H-B, Pamin J. Fundamental issues in finite element analyses of localization of deformation. Engng Comput
1993;10(2):99–121.
[5] Beremin F. A local criterion for cleavage fracture of a nucear pressure vessel steel. Metall Trans A 1983;14A:2277–87.
[6] Mudry F. A local approach to cleavage fracture. Nucl Engng Des 1987;105:65–76.
[7] Ritchie RO, Knott JF, Rice JR. On the relationship between critical tensile stress and fracture toughness in mild steel. J Mech Phys Solids
1973;21(6):395–410.
[8] Jirásek M. Nonlocal models for damage and fracture: comparison of approaches. Int J Solids Struct 1998;35(31–32):4133–45.
[9] Peerlings R, Geers M, de Borst R, Brekelmans W. A critical comparison of nonlocal and gradient-enhanced softening continua. Int J Solids Struct
2001;38(44–45):7723–46.
[10] Bazant ZP, Jirásek M. Nonlocal integral formulations of plasticity and damage: survey of progress. J Engng Mech 2002:1119–49.
[11] Bazant ZP. Imbricate continuum and its variational derivation. J Engng Mech 1984;110(12):1693–712.
[12] Pijaudier-Cabot G, Bazant ZP. Nonlocal damage theory. J Engng Mech 1987;113(10):1512–33.
[13] de Borst R, Mühlhaus H-B. Gradient-dependent plasticity: formulation and algorithmic aspects. Int J Numer Methods Engng 1992;35:521–39.
[14] Peerlings R, de Borst R, Brekelmans W, de Vree J. Gradient enhanced damage for quasi-brittle materials. Int J Numer Methods Engng
1996;39(19):3391–403.
[15] Engelen RAB, Geers MGD, Baaijens FPT. Nonlocal implicit gradient-enhanced elasto-plasticity for the modelling of softening behaviour. Int J Plast
2003;19(4):403–33.
[16] Jirásek M, Rolshoven S. Localization properties of strain-softening gradient plasticity models. Part i: strain-gradient theories. Int J Solids Struct
2009;46:2225–38.
[17] Forest S. Micromorphic approach for gradient elasticity, viscoplasticity, and damage. J Engng Mech 2009;135(3):117–31.
[18] Eringen AC, Suhubi ES. Nonlinear theory of simple micro-elastic solids–i. Int J Engng Sci 1964;2(2):189–203.
[19] Tvergaard V, Needleman A. Effects of nonlocal damage in porous plastic solids. Int J Solids Structures 1995;32(8–9):1063–77.
[20] Leblond J, Perrin G, Devaux J. Bifurcation effects in ductile metals with nonlocal damage. J Appl Mech 1994;61:236.
[21] Jackiewicz J, Kuna M. Non-local regularization for fe simulation of damage in ductile materials. Comput Mater Sci 2003;28(3–4):684–95.
[22] Reusch F, Svendsen B, Klingbeil D. A non-local extension of Gurson-based ductile damage modeling. Comput Mater Sci 2003;26:219–29.
[23] Samal M, Seidenfuss M, Roos E, Dutta B, Kushwaha H. Finite element formulation of a new nonlocal damage model. Finite Elem Anal Des 2008;44(6–
7):358–71.
[24] Geers MGD. Finite strain logarithmic hyperelasto-plasticity with softening: a strongly non-local implicit gradient framework. Comput Methods Appl
Mech Engng 2004;193(30–32):3377–401.
[25] Aslan O, Quilici S, Forest S. Numerical modeling of fatigue crack growth in single crystals based on microdamage theory. Int J Damage Mech
2011;20(5):681–705.
[26] Needleman A, Tvergaard V. Dynamic crack growth in a nonlocal progressively cavitating solid. Eur J Mech – A/Solids 1998;17(3):421–38.
28 T. Linse et al. / Engineering Fracture Mechanics 95 (2012) 13–28
[27] Reusch F, Hortig C, Svendsen B. Nonlocal modeling and simulation of ductile damage and failure in metal matrix composites. J Engng Mater Technol
2008;130(2):021007–9.
[28] Bargellini R, Besson J, Lorentz E, Michel-Ponnelle S. A non-local finite element based on volumetric strain gradient: application to ductile fracture.
Comput Mater Sci 2009;45(3):762–7.
[29] Tvergaard V. Influence of voids on shear band instbilities under plane strain conditions. Int J Fract 1981;17:389–407.
[30] Tvergaard V, Needleman A. Analysis of the cup-cone fracture in a rod tensile bar. Acta Metall 1984;32:157–69.
[31] Chu C, Needleman A. Void nucleation effects in biaxially stretched sheets. J Engng Mater Technol 1980;102(3):249–56.
[32] Forest S, Sievert R. Nonlinear microstrain theories. Int J Solids Struct 2006;43(24):7224–45.
[33] Jirásek M, Rolshoven S. Comparison of integral-type nonlocal plasticity models for strain-softening materials. Int J Engng Sci 2003;41(13–
14):1553–602.
[34] Benallal A, Tvergaard V. Nonlocal continuum effects on bifurcation in the plane strain tension-compression test. J Mech Phys Solids
1995;43(5):741–70.
[35] Reusch F. Entwicklung und Anwendung eines nicht-lokalen Materialmodells zur Simulation duktiler Schädigung in metallischen Werkstoffen. Ph.D.
thesis, Lehrstuhl für Mechanik, Universität Dortmund; 2003.
[36] Rudnicki JW, Rice JR. Conditions for the localization of deformation in pressure-sensitive dilatant materials. J Mech Phys Solids 1975;23(6):371–94.
[37] Needleman A, Rice J. Limits to ductility set by plastic flow localization. In: Koistinen DP, Wang N-M, editors. Mech Sheet Met Form. New York: Plenum
Publishing Corporation;1978. p. 237–67.
[38] Yuan H. Numerical assessments of cracks in elastic–plastic materials. Berlin Heidelberg New York: Springer; 2002.
[39] Bathe K-J, Ramm E, Wilson EL. Finite element formulations for large deformation dynamic analysis. Int J Numer Methods Engng 1975;9(2):353–86.
[40] Zienkiewicz O, Taylor R. The finite element method – 2: solid mechanics. 5th ed. Butterworth-Heinemann; 2000.
[41] Hughest T, Winget J. Finite rotation effects in numerical integration of rate constitutive equations arising in large-deformation analysis. Int J Numer
Methods Engng 1980;15:1862–7.
[42] Weber G, Lush A, Zavaliangos A, Anand L. An objective time-integration procedure for istropic rate-independent and rate-dependent elastic-plastic
constitutive equations. Int J Plast 1990;6(6):701–44.
[43] Aravas N. On the numerical integration of a class of pressure-dependent plasticity models. Int J Numer Methods Engng 1987;24:1395–416.
[44] Zhang Z. Explicit consistent tangent moduli with a return mapping algorithm for pressure-dependent elastoplasticity models. Comput Methods Appl
Mech Engng 1995;121:29–44.
[45] Linse T, Kuna M, Schuhknecht J, Viehrig H-W. Usage of the small-punch-test for the characterisation of reactor vessel steels in the brittle-ductile
transition region. Engng Fract Mech 2008;75:3520–33.
[46] Linse T, Kuna M, Schuhknecht J, Viehrig H-W. Application of the small-punch-test to irradiated reactor vessel steels in the brittle-ductile transition
region. J ASTM Int 2008;5(4).
[47] Rice JR, Johnson MA. The role of large crack tip geometry changes in plane strain fracture. In: Kanninen MF, editor. Inelast Behav Solids. New
York: McGraw-Hill; 1970. p. 641–72.
[48] Mediavilla J, Peerlings RHJ, Geers MGD. Discrete crack modelling of ductile fracture driven by non-local softening plasticity. Int J Numer Methods
Engng 2006;66(4):661–88.
[49] Geers MGD, de Borst R, Brekelmans WAM, Peerlings RHJ. Strain-based transient-gradient damage model for failure analyses. Comput Methods Appl
Mech Engng 1998;160(1–2):133–53.
	Simulation of crack propagation using a gradient-enriched ductile damage model based on dilatational strain
	1 Introduction
	2 Gradient-enriched ductile damage model
	2.1 Ductile damage model
	2.2 Non-local modification
	2.3 Stability analysis
	2.3.1 General non-local material
	2.3.2 Non-local GTN-model
	2.4 Numerical implementation
	2.4.1 Incrementally objective integration algorithm
	2.4.2 Consistent tangent moduli
	3 Application and results
	4 Summary
	Acknowledgements
	Appendix A Derivation of the consistent tangent moduli
	References