Engineering Fracture Mechanics 73 (2006) 2210–2228

[Link]/locate/engfracmech

Simulation of dynamic ductile crack growth using

strain-rate and triaxiality-dependent cohesive elements
a,*
M. Anvari , I. Scheider b, C. Thaulow a

a
Department of Engineering Design and Materials, NTNU, Trondheim, Norway
b
GKSS Research Center, Geesthacht, Germany

Received 28 June 2005; received in revised form 15 March 2006; accepted 26 March 2006
Available online 12 June 2006

Abstract

Rate-sensitive and triaxiality-dependent cohesive elements are used to simulate crack growth under quasi-static and
dynamic loading conditions. The simulations are performed for a middle-cracked tension M(T) specimen made of an alu-
minum alloy (6XXX series). To consider the effect of stress triaxiality and strain rate on the cohesive properties, a single
plane strain element obeying the constitutive equations of a rate-dependent Gurson type model has been used. The single
element is loaded under various stress biaxiality ratios and strain rates and the obtained stress–displacement curves are
considered as traction separation law for the cohesive elements. These curves are used for analyzing the aluminum
M(T) specimen. The qualitative effects of constraint, strain rate, inertia and stress waves on the energy absorption of
the specimen and crack growth are discussed.
 2006 Elsevier Ltd. All rights reserved.

Keywords: Cohesive zone modeling; Dynamic loading; Rate dependency; Stress triaxiality; Rate-dependent Gurson type model

1. Introduction

When a structure is subject to dynamic loading, some additional phenomena compared to quasi-static cases
have to be considered. The high speed of deformation introduces stress waves into the structure. Inertia effects
become significant and can affect the energy absorption of the structure. Strain rate increases the flow stress
and influences the strain hardening of metals. And last, but not least, high local temperature increase due to
adiabatic heating reduces the flow strength and counteracts the strain-rate hardening. Whether it is a station-
ary crack subject to impact loading, called ‘‘impact fracture’’, or it is a fast motion of the crack tip, named
‘‘fast fracture’’ [1], the mentioned phenomena can affect the fracture toughness significantly and it is important
to have an efficient analysis tool for simulating damage and crack growth including these influences.
Analytical and macroscopic approaches used in fracture mechanics have some limitations with respect to
the amount of plasticity allowed at the crack tip, constraint and geometry dependency. As Siegmund and

*
Corresponding author.
E-mail address: [Link]@[Link] (M. Anvari).

0013-7944/$ - see front matter  2006 Elsevier Ltd. All rights reserved.
doi:10.1016/[Link].2006.03.016
 M. Anvari et al. / Engineering Fracture Mechanics 73 (2006) 2210–2228 2211

Brocks [2] point out, to the present, ‘‘local approaches’’ are the only really successful methods for prediction of
crack growth resistance. In a local approach, in principle, the parameters of the respective model depend only
on the material, and not on the geometry. In this kind of approach, one can simulate ductile fracture either by
employing a micromechanical model of damage, which represents the micromechanism of void initiation,
growth and coalescence, or by using a phenomenological model for material separation and coupling it to
the surrounding undamaged elastic–plastic material. In the former approach, a representative volume element
(RVE) or ‘‘unit cell’’ is considered to study the respective mechanism. The most widely known model of this
kind for ductile damage calculation is the Gurson–Tvergaard–Needleman or GTN model [3]. Zhang et al. [4]
proposed a model called ‘‘complete Gurson model’’ as a combination of GTN and Thomason’s coalescence
criteria [5]. This model has been used for unit cell calculations in the present article.
A phenomenological local approach used for the numerical simulation of crack propagation is known as
the cohesive zone model. The idea is based on the pioneering work by Dugdale [6] and Barenblatt [7] who
introduced the strip-yield model. Both authors divided the crack into two parts: One corresponds to the phys-
ical length of the crack, which is stress free, and the other one is the fracture process zone, where yielding and
degradation of the material occur and which is loaded by a finite stress named cohesive stress. Later develop-
ments of cohesive models, particularly in combination with finite element method (e.g. [8–11]), considered the
cohesive stress as a function of material separation and not depending on the distance from the crack tip as
Barenblatt did. In a finite element representation of cohesive zone models, cohesive elements are introduced as
interface between continuum elements and damage occurs only in the interface elements which obey a consti-
tutive equation named traction separation law (TSL) (Fig. 1). Separation in these elements is calculated from
the difference of the displacements of the continuum elements adjacent to them. The maximum opening at
which the cohesive element completely fails is called critical separation, d0, and is one of the fracture para-
meters. The other fracture parameter is the maximum traction or cohesive strength, S. The area under the
TSL is the energy absorbed by the cohesive element and is known as the cohesive energy, C0:
Z d0
C0 ¼ T ðdÞ dd ð1Þ
0

If the shape of a TSL is known or presumed, having two of the aforementioned parameters is enough to define
the cohesive law.
By using cohesive zone modeling in a finite element analysis, it is possible to split the total dissipated energy
into energy dissipated by plastic deformation in the process zone and the energy of separation. In this way, it is
possible to evaluate the damage and deformation processes separately, but coupled. Besides, by entrusting the
nucleation, propagation, branching and other aspects of the fracture behavior of materials to a master cohe-
sive law, the amount of phenomenology is considerably reduced compared to theories of distributed damage.
Since a characteristic length is included in the TSL (cohesive energy is the work of separation per unit area),
the models with cohesive elements are not mesh dependent, which is always an issue for the GTN models [12].
Cohesive elements used in simulating ductile fracture are supposed to represent the mechanism of nucle-
ation, growth and coalescence of microscopic voids that initiate at the inclusions and second phase particles.

Fig. 1. Cohesive zone obeying a TSL and the surrounding undamaged elastic–plastic material.
2212 M. Anvari et al. / Engineering Fracture Mechanics 73 (2006) 2210–2228

The idea of the present contribution is to obtain the cohesive properties by studying the mechanical response of
a single element obeying a rate-dependent Gurson type constitutive equation. This sort of mechanism-based
cohesive model has also been proposed by other authors, e.g. [12–16]. It is well known that the cohesive strength
increases with strain rate [17–19] and this implies an increase in stress triaxiality [20,21]. Using a Gurson type
model provides a mean for exploring the effect of stress triaxiality on TSL parameters [12]. Rate sensitivity of
the cohesive zone is determined by applying different values of loading speed on the rate-sensitive single element
and calculating the cohesive parameters from the mechanical response. A similar procedure is performed by
applying different stress ratios on the single element to study the effect of stress trixiality on the TSL parameters.
The shape of the TSL, cohesive strength and critical displacement are fitted to meet the results of the single ele-
ment calculations for various strain rates and triaxialities. The parameters are then used for crack growth sim-
ulations of an aluminum middle-cracked tension M(T). In the model, the cohesive elements are both triaxiality
and rate dependent. Triaxiality-dependent cohesive elements have already been introduced to cohesive elements
by Siegmund and Brocks [2], whereas the effect of dynamics on cohesive parameters is usually written by a sep-
aration rate dependency (see e.g. [18]). Even though the implementation of a strain-rate effect is a bit more com-
plicated than a separation rate, it has the advantages that are described in Section 2.3.
Since the measures of strain rate and stress triaxiality are not available in the cohesive elements, they are
calculated and transferred from the solid elements adjacent to them, which is performed in a similar way as
proposed in [2] for only stress triaxiality values. The technique has also been used for strain-dependent cohe-
sive elements by Tvergaard and Hutchinson [22]. The analyses are performed in quasi-static and dynamic
cases. The effects of rate dependency, inertia, constraints and waves on mechanical response of the structure
are shown by the load–displacement diagrams and the results are discussed. The influence of the pointed out
phenomena on plastic energy dissipation and crack growth are illustrated, too. The influence of adiabatic heat-
ing has not been considered in the analyses.
The static and transient dynamic analyses are performed in ABAQUS/Standard non-linear finite element
code [23]. The rate-dependent complete Gurson model developed and implemented into ABAQUS as a user
defined material (UMAT) subroutine [24] is used for the single element calculations. The cohesive element cal-
culations are performed by a user defined element (UEL) subroutine developed by Scheider [25] and expanded by
the present authors for rate and triaxiality-dependent cohesive elements to be used in transient dynamic analysis.

2. Calculations

2.1. Material

Smooth round bars made of aluminum alloy 6XXX series were tested by split-Hopkinson tension bar in
SIMLab of NTNU [26]. The experimental stress–strain curves at different strain rates are approximated by
the following equation [27] as shown in Fig. 2:
   m
 N _
 ¼ r0 1 þ
r ð2Þ
0 _ 0
where r , , r0 and 0 stand for true stress, true plastic strain, yield stress, and a reference strain respectively. _
and _ 0 are the strain rate and a reference strain rate, respectively. N is known as strain hardening and m as
strain-rate hardening exponent. In this formulation, the strength of the material increases with the increase
of strain rate. The constitutive bahavior is based on von Mises plasticity with pure isotropic hardening. Vis-
cous, i.e. time dependent effects are not considered. The values considered for the material tested are:
r0 = 217 MPa, 0 = 0.002, _ 0 ¼ 150 s1, N = 0.0526, m = 0.05. Mass density of aluminum, q = 2700 kg/m3
has been considered in all of the dynamic calculations.

2.2. Single element calculations

For ductile metals, the effective behavior of a unit cell containing a void can be captured by using the
Gurson model [28] developed by Tvergaard [29] and Tvergaard and Needleman [3]. The latter, known as
 M. Anvari et al. / Engineering Fracture Mechanics 73 (2006) 2210–2228 2213

300

1000 1/s
290 800 1/s
500 1/s
280
QS
True Equivalent Stress (MPa)

270

260

250

240

230

220

210

200
0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09
True Plastic Strain

Fig. 2. The experimental stress–strain curves for different strain rates approximated by Eq. (6).

GTN uses a fixed critical void volume fraction, i.e. the point where coalescence starts, as a material parameter.
Koplik and Needleman [30] showed that this assumption works well only for low stress triaxiality cases, and
slightly worse for high stress triaxiality cases, which is the case in front of a crack tip. Zhang et al. [4] intro-
duced the Thomason’s coalescence criteria [5] into the Gurson model so that it can predict void coalescence. In
their approach which they called ‘‘complete Gurson model’’, they also proposed final void volume fraction to
be a function of initial void volume fraction. In the strain-rate-sensitive version of the model, the flow stress is
a function of local strain rate with the same strain-rate hardening exponent as that of the bulk material.
Details on the rate-dependent Gurson model is found elsewhere (e.g. [31,32]).
The question that arises is which of the phenomena influencing ductility because of dynamic loading
conditions is more important to be considered in modeling a single Gurson type element. Curran et al. [33]
considered rate-sensitive material for a hollow sphere of internal radius, a, subjected to internal constant pres-
sure, P0. The material was considered to be rigid, perfectly plastic with a viscosity coefficient. Their numerical
investigation for aluminum showed that the inertia effect becomes important for quite large voids (larger than
10 lm) and for a wide range of radii a, the maximum rate is governed by the viscous regime only. In other
words, as far as void size is smaller than a certain value, inertia can be ignored in modeling void growth under
dynamic loading. Considering the work done by Johnson [34], it is concluded
pffiffiffiffiffiffiffiffiffiffi that the inertia effect is impor-
tant in void growth only if the speed of crack growth is higher than 4 r0 =q where r0 and q are the flow stress
and the density of the material, respectively. This value is around 1100 m/s for a typical aluminum alloy. Liu
[35] also showed that increasing strain rate will increase the inertia effect, but this effect is significant only when
the strain rate is over 10,000 s1.
Based on these discussions, only rate sensitivity of the material has been considered. The effect of the adi-
abatic heating is ignored in all of the calculations.
In the present study, a single four node plane strain element with rate-dependent complete Gurson model as
constitutive equation is used. The initial void volume fraction of f0 = 0.002 is considered in the calculations as
suggested in [36], since no experimental data on micromechanical information is available for this parameter.
2214 M. Anvari et al. / Engineering Fracture Mechanics 73 (2006) 2210–2228

The analyses are performed for an element with the initial size of h · h, h = 1 mm, for different stress biaxial-
ities b, as shown in Fig. 3:
r11
b¼ ð3Þ
r22
If the material is incompressible, the stress triaxiality which is the ratio between mean normal stress and von-
Mises equivalent stress is given by [12]
1þb
H  pffiffiffi ð4Þ
3ð1  bÞ
The constant stress ratios in static loadings were obtained using modified Riks method implemented in
ABAQUS. The reason for using Riks method is the instability in load controlled calculations due to softening.
Fig. 3 shows the normalized traction separation obtained for different values of triaxiality. Fig. 4 shows the
variation of normalized strength and energy obtained from Fig. 3 versus triaxiality. The curves are fitted to
the calculated values.
In order to investigate the effect of strain rate on the stress–elongation behavior, quasi-static analyses were
performed with a rate-dependent model and high speed loading. Since strain rate was considered, Riks method
was not applicable any more. Thus, two linear springs were added to the model and the multi-point constraint
(MPC) subroutine in ABAQUS was used to obtain constant stress ratios. For low speed of applied loads, the
responses were the same as the Riks method. Fig. 5 shows the normalized traction separation behavior for
different triaxialities for a loading rate of 500 mm/s or, in other words, an initial strain rate of 500 s1 for
an element of initial length h = 1 mm. The figure shows that the effect of triaxiality is similar to the case with
rate-independent material behavior, but the values of cohesive strength and consequently cohesive energy are
different. To check this, various loading speeds were applied for the same triaxiality of 1.5 (b = 0.44). Fig. 6
shows that the curves are parallel and only the value of strength changes for the same maximum displacement.
Fig. 7 shows the curves fitted to the calculated values of strength and energy for the rate-dependent case at
constant loading rate (V = 500 mm/s) and different triaxiality values. Comparing this figure with Fig. 4, it
is observed that the fitted curves are very similar to the rate-insensitive case.

Fig. 3. The effect of triaxiality on the traction separation behavior of a single element obeying complete Gurson model.
 M. Anvari et al. / Engineering Fracture Mechanics 73 (2006) 2210–2228 2215

4.5 1.8

4 1.6
Strength

3.5 1.4

3 y=1.08 ln(x)+2.09 1.2

2.5 1

Γ /σ h
S/σ0

0
0
2 Energy 0.8

1.5 0.6

–1.3532
y=1.3654 x
1 0.4

0.5 0.2

0 0
0 1 2 3 4 5 6
Triaxiality

Fig. 4. Variation of normalized cohesive strength and cohesive energy versus applied triaxiality for rate-insensitive material.

4
H=5
H=4
H=3
3.5 H=2
H=1.5
H=1
3

2.5
0
σ /σ

2
22

1.5

0.5

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

U2/h

Fig. 5. Normalized traction separation behavior in different triaxialities for a load speed of 500 mm/s (initial strain rate of 500 s1).
2216 M. Anvari et al. / Engineering Fracture Mechanics 73 (2006) 2210–2228

3
V=1500 mm/s
V=1000
V=500
V=200
2.5 V=100
Rate independent

2
0
σ /σ

1.5
22

0.5

0
0 0.1 0.2 0.3 0.4 0.5 0.6

U2/h

Fig. 6. Normalized traction separation behavior at different load speeds (initial strain rate) and constant triaxiality (H = 1.5).

4.5 1.8

Strength
4 1.6

3.5 1.4
y=1.1571 ln(x)+2.1962
3 1.2

2.5 1
Γ /σ h
0

0
S/σ

Energy
2 0.8

1.5 0.6

y=1.497 x–1.3666
1 0.4

0.5 0.2

0 0
0 1 2 3 4 5 6
Triaxiality

Fig. 7. Variation of normalized cohesive strength and normalized cohesive energy versus different triaxialities in constant load speed of
500 mm/s.
 M. Anvari et al. / Engineering Fracture Mechanics 73 (2006) 2210–2228 2217

Fig. 8 shows how the values of cohesive strength and cohesive energy change with loading rate for a con-
stant triaxiality of 1.5. As it is shown, the value of the maximum traction changes as a function of the rate-
independent case for the same triaxiality multiplied by the relative strain rate to the power of m (strain-rate
hardening exponent).

2.3. Traction separation law

Since the cohesive model is a phenomenological model, authors use various formulations for defining the
shape of TSL and the cohesive values (e.g. [8,10,37,13,38,39]). Among these definitions, the one introduced
by Scheider and Brocks [39] has been used to reproduce the calculated damage and failure of the unit cell by
a cohesive element. This TSL has two shape parameters, d1 and d2, in addition to the cohesive parameters
already introduced and consists of three parts: (a) increasing traction or cohesion, (b) constant traction and
(c) decreasing traction or decohesion (see Fig. 9). The traction as a function of the separation in this model is
8    2
>
> 2 d  dd1 ; 0 < d < d1
>
< d1
T ¼ S 1; d1 < d < d2 ð5Þ
>
>  3  2
>
: 2 dd2  3 dd2 þ 1; d < d < d
d0 d2 d0 d2 2 0

By changing d1 and d2, one can have a variety of TSL shapes. In this formulation, the maximum separation is
2C0 1
d0 ¼ ð6Þ
S 1  23 dd1 þ dd2
0 0

Although it has been claimed that the shape of the TSL hardly influences the crack growth behavior
[13,37,40,41], there are a few investigations that show higher effects of the shape [19,42,43]. For example,

3 10

Energy 0.04
y=8.25 (x/150)
2.9 9

2.8 8

y=2.56 (x/150)0.04
Γ /σ h
S/σ0

2.7 7
0

2.6 6

Strength

2.5 5

2.4 4
0 200 400 600 800 1000 1200 1400 1600
Loading Rate (mm/s)

Fig. 8. Cohesive strength and cohesive energy versus load speed for a constant triaxiality of 1.5.
2218 M. Anvari et al. / Engineering Fracture Mechanics 73 (2006) 2210–2228

ΤN

Γ0

δ1 δ2 δ0 δN

Fig. 9. TSL shape proposed by Scheider and Brocks [39].

Scheider and Brocks [42] showed numerically that not only the shape of the TSL can affect the load–displace-
ment behavior, but also the TSL shapes that make the same results in a fracture specimen, might make dif-
ferent results in another specimen. Falk et al. [43] showed that if the TSL is initially rigid, the model
cannot capture crack branching by using a regular finite element mesh and it is possible when the TSL has
initial elasticity.
Regarding uncertainties like these, rate-dependent Gurson type model has been used as the basis for obtain-
ing the TSL.
It is possible to obtain the following expressions for the cohesive strength and energy using the mathemat-
ical approximations shown in Figs. 4, 7 and 8 as
 m
C0 _
 1:43H 1:36 r0 ð7Þ
h _ 0
 m
_
S  ð1:1 lnðH Þ þ 2:1Þr0 ð8Þ
_ 0
The variable h is the characteristic length scale of the Gurson model, which is usually based on the microme-
chanical structure of the material. A value of h = 0.1 mm is chosen in accordance to [44] who also investigated
an aluminum alloy of 6XXX series. H is the stress triaxiality. The expression in Eq. (8) is similar to the results
obtained analytically by Zhang et al. [19] which show that the curves of normalized traction along the cohesive
zone are parallel at various crack speeds and the cohesive strength increases with crack growth rates. It has
also been pointed out by Freund [45] that as the crack moves more rapidly, the material is deformed more
rapidly and a larger cohesive stress is required in order to achieve the requisite crack tip opening displacement.
Considering e.g. Fig. 5, it is recognizable that the piecewise TSL defined by Eq. (5) and shown in Fig. 5 fits
properly to the curves if the shape parameters d1 and d2 are adjusted accordingly. Combining Eqs. (6)–(8), the
critical separation reads
d0 2:8H 1:36 1
 ð9Þ
h 1:1 lnðH Þ þ 2:1 1  23 dd1 þ dd2
0 0

The comparison between the stress elongation behavior obtained from rate-dependent Gurson model calcu-
lations and the TSL approximation of Eq. (5) with the cohesive parameters from Eqs. (8) and (9) and
d1
d0
¼ 0:07 and dd20 ¼ 0:35 are shown in Fig. 10 for two different triaxiality values. Since the curves obtained from
the single element calculations and those based on Eq. (9) are not exactly the same, the critical displacements
are a slightly different for the same cohesive strength and energy, but the differences are reasonable.
It is worth noting that the initial stiffness of the cohesive elements is relatively high. Therefore, the separa-
tions are very small in the beginning until the cohesive strength is reached. This is an important indication for
the use of a strain rate instead of a separation rate-dependent cohesive model. The cohesive strength is
strongly influenced by the rate, and in the beginning most of the deformation takes place in the adjacent
 M. Anvari et al. / Engineering Fracture Mechanics 73 (2006) 2210–2228 2219

800

TSL (H=3)
Gurson (H=3)
700 TSL (H=1.5)
Gurson (H=1.5)

600

500
Traction (MPa)

400

300

200

100

0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
δ (mm)

Fig. 10. A comparison between traction separation behavior obtained from Gurson type model and the TSL proposed in Eq. (5).

continuum elements whereas the cohesive elements are still closed. On the other hand, the effect of the rate on
the softening branch of the traction separation law is of minor importance since the characteristic of the dam-
age behavior is then more affected by the cohesive energy parameter which is not much dependent on strain
rate (see Figs. 6 and 8).

3. Crack growth simulations and results

An M(T) specimen of unit thickness and dimensions of 100 · 100 mm2 has been modeled using four node
plane strain elements. The ratio of the initial crack length to the specimen’s width is a/W = 0.5. One row of
cohesive elements with initial zero height has been used at the ligament. Because of symmetry, one fourth of
the specimen has been modeled. Fig. 11 shows the specimen and the detailed mesh at the crack tip and a part
of the ligament. The smallest continuum element size belongs to the elements adjacent to the ligament and it is
0.1 · 0.1 mm2.
The surrounding continuum elements behave according to Eq. (2) and the cohesive elements based on Eqs.
(5), (7), (8) and (9). Each cohesive element possesses four nodes and two integration points. Because of the
symmetry, only the two upper nodes of cohesive elements are shared with the adjacent continuum element
and the other two obey the displacement conditions of symmetry. For each cohesive element, the actual values
of nodal displacement, triaxiality and strain rate are used to compute the current value of traction and element
stiffness. The values of triaxiality and strain rate are calculated in the continuum elements along the ligament
at all of the integration points using UVARM subroutine in ABAQUS. The averages of these values are then
calculated after the load increment and provided to the respective adjacent cohesive elements in the next incre-
ment as shown in Fig. 12. The error of this kind of explicit scheme is accepted since the time increments are
always chosen to be small enough.
Material separation process involving void growth and plastic deformation at microlevel are irreversible by
nature. This property has been implemented in the cohesive element formulation in a sense that a separation
2220 M. Anvari et al. / Engineering Fracture Mechanics 73 (2006) 2210–2228

Fig. 11. The M(T) specimen and a detailed finite element model of the crack tip.

Fig. 12. The triaxiality and strain-rate values are transferred to cohesive elements from the adjacent continuum elements.

already happened to an element remains after local unloading and relaxing the stress [25]. Local unloading can
happen because of global unloading, stress waves or crack branching.
The simulations performed are described in Table 1. The aim of the analyses is to check the influence of
constraint, rate sensitivity, inertia and elastic waves on the energy absorption of a typical M(T) specimen.
In the static simulations, the effect of triaxiality on the energy absorption is studied. In the dynamic analyses,
denoted by DYN-xxx, the effects of inertia and elastic waves are investigated as well as the influence of stress
triaxiality and rate sensitivity.
Fig. 4 shows that the change of cohesive energy (critical displacement) is very high for triaxiality values of
H < 1.5. The reason is that for small values of triaxiality, the initial void grows very little and therefore the
energy absorbed is mostly related to the plastic deformation and not damage in the cell model. For more
detailed discussions on the limitations of the application of cell model and the competition between cohesive

Table 1
Various simulations of M(T) specimen
Case Analysis Plasticity Cohesive zone
QS Quasi-static Rate independent Constant
QS-CTD Quasi-static Rate independent Triaxiality dependent
DYN Transient dynamic Rate independent Constant
DYN-MRD Transient dynamic Rate dependent Constant
DYN-MCRD Transient dynamic Rate dependent Rate dependent
DYN-MCRD-CTD Transient dynamic Rate dependent Rate and triaxiality dependent
 M. Anvari et al. / Engineering Fracture Mechanics 73 (2006) 2210–2228 2221

process zone and plastic flow (see [20]). It is also known that the M(T) specimen is a low constraint fracture
specimen. The analyses show that the stress triaxiality at the crack tip increases while the crack is growing.
Regarding these discussions, H = 1.2 was considered as an initial stress triaxiality value used in the simula-
tions. Below this value, the cohesive parameters were considered to be constant, i.e.:
(  m
C0 237:6 __0 ; H < 1:2
¼ ð10Þ
h Eq: 7; H P 1:2
(  m
474 __0 ; H < 1:2
S¼ ð11Þ
Eq: 8; H P 1:2
The cohesive parameters for H = 1.2 were also used for the triaxiality independent analyses.
The load was applied as a prescribed displacement at the upper edge of the specimen. In dynamic cases, the
speed of the load was 3.3 m/s. The load and displacements presented in the following diagrams are the values
calculated on the boundary of the model.
Fig. 13 shows the load–displacement curves for quasi-static solutions. The figure shows that ignoring tri-
axiality dependency of interface elements can lead a static crack growth analysis to be highly conservative,
because the prediction of the energy absorption is much lower than the constraint (stress triaxiality) dependent
case. This shows that the values of triaxiality increase to more than the initial value considered in the analyses.
Fig. 14 shows the load–displacement curves obtained from dynamic analyses. Since transient dynamic solu-
tion has been employed, inertia and elastic waves are inherent in the analyses. The existence of oscillations in
the curves is inevitable because of the existence of elastic waves. The time that elastic waves travel the length of
the specimen is around 20 ls; that is much shorter than the total time of fracture (the shortest time of fracture
is 131.8 ls, which is for the case of rate-dependent bulk material). The elastic waves can have two sources, one
is the dynamically applied load and the other is the wave induced by broken cohesive elements or, in other

16000

14000

12000

10000
Load (N)

8000

6000

Quasi–static, Cohesive triaxiality dependent

4000 Quasi–static

2000

0
0 0.2 0.4 0.6 0.8 1 1.2
Displacement (mm)

Fig. 13. The effect of constraint on the load–displacement behavior (static analysis).
2222 M. Anvari et al. / Engineering Fracture Mechanics 73 (2006) 2210–2228

Dynamic
Dynamic–Bulk material and cohesive rate dependent, Cohesive triaxiality dependent
16000 Dynamic–Bulk material and cohesive rate dependent
Dynamic–Bulk material rate dependent

14000

12000

10000
Load (N)

8000

6000

4000

2000

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Displacement (mm)

Fig. 14. The effect of strain rate and constraint on the load–displacement behavior (dynamic analysis).

words, the waves induced because of the relatively high speed of the crack growth. Because of the implicit inte-
gration used in the dynamic analyses, there is no severe limitation on the time increment chosen. Repeating
one of the simulations with very short maximum allowable time increments (0.02 ls instead of 24 ls) showed
that the averaged results are similar, but the rate of convergence is very low and the simulation is extremely
time consuming.
Dynamic simulation considering the rate dependency only for bulk material leads to the most conservative
assumption. Ignoring the effect of strain rate and triaxiality on the behavior of interface elements makes the
results unrealistic. Using rate- and triaxiality-dependent cohesive elements results in more energy absorption,
although it is less than the simulation which ignores all of these influences.
Fig. 15 shows the crack growth versus time for dynamic simulations. In all of the cases, the crack speed is
low in the beginning and then it increases to a somewhat steady state speed. The case with no rate dependency
is an exception in it seems that the crack speed is changing during the growth. The figure shows clearly that the
steady state speed of the crack growth is highest, about 1430 m/s, for the case in which only bulk material is
rate dependent. If an average speed is defined for the case with no rate dependency, the lowest crack growth
speed of around 290 m/s is obtained for this case. It should be noted that for the case where both bulk material
and cohesive zone are rate dependent, the crack growth initiation happens almost at the same time as for the
case with no rate dependency. This is in contrast to the investigations of Basu [46] who numerically showed
that strain-rate sensitivity plays a beneficial role on dynamic ductile fracture initiation. This contradiction is
related to ignoring triaxiality change even at the first stages of deformation. It is observed in Fig. 15 that when
triaxiality dependency is considered, crack growth initiation happens at 69 ls which is higher than the rate-
insensitive case. To check if this effect is due to the rate sensitivity of the material, a dynamic case was sim-
ulated with triaxiality-dependent cohesive elements while no rate sensitivity was considered. The analysis
showed that the crack growth initiation happened at 61.8 ls. This proves that rate sensitivity has postponed
crack growth initiation, although after a short time, crack growth speed is higher than the rate-insensitive case.
 M. Anvari et al. / Engineering Fracture Mechanics 73 (2006) 2210–2228 2223

Dynamic–Bulk material rate dependent ,

Dynamic–Bulk material and cohesive rate dependent
30 Dynamic
Dynamic–Bulk material and cohesive rate dependent, Cohesive triaxiality dependent

25

20
Crack Growth (mm)

15

10

0
40 60 80 100 120 140 160
Time (μs)

Fig. 15. The effect of strain rate and constraint on crack growth rate (dynamic loading).

It should also be noted that in the case where stress triaxiality change has been ignored, the crack growth
velocity becomes so high that the assumption of ductile mechanism of crack growth might be inaccurate. It
is worth noting that for the case in which only rate dependency of bulk material has been considered, the crack
speed is around 1400 m/s. In this case, the inertia effect on fracture energy is important, but it has been ignored
in the calculation. Regarding these discussions and back to Fig. 14, it is well understood why the load in the
triaxiality- and rate-dependent case is higher in the beginning and drops very fast after some time. This is qual-
itatively compatible with Freund’s finding [45]. He used a rate-sensitive cohesive zone and showed analytically,
with some assumptions on the amount of plasticity and material behavior, that for high speed of crack growth
(more than about 20% of shear wave propagation velocity), ‘‘as the crack tip speed increases for some fixed
level of viscosity (rate sensitivity), the local stress is elevated and a smaller applied stress intensity factor is
required for crack growth’’. He also showed that for lower crack speeds or smaller viscosities, however, the
criteria is reversed.
The reason why crack speed is much higher in the case of rate-dependent plasticity is that the material hard-
ens due to high loading rate and opens the crack (cohesive elements) more easily. Fig. 16 shows the ratio of
dissipated plastic energy to the external work versus crack growth. This figure shows that for the rate-sensitive
material, the percentage of the energy absorbed by the plastic deformation is less compared to the rate-insen-
sitive case and the difference increases with crack growth. When the cohesive zone is rate dependent, the cohe-
sive strength is increased at high loading rates and this causes more energy dissipated by plasticity in the
surrounding material. Although the stresses needed for crack propagation are higher compared to the case
with rate-independent cohesive parameters, the cohesive strength, S, is not high enough to decrease the crack
growth rate as much as for the rate-insensitive case. The reason is that with increasing plasticity, the crack
opening speed decreases, so the local strain rate and consequently the cohesive strength decrease again.
Besides, stress triaxiality increases initially, but because of inertia it decreases while the crack is growing
[47]. This means that the decrease of the cohesive strength is related to the inertia. Fig. 17 shows the change
2224 M. Anvari et al. / Engineering Fracture Mechanics 73 (2006) 2210–2228

Dynamic
90 Dynamic–Bulk material and cohesive rate dependent
Dynamic–Bulk material rate dependent

80

70
Plastic Energy Ratio (%)

60

50

40

30

20

10

0
0 5 10 15 20 25 30 35 40 45 50
Crack Growth (mm)

Fig. 16. The effect of strain-rate sensitivity on percentage of energy dissipated by plastic deformation during crack growth.

1.8
Maximum Stress Triaxiality

1.7

1.6

1.5

1.4

1.3
0 2 4 6 8 10 12 14 16
Crack Growth (mm)

Fig. 17. The change of stress triaxiality during crack growth under dynamic loading conditions.
 M. Anvari et al. / Engineering Fracture Mechanics 73 (2006) 2210–2228 2225

Quasi–static–Cohesive triaxiality dependent

Dynamic–Bulk material and cohesive rate dependent, Cohesive triaxiality dependent
16000 Dynamic

14000

12000

10000
Load (N)

8000

6000

4000

2000

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Displacement (mm)

Fig. 18. A general comparison between quasi-static and dynamic cases.

2
Element size: 0.05 x 0.05 mm
2
Element size: 0.1 x 0.1 mm
16000 22

20
14000
18
12000
16
Crack Growth (mm)

10000 14
Load (N)

12
8000
10

6000 8

6
4000
4
2000
2

0 0
0 0.2 0.4 0.6 0.8 0 50 100 150 200
(a) Displacement (mm) (b) Time (μs)

Fig. 19. The effect of element size on (a) load–displacement curve, (b) crack growth.
2226 M. Anvari et al. / Engineering Fracture Mechanics 73 (2006) 2210–2228

of triaxiality value during crack growth for the case in which rate dependency and triaxiality have been con-
sidered for the cohesive elements.
Fig. 18 shows how the phenomena mentioned in the introduction part affect the mechanical response of the
structure and makes it clear that the effect of the combination of them is not easily predicted. It demonstrates
that although rate dependency and inertia can increase the maximum load in a displacement-controlled load-
ing, rate dependency decreases the total energy absorption of the specimen. This is in contrast to an uncracked
specimen in which rate-dependent plasticity and inertia increase the energy absorption under dynamic loading.
It was already mentioned that models with cohesive elements are mesh independent, but since an average of
field variables is calculated in the adjacent continuum element, this average may depend on the element size.
To check the influence of rate-sensitive cohesive elements on the mesh dependency of the model, the specimen
was also modeled by element size of 0.05 · 0.05 mm2 for the continuum elements adjacent to the ligament.
Fig. 19(a) shows the load–displacement curves for the two models. It is observed that the effect is insignificant
and only exists in the last part of the deformation. Fig. 19(b) shows that the crack growth rates are similar up
to some time and then the crack growth speed is slightly higher in the model with finer mesh. With the element
sizes considered in the present models, mesh density is not very important, especially for the global energy
absorption predicted. It should be noted that although cohesive zone models are theoretically mesh indepen-
dent, the mesh must not be so coarse that the calculation of local plastic strain field is imprecise or mesh
dependent [48].

4. Conclusion

The effect of triaxiality and strain rate on cohesive properties has been considered simultaneously in a finite
element dynamic crack growth simulation. Cohesive strength increases with strain rate and this implies an
increase in stress triaxiality. Taking into account the interaction between these influences to solve dynamic
crack growth problems in a closed form is not an easy task to do. The change of cohesive properties with
these variables has been obtained through calculations on a cell model obeying rate-dependent Gurson con-
stitutive equation. The results of the calculations have been used in crack growth simulations of a center-
cracked specimen under static and dynamic loading. The procedure has the ability to take into account
the effect of stress triaxiality, strain rate and inertia in a proper and relatively simple way. The general results
are:

(1) In a quasi-static loading, triaxiality increases with crack growth, but in dynamic cases it increases ini-
tially whereas inertia leads to its decrease during crack growth.
(2) Although strain-rate sensitivity makes a ductile crack initiation to be postponed, it leads to a faster crack
growth due to a decrease in the amount of plasticity at the crack tip area. In other words, although the
energy absorption increases initially, it drops very fast after a short crack growth. This is in contrast with
a dynamically loaded uncracked specimen, where positive strain-rate sensitivity always makes the spec-
imen absorb more energy.
(3) In the dynamic analyses performed, the presumption of ductile fracture will be inaccurate if the change
of triaxiality is ignored.
(4) Considering strain rate in the dynamic simulations while ignoring the stress triaxiality leads to a high
underestimation of the toughness.
(5) The global energy absorption of the specimen with rate-sensitive cohesive elements is not highly mesh
dependent. This is valid as long as the mesh is fine enough to calculate plastic stress distribution
correctly.

These results are compatible with other studies in this field which have used different approaches.

Acknowledgements

The authors from NTNU would like to acknowledge Hydro Aluminium and the Research Council of Nor-
way for their support through the NorLight project. The first author appreciates the scientific and friendly
 M. Anvari et al. / Engineering Fracture Mechanics 73 (2006) 2210–2228 2227

atmosphere in WMS department at GKSS research center during his stay. He also gives his special thanks to
professor Wolfgang Brocks for his precious support and advises.

References

[1] Nishioka T. Computational dynamic fracture mechanics. Int J Fract 1997;86:127–59.

[2] Siegmund T, Brocks W. A numerical study on the correlation between the work of separation and the dissipation rate in ductile
fracture. Engng Fract Mech 2000;67:139–54.
[3] Tvergaard V, Needleman A. Analysis of the cup-cone fracture in a round tensile bar. Acta Metall 1984;32:157–69.
[4] Zhang ZL, Thaulow C, Ødegård J. A complete Gurson model approach for ductile fracture. Engng Fract Mech 2000;67:155–68.
[5] Thomason PF. Ductile fracture of metals. Oxford: Pergamon Press; 1990.
[6] Dugdale DS. Yielding of steel sheets containing slits. J Mech Phys Solids 1960;8:100–4.
[7] Barenblatt GI. The mathematical theory of equilibrium cracks in brittle fracture. Adv Appl Mech 1962;7:55–129.
[8] Hillerborg A, Modéer M, Petersson PE. Analysis of crack formation and crack growth in concrete by means of fracture mechanics
and finite elements. Cem Concr Res 1976;6:773–82.
[9] Petersson PE. Crack growth and development of fracture zones in plain concrete and similar materials. Technical report, Report
LUTVDG/TVBM-1006, Lund Institute of Technology. 1981.
[10] Needleman A. A continuum model for void nucleation by inclusion debonding. J Appl Mech 1987;54:525–31.
[11] Tvergaard V. Material failure by void growth to coalescence. Adv Appl Mech 1990;27:83–151.
[12] Siegmund T, Brocks W. Predictions of the work of separation and implications to modeling. Int J Fract 1999;99:97–116.
[13] Tvergaard V, Hutchinson JW. The relation between crack growth resistance and fracture process parameters in elastic–plastic solids.
J Mech Phys Solids 1992;40:1377–97.
[14] Brocks W, Sun DZ, Hönig A. Verification of the transferability of micromechanical parameters by cell model calculations for visco-
plastic materials. Int J Plast 1995;11:971–89.
[15] Broberg KB. The cell model of materials. Comput Mech 1997;19:447–52.
[16] Tvergaard V. Crack growth predictions by cohesive zone model for ductile fracture. J Mech Phys Solids 2001;49:2191–207.
[17] Costanzo F, Walton JR. Numerical simulations of a dynamically propagating crack with a nonlinear cohesive zone. Int J Fract
1998;91:373–89.
[18] Corigliano A, Ricci M. Rate-dependent interface models: formulation and numerical applications. Int J Solids Struct 2001;38:547–76.
[19] Zhang X, Mai YW, Jeffrey RG. A cohesive plastic and damage zone model for dynamic crack growth in rate-dependent materials. Int
J Solids Struct 2003;40:5819–37.
[20] Broberg KB. Influence of T-stress, cohesive strength and yield strength on the competition between decohesion and plastic flow in a
crack edge vicinity. Int J Fract 1999;100:133–42.
[21] Wnuk MP, Legat J. Work of fracture and cohesive stress distribution resulting from triaxiality dependent cohesive zone model. Int J
Fract 2002;114:29–46.
[22] Tvergaard V, Hutchinson JW. Effect of strain-dependent cohesive zone model on predictions of crack growth resistance. Int J Solids
Struct 1996;33(20–22):3297–308.
[23] ABAQUS. ABAQUS version 6.4. H.K.S. Inc. Pawtucket, USA. 2003.
[24] Zhang ZL. A practical micro-mechanical model-based local approach methodology for the analysis of ductile fracture of welded T-
joints. PhD thesis. Lappeenranta University of Technology, Finland. 1994.
[25] Scheider I. Cohesive model for crack propagation analyses of structures with elastic–plastic material behavior. Foundations and
implementation. Technical report, GKSS internal report no. WMS/2000/19. 2000.
[26] Clausen AH, Auested T. Split-Hopkinson tension bar, experimental set-up and theoretical considerations. Technical report, SIMLab
internal report no. R-16-02. 2002.
[27] Pan J, Saje M, Needleman A. Localization of deformation in rate sensitive porous plastic solids. Int J Fract 1983;21:261–78.
[28] Gurson J. 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 1977;99:2–15.
[29] Tvergaard V. On localization in ductile materials containing spherical voids. Int J Fract 1982;18:237–52.
[30] Koplik J, Needleman A. Void growth and coalescence in porous plastic solids. Int J Solids Struct 1988;24:835–53.
[31] Needleman A, Tvergaard V. An analysis of dynamic, ductile crack growth in a double edge cracked specimen. Int J Fract
1991;49:41–67.
[32] Sun DZ, Hönig A, Böhme W, Schmitt W. Application of micromechanical models to the analysis of ductile fracture under dynamic
loading. In: Erdogan F, Hartranft RJ, editors. Fracture mechanics: 25th syp. ASTM STP, vol. 1220. Philadelphia: American Society
for Testing and Materials; 1994. p. 343–57.
[33] Curran DR, Seaman L, Shockey DA. Dynamic failure of solids. Phys Rep 1987;147:253–388.
[34] Johnson JW. Dynamic fracture and spallation in ductile solids. J Appl Phys 1981;52(4):2812–25.
[35] Liu J. A new dynamic void growth model. Doctoral thesis. Norwegian University of Science and Technology, Norway. 2004.
[36] Bernauer G, Brocks W. Micro-mechanical modelling of ductile damage and tearing – results of a european numerical round robin.
Fatigue Fract Engng Mater Struct 2002;25:384–463.
[37] Needleman A. An analysis of decohesion along an imperfect interface. Int J Fract 1990;42:21–40.
[38] Guinea GV, Planas J, Elices M. A general bilinear fit for the softening curve of the concrete. Mater Struct 1994;27:99–105.
2228 M. Anvari et al. / Engineering Fracture Mechanics 73 (2006) 2210–2228

[39] Scheider I, Brocks W. Simulation of cup-cone fracture using the cohesive model. Engng Fract Mech 2003;70:1943–61.
[40] Siegmund T, Needleman A. A numerical study of dynamic crack growth in elastic-viscoplastic solids. Int J Solids Struct
1997;34:769–87.
[41] Lin G, Cornec A, Schwalbe KH. Three-dimensional finite element simulation of crack extension in aluminum alloy 2024-FC. Fatigue
Fract Engng Mater Struct 1998;21:1159–73.
[42] Scheider I, Brocks W. The effect of the traction separation law on the results of cohesive zone crack propagation analyses. Key Engng
Mater 2003;251–252:313–8.
[43] Falk M, Needleman A, Rice JR. A critical evaluation of cohesive zone models of dynamic fracture. J Phys IV France 2001;11:43–50.
[44] Négre P, Steglich D, Brocks W. Crack extension in aluminium welds: a numerical approach using the Gurson–Tvergaard–Needleman
model. Engng Fract Mech 2004;71(71):2365–83.
[45] Freund LB. Dynamic fracture mechanics. Cambridge University Press; 1998.
[46] Basu S, Narashimhan R. Finite element simulation of mode I dynamic ductile fracture initiation. Int J Solids Struct
1996;33(8):1191–207.
[47] Basu S, Narashimhan R. A numerical investigation of loss of crack tip constraint in a dynamically loaded ductile specimen. J Mech
Phys Solids 2000;48:1967–85.
[48] Needleman A. Numerical modeling of crack growth under dynamic loading conditions. Comput Mech 1997;19:463–9.