3D CAFE Modelling of Transitional Ductile - Brittle Fracture in Steels

Department of Mechanical Engineering
The University of Sheeld

3D CAFE modelling of transitional
ductile brittle fracture in steels
A thesis submitted for the degree of Doctor of
Philosophy
by
Anton Shterenlikht
September 2003
Ingenious modications. . . cannot change the basic error of the Berg-Gurson
approach. . .
P F Thomason. Ductile Fracture of Metals, Pergamon Press, 1990
3
Contents
Acknowledgements 7
Summary 8
Nomenclature 10
1 The problem 13
2 Solutions 15
2.1 Microanalysis of ductile fracture . . . . . . . . . . . . . . . . . . 16
2.1.1 McClintock model . . . . . . . . . . . . . . . . . . . . . . 17
2.1.2 Rice-Tracey model . . . . . . . . . . . . . . . . . . . . . . 18
2.1.3 Argon-Im-Safoglu model . . . . . . . . . . . . . . . . . . . 19
2.1.4 Berg-Gurson-Tvergaard-Needleman model . . . . . . . . . 19
2.1.5 Lemaitre model . . . . . . . . . . . . . . . . . . . . . . . . 22
2.1.6 Rousselier model . . . . . . . . . . . . . . . . . . . . . . . 23
2.1.7 Thomason model . . . . . . . . . . . . . . . . . . . . . . . 24
2.1.8 Cavitation models . . . . . . . . . . . . . . . . . . . . . . 26
2.2 Microanalysis of brittle fracture . . . . . . . . . . . . . . . . . . . 26
2.2.1 Crack initiation models . . . . . . . . . . . . . . . . . . . 26
2.2.2 Weakest link models . . . . . . . . . . . . . . . . . . . . . 28
2.2.3 Crack arrest . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.3 Coupled ductile-brittle fracture modelling . . . . . . . . . . . . . 30
2.3.1 Size scales . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.3.2 Brittle fracture as a postprocessing operation . . . . . . . 32
2.3.3 Folch model . . . . . . . . . . . . . . . . . . . . . . . . . . 32
5
6 CONTENTS
2.4 Model calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3 The CAFE solution 37
3.1 A CAFE model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.2 The full model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.2.1 The ductile CA array . . . . . . . . . . . . . . . . . . . . 44
3.2.2 The brittle CA array . . . . . . . . . . . . . . . . . . . . . 44
3.2.3 The FE part . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.2.4 How the model works . . . . . . . . . . . . . . . . . . . . 48
3.2.5 Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.3 The simplied model . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.3.1 How the model works . . . . . . . . . . . . . . . . . . . . 55
3.4 Important properties . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.5 The list of model parameters . . . . . . . . . . . . . . . . . . . . 57
4 Results 59
4.1 The full CAFE model . . . . . . . . . . . . . . . . . . . . . . . . 59
4.1.1 Single FE, tension compression . . . . . . . . . . . . . . 59
4.1.2 Single FE, forward tension simulation of scatter . . . . 69
4.1.3 The Charpy test . . . . . . . . . . . . . . . . . . . . . . . 79
4.2 The simplied CAFE model . . . . . . . . . . . . . . . . . . . . . 95
4.2.1 The Charpy test . . . . . . . . . . . . . . . . . . . . . . . 95
5 Discussion 111
5.1 Unresolved problems and future work . . . . . . . . . . . . . . . 115
5.2 Overall conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . 116
A The CA cell neighbourhood 119
B Rousselier model integration 121
Bibliography 127
Acknowledgements
I would like to acknowledge my academic advisor Professor Ian C Howard for
encouragement, support and mind-stimulating discussions.
I would also like to acknowledge Dr C Davis and Dr D Bhattacharjee, both
from The University of Birmingham, Metallurgy and Materials for the permis-
sion to use their unpublished data (Dr D Bhattacharjee is now with Tata Iron
& Steel Co., Jamshedpur, India).
I would also like to acknowledge Corus UK Ltd for the provision of broken
Charpy samples and corresponding test data.
7
Summary
A coupled Cellular Automata Finite Element (CAFE) three-dimensional multi-
scale model was applied in this work to the simulation of transitional ductile-
brittle fracture in steels. In this model material behaviour is separated from the
representation of structural response and material data is stored in an appropri-
ate number of cellular automata (CA). Two CA arrays, the ductile and the
brittle, are created, one is to represent material ductile properties, another is
to account for the brittle fracture. The cell sizes in both arrays are independent
of each other and of the nite element (FE) size. The latter is chosen only to
represent accurately the macro strain gradients. The cell sizes in each CA array
are linked to a microstructural feature relevant to each of the two fracture mech-
anisms. Such structure of the CAFE model results in a dramatic decrease of the
number of nite elements required to simulated the damage zone. Accordingly
the running times are cut down signicantly compared with the conventional FE
modelling of fracture for similar representation of microstructure. The Rousse-
lier continuing damage model was applied to each cell in the ductile CA array.
The critical value of the maximum principal stress was used to assess the failure
of each cell in the brittle CA array. The model was implemented through a
user material subroutine for the Abaqus nite element code. Several examples
of model performance are given. Among them are the results of the modelling
of the Charpy test at transitional temperatures. For a laboratory rolled TMCR
steel the model was able to predict the transitional curve in terms of the Charpy
energy and the percentage of brittle phase, including realistic levels of scatter,
and the appearance of the Charpy fracture surface. The ways in which material
data can be tted into the model are discussed and particular attention is drawn
upon the signicance of the fracture stress distribution.
9
Nomenclature
In this work tensor analysis is used whenever possible. The tensor quantities
are given as in Kachanov (1971).
Many symbols might have various sub- and superscripts. These are described
in the text.
grain orientation angle
damage variable (Rousselier model)
cell solution-dependent variable
change in variable during one time increment
ij
Kronecker delta
ij
strain tensor
e
ij
elastic strain tensor
p
ij
plastic strain tensor,
p
ij
= e
p
ij
+
p
m
ij
e
p
ij
plastic strain deviator
p
m
mean plastic strain,
p
m
=
1
3
ii
p
eq
equivalent plastic strain,
p
eq
=
_
2
3
e
p
ij
e
p
ij
the fraction of the brittle CA cells which have a grain boundary carbide
F
misorientation threshold
cell property
Poissons ratio
CA to FE transition function
ij
stress tensor,
ij
= S
ij
+
m
ij
m
mean stress,
m
=
1
3
ii
eq
equivalent stress,
eq
=
_
3
2
S
ij
S
ij
I
maximum principal stress
F
fracture stress
11
12 NOMENCLATURE
Y
yield stress
Y0
rst yield stress
cell state
cell transitional rule
A total number of state variables per FE integration point
C
v
the total energy absorbed in the Charpy V-notch impact test
c concentration factor for a CA array
d
g
grain size
d
k
direction cosines
E Youngs modulus
E
ijkl
isotropic elastic modulus tensor, E
ijkl
= 2G
ik
jl
+
_
K
2
3
G
_
ij
kl
f a probability density function
f
0
initial void volume fraction (Rousselier model)
G shear modulus, G =
E
2(1+)
K compression modulus, K =
E
3(12)
L damage cell size
L
FE
nite element size
M mapping nction
M total number of cells per CA
N the set of natural numbers
N total number of cell properties
n hardening exponent
Q total number of cell state variables
R total number of integration points per nite element
S
ij
stress deviator
t time
T temperature, in

C
W
shape parameter of Weibull dustribution

W
location parameter of Weibull dustribution

W
scale parameter of Weibull dustribution

X
max
the maximum number of dead cells allowed per CA
Y nite element solution-dependent variable
Chapter 1
The problem
There are two fundamental problems in modelling transitional ductile-brittle
fracture with nite element analysis. Both problems have their roots in the
complex inhomogeneous nature of materials such as steels and in the limitations
of the nite element approach. The rst problem is the high computational cost
due to large numbers of small nite elements. Conicting demands for the
mesh size due to the dierent physical nature of ductile and brittle fracture is
the second.
The local approach to fracture is a technique suitable for fracture propaga-
tion modelling because it takes into account only a small area ahead of the crack
tip. Therefore this approach is geometry-independent as opposed to single- and
two-parameter methods of fracture mechanics.
Exactly how small this area should be is determined by the need to correctly
represent the stress and strain gradients ahead of the notch tip. The stress and
strain elds there are the result of a complex interaction of dierent microstruc-
tural features. These can be grains, grain clusters, lath packets (in martensitic
and bainitic steels), large and small inclusions, grain boundary carbides, larger
precipitates, microcracks and microvoids etc. One common feature of all entries
in the above list is their size they are all small compared to any structure
of engineering interest. Thus a nite element mesh of a structure with a crack
must have a highly rened region extending long enough ahead of the crack
tip to allow for modelling of the desired crack advance. In practice meshes
13
14 CHAPTER 1. THE PROBLEM
with tens of thousands of nite elements are not uncommon. The analysis of
such meshes takes weeks or months and is very unstable due to ill-conditioned
stiness matrices.
At the same time the microstructural objects themselves can dier in size,
e.g. a grain is typically tens of times larger than a grain boundary carbide
and tens of times smaller than a lath packet. This has a profound inuence
on the ruling mesh size designed for the analysis of brittle or ductile fracture
because the fracture progresses in microstructurally sensitive steps. In the case
of ductile fracture these steps will usually be of the order of spacing between
the microvoids or large inclusions. Grains, lath packet or a group of grains with
small misorientation angles are the objects whose sizes are usually taken as a
basis for the steps of brittle fracture advance. As the above step sizes might
dier tens of times so do the mesh sizes required to simulate the propagation of
brittle or ductile fracture. The only way these conicting requirements can be
satised within a single nite element mesh is by choosing a compromise mesh
size. The accuracy of the solution is then a question.
The above two fundamental problems exist because in conventional nite
element analysis a nite element is a material and a structural unit simultane-
ously. The structure and material are thus merged into an inseparable entity.
This approach can be very ineective.
The Cellular Automata Finite Element (CAFE) approach used in this work
oers solution to both problems mentioned above. In this approach material
properties are moved away from the nite element mesh and distributed across
the appropriate number of cellular automata arrays. Thus a nite element
mesh is designed only to represent the macro strain gradients adequately. This
is now a solely structural entity. A number of cellular automata arrays, in which
cell sizes can be chosen independently, provide the means to analyse material
properties at each size scale separately. So a CAFE model can accommodate
as many size scales as necessary to address all material properties of interest.
However only two cellular automata arrays are required to model the transitional
ductile-brittle fracture.
The following chapter leads to the formulation of the CAFE model starting
with a review of major models for ductile, brittle and ductile-brittle fracture
proposed during the last half-century or so.
Chapter 2
Solutions
The fact that materials have a complex microstructure has long been recognised
by materials engineers and scientists (Czochralski, 1924; Nadai, 1950; Cottrell,
1967; Gilman, 1969). In fact, had a material been homogeneous, it would be
perfectly elastic until the nal rupture by the separation of atoms (Knott, 1973;
Thompson and Knott, 1993; Hertzberg, 1996). This case would be perfectly
described by a single critical parameter, fracture toughness. It is the existence
of grains, grain boundaries, inclusions or, on even lower level, dislocations, that
demands the use of more complicated approaches to fracture analysis.
Extensive experimental studies of macro and micro fracture mechanisms
resulted in understanding of two distinctive failure physical processes. The rst
is broadly called ductile and is characterised by relatively high energy needed
for fracture to take place, high level of macro plasticity and dull appearance
of the fracture surface. The fracture process that requires much less energy,
produces bright, light-reective fracture surfaces and accompanied by little or
no plasticity is commonly called brittle. This is the second type of fracture.
Exactly how these two processes take place on a micro scale has been one
of the main issues of experimental research in fracture mechanics for the last
three decades. Simultaneously a number of material models describing the ex-
perimental ndings have been developed.
15
16 CHAPTER 2. SOLUTIONS
2.1 Microanalysis of ductile fracture
A number of authors have observed regions of increased porosity next to the
fracture surfaces in ductile metals (Tipper, 1949; Puttick, 1959; Rogers, 1960;
Beachem, 1963; Gurland and Plateau, 1963; Bluhm and Morrissey, 1966; Liu
and Gurland, 1968; Hayden and Floreen, 1969; Gladman et al., 1970, 1971;
Gurland, 1972; Goods and Brown, 1979). Rhines (1961) was able to reproduce
the observed porosity in plasticine using polystyrene spheres as inclusions.
Therefore it was proposed that ductile fracture in steels is fracture by the
growth of holes (McClintock, 1968), ductile fracture by internal necking of
cavities (Thomason, 1968), is caused by the large growth and coalescence
of microscopic voids (Rice and Tracey, 1969) and is via the nucleation and
growth of voids (Gurson, 1977a). Long before, Bridgman (1952) came to sim-
ilar conclusions analysing the inuence of hydrostatic pressure on the necking
behaviour in tensile tests. He found that ductility is increasing with increased
pressure up to a point where no cup-and-cone fracture can be observed and
the diameter of the neck is approaching zero. Bridgman explained this by the
closure of voids under very high pressure. Beachem (1975) reported eight (and
predicted another possible six) types of dimple shapes tied to a fracture mode.
Ductile fracture by void growth and coalescence involves three stages: mi-
crovoid nucleation, void growth and void coalescence (Bates, 1984; Thomason,
1990; Gladman, 1997; Thomason, 1998).
Voids might nucleate at cleaved particles (Gladman et al., 1971; Cox and
Low, 1974) or by decohesion of the interfaces of the second phase particles
(Beachem, 1975; Argon et al., 1975; Argon and Im, 1975). Smaller particles
require higher applied stresses for decohesion than larger ones. Based on this
Bates (1984) showed that although carbides play secondary role in tensile test
fracture, they might dominate the fracture process in fracture toughness test.
Void growth can be dilatational (volumetric) or by shape change. Stress
triaxiality has a dramatic eect on void growth type and therefore on strain to
fracture. In a tensile test voids grow in the direction of tensile stress prior neck-
ing. The onset of necking changes the uniaxial stress state to triaxial (Bridgman,
1952) which causes some volumetric growth (Gladman, 1997; Thomason, 1998)
and therefore signicantly lowers the strain to rupture.
2.1. MICROANALYSIS OF DUCTILE FRACTURE 17
Void coalescence is a process involving a localised internal necking of the
intervoid material (Thomason, 1981) and was observed in dierent materials
by Puttick (1959); Rhines (1961); Bluhm and Morrissey (1964) (very impressive
photographs from these works were reprinted by McClintock (1968) and Thoma-
son (1990, 1998)). The nal stages of this process are associated with the failure
of the submicron intervoid ligament by shearing along crystallographic planes
or by microcleavage (Rogers, 1960; Cox and Low, 1974).
Development of theoretical models was slow due to the complex nature of
ductile fracture phenomena. Three stages (nucleation, growth and coalescence
of voids) have dierent nature and require separate physical models. As noted
by McClintock (1968), contrary to the initial yielding or brittle fracture, where
only the present stress state is needed for analysis, the size, shape and spacing
of holes are a result of the whole history of straining. Some of the major models
for ductile fracture are described below.
2.1.1 McClintock model
McClintock (1968) proposed a model for void growth and derived a criterion for
ductile fracture. He assumed a material containing a regular three-dimensional
array of cylindrical voids of elliptical section. The main axes of this array are
parallel to the principal stress axes. The condition for fracture was that each
void touches the neighbouring one. If the voids have the cylindrical axes parallel
to the z direction and two semi axes are designated as a and b, and if the voids
grow in the b direction then the approximate expression for the onset of fracture
takes the form:
d
zb
d
eq
=
1
lnF
f
zb
_

3
2(1 n)
sinh
_
3(1 n)
2
a
+
b
eq
_
+
3
4
eq
_
(2.1)
where
d
zb
deq
is a damage rate (
eq
- equivalent strain, d
zb
- damage increment),
F
f
zb
is a critical value of the relative growth factor, n is a hardening exponent,
a
and
b
are two of the principal stresses at innity and
eq
is the equivalent
stress.
The over-simplied nature of this model leads to unrealistic results. Most
important is that according to this model void growth is a smooth process
until the nal rupture, whereas, as argued by Thomason (1968, 1981, 1998) and
observed by Liu and Gurland (1968) and Hayden and Floreen (1969), the onset
of failure by void coalescence is essentially due to a loss of stability.
Nevertheless even this simple model demonstrates some fundamental fea-
tures of ductile fracture, e. g. very strong decrease of failure strain with increase
of stress triaxiality and a size eect, the need to know the stress history over
a region of the order of the void spacing.
2.1.2 Rice-Tracey model
The approach undertaken by Rice and Tracey (1969) is based on variational
analysis and the principle of maximum plastic work (Hill, 1983; Prager, 1959)
or Druckers stability postulate (Drucker, 1951, 1959; Khan and Huang, 1995).
The authors analysed a case of dilatational growth of a single spherical void in a
material under uniform stress state applied at innity. They derived a classical
equation for void enlargement under a high triaxiality stress state:
D = 0.283 exp
_
1.5
eq
_
(2.2)
where D is the ratio of the strain rate on the surface of a void to the strain rate
at innity,
m
is the mean stress and
eq
is the equivalent stress.
The simplicity of the resulting equation is the major advantage of this model.
Probably it is for simplicity that it is by far the most famous void growth related
equation.
The practical use of this equation however is quite limited because the model
does not address void interaction, it cannot predict the fracture strain and
cannot explain ductile failure in pure shear. Indeed according to the equation
(2.2) if
m
= 0 then the void acts merely as a stress concentrator with a constant
concentration factor.
Finally as pointed out by Thomason (1990) and Gladman (1997) void exten-
sion can be found only at very high levels of negative hydrostatic pressure. Void
shape distortion has a much bigger contribution in the process of void growth
(Liu and Gurland, 1968; Hayden and Floreen, 1969).
Needleman (1972) applied similar variational analysis to a doubly periodic
square array of circular cylindrical voids under plain strain conditions. He used
FEA to minimise the resulting functional. This model is not particularly famous,
however, it inspired Gurson, section 2.1.4.
2.1.3 Argon-Im-Safoglu model
The authors analysed only the rst stage of the process, i.e. void nucleation.
They proposed that the criterion for separation of large particles is reaching lo-
cally of a critical interfacial tensile strength. For the case of spherical inclusions
they derived the following equations for the radial stresses on the inclusion-
matrix interface (Argon et al., 1975).
For non-interacting inclusions:
rr
= k
0
_
_
_

y
_1
n
+
3
_
6 (n + 1)
m
y
_ 1
n+1
_
_
(2.3)
For interacting inclusions:
rr
= k
0
_
_
3
_
y
_1
n
+
6
m
+
_

y
_1
n
_
_
(2.4)
where k
0
is the yield stress in shear, and
y
are the present and yield shear
strains, n is a hardening exponent, m is the Taylor factor, generally taken as
3.1, is the inter-particle spacing and is the particle radius.
According to this approach decohesion occurs when
rr

c
(2.5)
In a companion paper (Argon and Im, 1975) the authors obtained
c
experi-
mentally for several materials and found values from
c
= 990 MPa (Cu-0.6Cr
alloy) to
c
= 1820 MPa (martensitic steel).
2.1.4 Berg-Gurson-Tvergaard-Needleman model
Berg (1970) suggested that localisation occurs when the hardening behaviour of
the matrix material is overweighed by softening due to the dilation of voids.
Inspired by the works of Berg (1970), McClintock (1968), Rice and Tracey
(1969) and Needleman (1972), Gurson (1977a) proposed a methodology for
obtaining an approximate yield surface for a material containing voids. He
applied a maximum plastic work principle to kinematically admissible velocity
elds (Nadai, 1950; Drucker, 1959; Prager, 1959; Hill, 1983; Kachanov, 1971)
for long circular cylindrical and spherical voids. For the latter case the yield
function had the form:
=
_
eq
y
_
2
+ 2fcosh
_
3
m
2
y
_
1 f
2
= 0 (2.6)
where f is a void volume fraction. This condition is reduced to the classical
Mises yield criterion if f = 0.
The change in void volume fraction was described as:
f =

f
g
+

f
n
(2.7)
where

f is the void volume fraction rate,

f
g
is the rate of growth of existing
voids and

f
n
is the void nucleation rate. For the growth of existing voids Gurson
(1977b) proposed only dilation:
f
g
= (1 f)
p
ij
I
ij
(2.8)
where
p
ij
is a plastic strain rate and I
ij
is the second-order unit tensor.
Various nucleation models have been proposed by Gurson (1977b), e.g.:
f
n
= f
(
p
eq
)
p
eq
(2.9)
where f
- is the void nucleation intensity and

p
eq
- is equivalent plastic strain.
Dierent aspects of void nucleation are discussed in Zhang et al. (2000);
Tvergaard (1990). In materials containing large inclusions, e.g. MnS particles,
voids would typically nucleate from these particles at the beginning of the plastic
deformation. For modelling purposes it is reasonable to assume that all voids are
nucleated at the beginning of the simulation and their amount is described by
a single parameter the initial void volume fraction, f
0
, (Zhang et al., 2000).
If, however, the material is such that voids are mostly nucleated from small
particles, typically carbides, than a continuous nucleation mechanism is more
appropriate (HKS, 2001).
When Yamamoto (1978) applied the yield function of equation (2.6) to the
analysis of a shear band following a localisation theory of Rice (1977), he found
that for a body without imperfection this yield condition predicts unrealistically
large strain at localisation. He concluded that initial imperfections are necessary
in order to achieve localisation at reasonable strain.
In an attempt to reduce the above discrepancy Tvergaard (1981, 1982a,b)
introduced two adjustable parameters, q
1
and q
2
in Gursons yield function (2.6):
=
_
eq
y
_
2
+ 2q
1
fcosh
_
3q
2
m
2
y
_
_
1 +
_
q
1
f
2
_
= 0 (2.10)
If q
1
= q
2
= 1 then (2.10) is reduced to (2.6). After comparison of modelling
results with those obtained experimentally by Gladman et al. (1970) and Glad-
man et al. (1971), Tvergaard (1981, 1982b) concluded that q
1
= 1.5 and q
2
= 1
improve the performance of the yield condition (2.6) by approximately 50%.
Some authors went further and introduced the q
3
parameter (HKS, 2001):
=
_
eq
y
_
2
+ 2q
1
fcosh
_
3q
2
m
2
y
_
_
1 +q
3
_
f
2
_
= 0 (2.11)
Although some authors argued that q
1
, q
2
and q
3
are true material constants
(Tvergaard, 1982b, 1990; HKS, 2001), there is a growing experimental evidence
that they depend on the triaxiality level (Pardoen and Hutchinson, 2000; An-
drews et al., 2002).
The initial Gurson model can only simulate void nucleation and dilation.
However it does not account for void coalescence in any way. Tvergaard and
Needleman (1984) introduced the function f
(f) to model the rapid loss of

stress-carrying capacity. This was an attempt to account for void coalescence.
The function f
(f) was chosen as:

f
(f) =
_
_
_
f for f f
c
f
c
u
f
c
f
f
f
c
(f f
c
) for f > f
c
where f
c
is a critical value of void volume fraction, the value at which a rapid
loss of load-bearing capacity begins; f
f
is void volume fraction at nal fracture
and f
u
= 1/q
1
.
Based on experimental (Brown and Embury, 1973; Goods and Brown, 1979)
and numerical (Andersson, 1977) results Tvergaard and Needleman (1984) have
chosen the values f
c
= 0.15 and f
f
= 0.25.
This model, which is usually called Gurson-Tvergaard-Needleman or sim-
ply GTN, is probably used most frequently in engineering applications. It is
included in commercial nite element packages (HKS, 2001). However the na-
ture of the model still makes it dicult to achieve a good correlation with
experiment. The main problem is that the model predicts strain to fracture
that is much higher than that observed in experiments.
Thomason (1981, 1985b); Zhang et al. (2000) and others argued that the
excessive high strain to fracture is a direct consequence of the model taking into
account only the homogeneous deformation. In fact the voids were eectively
substituted by a continuous porosity eld. The only eect of voids in this model
is through the pressure-dependent yield surface.
There are other attempts, apart from the modications introduced by Tver-
gaard and Needleman, to extend the validity of the Gurson (1977a) model. A
model dealing with prolate and oblate ellipsoidal cavities was proposed (Golo-
ganu et al., 1993, 1994). Based on this model and ideas of Thomason (section
2.1.7) Pardoen and Hutchinson (2000) introduced an extended model. Another
combination of Gursons and Thomasons approaches produced a complete
Gurson model (Zhang et al., 2000). These modied Gurson models produce
more realistic results than the GTN model. However they are signicantly
more complex.
2.1.5 Lemaitre model
Based on concepts of a damage variable, D, and eective stress, =

1D
,
(Kachanov, 1971; Rabotnov, 1969), Lemaitre (1985, 1996) proposed a thermo-
dynamically consistent (Ziegler, 1977; Germain et al., 1983) theory of damage.
The model assumes a representative volume of material containing defects
(microcracks or microvoids). If the intersection of this volume with a plane
dened by a normal vector n is S and the area of intersection of voids and cracks
of the volume by this plane is S
D
(n), then the damage variable is dened as
D(n) =
SD(n)
S
. Since 0S
D
(n) S then 0 D(n) 1. D = 0 means
undamaged material whereas D = 1 means that material has no load-bearing
capacity.
If D does not depend on n that the damage is considered isotropic and
D =
SD
S
. In this case the damage variable has the meaning of eective density
of microdefects.
The major principle of Lemaitres work is that any strain constitutive equa-
tion for a damaged material may be derived in the same way as for a virgin ma-
terial except that the usual stress is replaced by the eective stress (Lemaitre,
1996). Based on this principle he derived the constitutive equation for ductile
damage. For the case of isotropic damage this condition has the form:
D =
_
A
2
2ES
0
_
2
3
(1 +) + 3(1 2)
_
eq
_
2
_
_
p
eq
_ 2
n
_
s0

p
eq
(2.12)
where A and n are material properties in the Ramberg-Osgood hardening law:
p
=
_

A
_
n
, (2.13)
S
0
and s
0
are parameters in the damage evolution law:
D =
_
y
S
0
_
s0+1

p
eq
, (2.14)
which is based on the normality rule of potential of dissipation, :
D =
y
. (2.15)
In the last two equations y is called the damage strain energy release rate
(Lemaitre, 1985).
The Lemaitre model is very powerful in the sense that it can be applied to
any damage process, not just ductile damage. Its weak side, however, is that
by the very nature of thermodynamically consistent theory it is a continuum
theory. Therefore the Lemaitre ductile damage model is essentially a continuum
softening one where the presence of voids or cracks is introduced via damage
variable, D.
Other models based on Continuum Damage Mechanics (CDM) have been
formulated over the years (McDowell, 1997).
2.1.6 Rousselier model
Another thermodynamically consistent ductile damage theory was introduced
by Rousselier (1981). The plastic potential in this model has the form:
eq
H
_
p
eq
_
+B() Dexp
_
1
_
= 0 (2.16)
where:
=
p
eq
Dexp
_
1
_
(2.17)
() =
1
1 f
0
+f
0
exp
(2.18)
B() =

1
f
0
exp
1 f
0
+f
0
exp
, (2.19)
is a scalar damage variable. Its evolution is determined by equation 2.17.
While material is within elasticity limits = 0,
B is the damage function,
is a dimensionless density. From equation 2.18 it follows that decreases
with increasing ,
D and
1
are material constants,
f
0
is the initial void volume fraction and
H
_
p
eq
_
is a term describing the hardening properties of material. Usually
this is equal to the yield stress of the undamaged material, H
_
p
eq
_
=
Y
_
p
eq
_
.
The Rousselier model has the same strong and weak sides as the previous
two, GTN and Lemaitre (Rousselier, 1987). All three are continuum damage
models and can therefore be used as constitutive models for material with mi-
crocavities. The models can be used for numerical simulation of fracture propa-
gation. Their weak side, however, is the inability to model shear fractures since
only a volumetric void growth is allowed.
2.1.7 Thomason model
Thomason (1968, 1981, 1982, 1985a,b, 1993) studied the details of the coales-
cence phenomenon. He formulated a sucient condition for the stability of
incompressible plastic ow in the presence of microvoids for two- and three-
dimensional cases. His models based on plasticity theory (Hill, 1983; Kachanov,
1971) and theorems of limit analysis (Prager, 1959) were summarised in a book
(Thomason, 1990). He later criticised the models proposed by Gurson and
Rousselier as being based on an erroneous criterion of microvoid coalescence
(Thomason, 1998). Indeed Yamamoto (1978) has shown that Berg-Gurson
model gives realistic critical strains only after signicant changes in void volume
fractions.
Thomasons analysis resulted in the concept of incipient void coalescence
leading to an instantaneous change from incompressible to dilational plasticity
(Thomason, 1981). The condition for the onset of coalescence in a plane strain
case has the following form:
(
Ic
I
)
Ic
= 0 (2.20)
where
Ic
is the plastic limit-load stress,
I
is the maximum principal stress and

Ic
is the maximum principal strain rate across an intervoid matrix neck.
If only cylindrical voids are assumed then
Ic
can be represented by the
following empirical equation (Thomason, 1998):
Ic
2k
= 1.43 f
1/6
0.91 (2.21)
where k is the maximum shear stress, k =
I

m
, and f is a void volume
fraction.
The Thomason model for void coalescence is probably the most physically
based and accurate to the present day. A strain to failure in a uniaxial test
predicted with this model is in a good agreement with the experimental results
(Thomason, 1982). However it can hardly be used for numerical modelling
because of two fundamental problems.
The limit-load model is not a constitutive one. It can only predict the onset
of void coalescence as a start of material and, what is more important, structural
instability. Thus a structure is considered instantly failed when the condition
of equation (2.20) is met.
The material rate of hardening in the intervoid matrix approaching duc-
tile fracture is reduced to a very low level. Plastic solids with low work-
hardening rate are described by second-order hyperbolic partial dierential
equations (Thomason, 1998). These equations cannot at present be solved with
nite element methods (Johnson, 1987; Belytschko et al., 2000) but with what
mathematicians call the method of characteristics (Smith, 1965; Johnson, 1987)
and engineers the slip-line technique (Hill, 1983).
A similar approach to void coalescence problem was explored by Szczepi nski
(1982) who argued that a theoretical analysis of a plane strain rigid-plastic
material model with cylindrical holes is very complicated because there exists
strong stress concentration at the edges of the holes even during purely elastic
response. He proposed an idealised conguration where voids are initially intro-
duced as slits. The theoretical analysis can be then performed based on slip-line
technique.
2.1.8 Cavitation models
In most cases microvoids do not show signicant growth before the onset of
coalescence (Thomason, 1998). However if constraint is very high and if very
few void nucleation cites are present then volumetric void growth can be very
strong.
Ashby et al. (1989) observed the enlargement of a single void by a factor of
more than 10
6
in tensile tests of highly constrained lead wires.
Huang et al. (1991) analysed a single spherical void in elastic-plastic ma-
terials under a remote stress eld. They showed that a complex interaction
of elasticity and plastic yielding can lead to a cavitation instability, if the
stresses in the material surrounding the void are suciently high so that the
work produced by these stresses to expand the void is less than the energy
released by such expansion.
It is easy to draw an analogy between the above analysis of the cavitation
instability and the energy condition of Grith (1924) for an unstable crack
growth.
Faleskog and Shih (1997) conducted a two-dimensional plane strain nite
element analysis of a square material cell containing a single cylindrical void
in its centre. Their results were very similar to those reported by Huang et al.
(1991), that the stored elastic energy can cause void expansion by several orders
of magnitude over a negligible macroscopic strain increment.
2.2 Microanalysis of brittle fracture
Signicant advance has been made in understanding of the brittle fracture phe-
nomenon since Griths days (Grith, 1921, 1924). At present there is a vast
amount of literature on the subject. A number of review books and papers have
been published (e.g. Averbach et al., 1959; Knott, 1973; Hahn, 1984; Thompson
and Knott, 1993).
2.2.1 Crack initiation models
It is generally agreed that a dislocation pile-up at an obstacle, such as a grain
carbide interface, can cleave a grain boundary carbide and thus initiate a
2.2. MICROANALYSIS OF BRITTLE FRACTURE 27
microcrack (McMahon and Cohen, 1965; Lin et al., 1987; Thompson and Knott,
1993). Thus some degree of plastic deformation in a ferrite grain is always
necessary to fracture a neighbouring carbide (Lin et al., 1987).
The stresses required to generate a microcrack can be written as follows:
= 4.4

na
, (2.22)
= K

na
, (2.23)
where and are shear and normal stresses accordingly, is an eective surface
energy, n is the number of dislocations piled up against a grain boundary, a is
the atomic spacing and K is a coecient depending on the arrangement of the
dislocation pile-up.
Equation (2.22) was obtained by Zener in 1948 (Hahn et al., 1959). The
coecient K in equation (2.23) can be K = 2.7 (Orowan model, 1954), K = 5.3
(Bullough model, 1956) or K = 2 (Cottrell, 1959). The rst two values are
taken from Hahn et al. (1959). The exact value of K depends on the assumed
dislocation model. Zener analysed a crack forming on a plane normal to the
operative slip plane. Orowan suggested that a polygonised array of dislocations
can generate a crack in the slip plane. In the Bullough model the crack occurs
in the slip plane. Cottrell assumed that two intersecting (110) slip planes in
b.c.c. materials produce a microcrack on the common (100) plane (Hahn et al.,
1959).
The shape of equations (2.22) and (2.23) suggests that the number of cracked
carbides increases with applied strain. Indeed the number and intensity of dis-
location pile-ups increases with plastic straining and hence the stresses required
to generate a microcrack decrease. This point is supported by experimental
observations (Gurland, 1972).
A microcrack in a cleaved carbide can advance if the following condition is
met:
n

F
(2.24)
where
n
is a normal stress acting across the grain-carbide interface and
F
is
a fracture or cleavage stress.
Smith (1966a,b) derived an equation for the fracture stress of a carbide
ferrite interface. Based on Smiths analysis Lin et al. (1987) obtained a similar
equation for the fracture stress of a ferrite ferrite interface. Both equations
are shown below.
cf
F
=
E
cf
(1
2
) d
c
, (2.25)
ff
F
=
E
ff
(1
2
) d
g
, (2.26)
where
cf
F
and
ff
F
are the fracture stresses of a carbide ferrite and a ferrite
ferrite interfaces accordingly,
cf
and
ff
are the eective surface energies of
a carbide ferrite and a ferrite ferrite interfaces accordingly, d
c
and d
g
are
carbide and ferrite grain sizes accordingly, E is the elasticity modulus and is
the Poissons ratio.
Ritchie et al. (1973) postulated that the condition of equation (2.24) has
to be satised over a distance of two grain sizes ahead of the crack tip for the
fracture advance to take place. This is commonly called the critical distance
idea (Thompson and Knott, 1993).
Later Curry and Knott (1978) proposed a statistical analysis of eligible
particles that can be found within the critical distance. An eligible particle
is a cracked particle with the crack length equal or greater than the critical
one. Their conclusion was that a very small percentage of large particles have
a disproportionate inuence on the fracture resistance.
2.2.2 Weakest link models
Beremin (1983) developed the idea of eligible particles into a weakest link
statistical model. According to this model a certain volume, V , of material
ahead of the crack tip (usually the volume of the plastic zone) is assumed to
have a distribution of microcracks of dierent lengths. Catastrophic failure is
assumed to take place if a crack of critical length is found in this volume. This
microcrack is a weakest link, hence the name of the model. The probability of
failure is the probability of nding such microcrack.
It is further assumed that the volume V can be divided into smaller volumes
V
0
, which must be big enough so that the probability of nding a microcrack
of critical length is not negligible. At the same time V
0
must not be too big so
that one can assume that the stress state is homogeneous over V
0
. Thus usually
2.2. MICROANALYSIS OF BRITTLE FRACTURE 29
V
0
is chosen to include several grains.
The failure probability takes the following form (Beremin, 1983; Lin et al.,
1987; Ruggieri, 1998):
= 1 exp
_
_
V
0
1
V
0
_

0
g(S)dS
_
, (2.27)
where g(S)dS is the number of microcracks per V
0
with stresses required to
propagate them between S and S +dS.
Usually a three-parameter Weibull probability distribution function (Wei-
bull, 1951) is used to express g(S)dS:
_

0
g(S)dS =
_
th
u
_
m
, (2.28)
where
I
is a maximum principal stress in V
0
, m is a shape parameter,
u
is a
scale parameter and
th
is an oset parameter, a threshold stress, required to
propagate the largest feasible microcrack.
By substituting equation (2.28) into (2.27) one can obtain:
= 1 exp
_
u
_
m
_
, (2.29)
where
w
=
_
1
V
0
_
V
0
(
I
th
) dV
_
1/m
(2.30)
is called Weibull stress (Beremin, 1983).
A progressive brittle fracture statistical model based on chain-of-bundles
statistics (G ucer and Gurland, 1962) was proposed by Ruggieri et al. (1995).
In this model several critical events are allowed before the catastrophic failure
takes place. The analysis leads to Weibull statistics and eectively to the same
relations as expressed by equations (2.29) and (2.30) (Ruggieri, 1998).
Other forms of equation (2.28) can be used. Kroon and Faleskog (2002)
introduced the inuence of applied strain on g(S)dS and used an exponential
distribution instead of Weibull:
_

0
g(S)dS = c
p
eq
_
exp
_
I
_
2
_
exp
_
th
_
2
__
, (2.31)
where
m
and c are material parameters;
m
corresponds to the stress needed
to propagate a mean size microcrack.
The authors claim that their model predicts the inuence of constraint on the
failure probability better than the model based on the three-parameter Weibull
distribution.
However, as pointed out by Wallin (1991), even though the models may
dier considerably in their basic assumptions of the microscopic fracture mech-
anism, macroscopically most of them still yield an identical result.
2.2.3 Crack arrest
The problem of crack arrest received somewhat less attention than the issue of
crack nucleation and growth. Perhaps this is due to fact that in many appli-
cations crack arrest does not happen, i.e. a brittle crack would not stop until
the end of a specimen or a structure component is reached. This situation is
described perfectly well by the critical event analysis.
Generally a running brittle crack can arrest if the applied stresses decrease
with increasing crack length (e.g. if the test is performed under displacement
control) or if a crack hits an area of ne grains. According to equation (2.26)
the fracture stress of a ferrite ferrite interface is inversely related to the grain
size, so the ne grain region will represent a signicant obstacle to an advancing
crack. There is some experimental evidence in support of this idea (Malik et al.,
1996; Jang et al., 2003).
There is also some experimental evidence that a high-angle misorientation
boundary can act as a crack arrester or at least retard or inhibit the crack
propagation (Zikry and Kao, 1996). Nohava et al. (2002) reported crack arrest
in A508 Class 3 steels at grain boundaries with 55
60
misorientation angles.
2.3 Coupled ductile-brittle fracture modelling
The main two problems of coupled ductile-brittle fracture modelling were dis-
cussed briey in Chapter 1. These are high computational costs and conicting
demands for the nite element mesh size.
2.3. COUPLED DUCTILE-BRITTLE FRACTURE MODELLING 31
2.3.1 Size scales
Rousselier et al. (1989) proposed an inter-inclusion spacing, l
c
, as a ductile
fracture propagation step. Their metallographic examinations resulted in the
value l
c
= 0.55 mm for A508 steel.
Tvergaard and Needleman (1984) introduced D
0
, the initial spacing between
particle centres. They reported D
0
0.1 0.14 mm for an unspecied high
strength steel.
Some modern steels contain very few or indeed no detectable larger inclusions
(typically MnS). Some authors suggested a spacing between larger precipitates
as a suitable measure of a ductile fracture advance step. For a high purity
laboratory rolled thermo mechanically controlled rolled (TMCR) steels Davis
(2003) suggested the values of around 0.01 mm.
A ferritic grain size for tempered bainitic microstructures and a lath packet
size for microstructures related to segregated bands were linked to brittle frac-
ture propagation step by Beremin (1983). The values reported for A508 steel
were 0.011 mm for a ferritic grain size and 0.067 mm for a lath packet size.
These values led to the choice of the appropriate reference material volume, V
0
,
(section 2.2.2). V
0
was taken as a cube with side 0.05 mm, which for A508
includes about 8 grains.
The concept of a damage cell or a computational cell (Xia and Shih, 1996;
Faleskog and Shih, 1997) is used to introduce the above microstructurally sig-
nicant size scales into the local approach to fracture mechanics methods.
Two ways of implementing the damage cell concept via FE methods have
been explored over the years.
The easiest approach is to associate a damage cell with each FE in or near
the damage zone. This assumes constructing the mesh of the damage zone with
damage cell sized FEs (Tvergaard and Needleman, 1984; Rousselier et al., 1989;
Howard et al., 1996; Xia and Shih, 1996; Koppenhoefer and Dodds, 1998).
In the other method FE sizes are not xed to that of the damage cell.
Instead an additional size parameter is introduced into the model. This results
in a mesh-independent, non-local use of the theory (Bilby, Howard and Li, 1994;
Howard et al., 2000).
The damage cell sizes reported in the literature are 0.1 0.5 mm for ductile
damage and 0.005 0.05 mm for brittle fracture. Whatever the exact values
for the ductile and brittle cell sizes are, they are substantially dierent. It is
therefore dicult to accommodate both damage cell sizes within one FE mesh.
The compromise approach is to use a unied damage cell for both types of
fracture. A damage cell of 0.125 mm was used by Sherry et al. (1998); Burstow
(1998); Howard et al. (2000) as a reasonable compromise between 0.05 mm
brittle and 0.25 0.5 mm ductile damage cells. Howard et al. (2000) reported
that performance of such a compromise model is virtually indistinguishable
from that of a more complicated mesh-independent model (Bilby, Howard and
Li, 1994).
2.3.2 Brittle fracture as a postprocessing operation
This approach involves two stages.
At the rst stage a nite element solution is obtained using the local ap-
proach model for ductile damage. The stress history of all FEs in the plastic
zone is saved during the analysis.
The second stage consists of applying the weakest link statistical model to
the stress evolution data. The result is a probability of brittle failure as a
function of crack advance or time.
Various combinations of ductile and brittle models described in sections 2.1
and 2.2 can be used. The most popular are the GTN + Beremin (sections 2.1.4
and 2.2.2) (Xia and Shih, 1996; Xia and Cheng, 1997; Koppenhoefer and Dodds,
1998) and Rousselier + Beremin (section 2.1.6 and 2.2.2) (Eripret et al., 1996;
Sherry et al., 1998; Burstow, 1998; Howard et al., 2000).
This approach assumes a loss of stability associated with rapid loss of sti-
ness somewhere in the damage zone as the onset of the catastrophic brittle
fracture. So this model can only predict the time of the critical event. Neither
the cleavage initiation site nor the shape of the crack can be predicted.
2.3.3 Folch model
In the model introduced by Folch the onset of cleavage of each damage cell is
assessed individually (Folch, 1997; Folch and Burdekin, 1999). In other words
the integration in equation (2.30) is performed over a volume of material within
2.4. MODEL CALIBRATION 33
individual cell. It is easy to see that if the reference volume, V
0
, is equal to the
cell volume and the threshold stress,
th
, is zero then the Weibull stress,
w
, is
just the maximum principal stress. Therefore the equation 2.29 will have the
following form:
= 1 exp
_
u
_
m
_
, (2.32)
so that the probability of cleavage is based only on the ratio of the maximum
principal stress to the scale parameter of a Weibull distribution. Such a con-
dition is very similar to the criterion for the onset of cleavage expressed by
equation 2.24.
In this approach the probability of cleavage of each cell is calculated at the
same time as its constitutive response. So both ductile and brittle fractures can
be modelled simultaneously. What is more important is that the progressive
element to element brittle fracture propagation can be simulated. The cleavage
initiation sites can now be identied and the brittle crack front can be obtained
explicitly. The authors reported good agreement with the results of Charpy and
the fracture toughness tests (Folch, 1997; Folch and Burdekin, 1999).
However the model is still limited by the compromise cell size.
2.4 Model calibration
Any continuous ductile damage or a statistical cleavage model has to be cali-
brated for a particular material, so that model parameters can be considered
true material properties.
A three-stage calibration of a GTN model was proposed by Faleskog et al.
(1998). In the rst stage the parameters of the constitutive equation of the
model, q
1
, q
2
and q
3
, are tuned so that GTN model predicts the same void
growth as the model of a discrete spherical cavity. The second stage consists
of tuning the critical and nal fracture values of void volume fraction, f
c
and
f
f
, using a coalescence mechanics (Faleskog and Shih, 1997). The last stage is
the fracture process calibration in which a ductile cell size, L
D
, and an initial
void volume fraction, f
0
, are tuned by reproducing the behaviour of a fracture
test. This can be a fracture toughness, a single-edge-notched (SEN) bending or
a SEN tensile test (Gao et al., 1998).
The above analysis is equally applicable to any other continuous ductile
damage model.
In practice however many people use a simplied calibration procedures. In
many cases L
D
is calculated with an empirical relationship, e.g. L
D
= 2N
1/3
v
(Rousselier, 1987) or L
D
= 5N
1/3
v
(Rousselier et al., 1989), where N
v
is the
average number of inclusions per unit volume. Empirical relationships are also
used for f
0
. Of these the most popular is Franklins formula based on the Mn
and S contents in a steel (Franklin, 1969). However other estimates are also
used, e.g. f
0
=

6
d
x
d
y
d
z
N
v
, where d
x
, d
y
and d
z
are the average inclusion sizes
in three perpendicular dimensions (Batisse et al., 1987). The rest of the model
parameters are tuned to reproduce the results of a fracture test.
To tune the shape, m, and the scale,
u
, parameters of a Weibull distribution
for the weakest-link statistical model the maximum likelihood method is usually
used (Khalili and Kromp, 1991; Burstow, 1998; Gao et al., 1999).
2.5 Conclusion
Although considerable success in the prediction of failure of engineering struc-
tures has been achieved over the last twenty years with the use of the local
approach to fracture, there are several important problems that demand fur-
ther investigation.
Among signicant achievements of the coupled ductile-brittle fracture mo-
delling and the local approach to fracture in general one can list successful
predictions of all four spinning cylinder tests designed by AEA Technology, Ris-
ley, UK (Lidbury et al., 1994; Bilby, Howard, Othman, Lidbury and Sherry,
1994; Howard et al., 1996) and of the NESC (Network for Evaluating Steel
Components) experiment (Sherry et al., 1998).
However all attempts to date to transfer cleavage results from notched tensile
tests to precracked specimens failed (Howard et al., 2000).
The other two problems are long computational times and the conicting
demands of the FE mesh size.
The author believes that the last two problems are rooted in the fact that
the local approach to fracture utilises nite element methods as its vehicle. This
2.5. CONCLUSION 35
vehicle might not be fully appropriate for use in microstructure-related fracture
analysis.
The authors thesis is that a combination of cellular automata and nite
element methods is more suitable for this task. The next chapter gives the
presentation of the proposed approach.
Chapter 3
The CAFE solution
3.1 A CAFE model
A combination of cellular automata and nite elements (CAFE) has been used
successfully for solidication (Gandin et al., 1999; Vandyousse and Greer,
2002), static recrystallisation (Raabe and Becker, 2000) or oxide scale failure
(Das et al., 2001; Das, 2002; Das et al., 2003) modelling.
The CAFE model proposed here is a logical continuation of works by Bilby,
Howard and Li (1994), Folch (1997), Folch and Burdekin (1999), Raabe and
Becker (2000) and Das (2002). The structure of this model was rst presented
by Beynon et al. (2002).
As opposed to pure nite element fracture modelling, where a nite element
is a structural and material unit simultaneously, the present model, as indeed all
CAFE models, separates the structure from the material. Separate independent
entities are used to carry structural and material information.
A nite element is completely dened by its stiness matrix, D
ijkl
, and by
the interpolation functions, N
p
(
k
):
D
ijkl
=

ij
kl
(3.1)
u
ij
(
k
) = N
p
(
k
)u
p
ij
(3.2)
where
ij
is the Cauchy (or true) stress tensor,
kl
is the logarithmic (or true)
strain tensor, u
p
ij
is the displacement tensor at the nite element node p, p =
1 . . . P, P is the total number of nodes in a nite element and
k
is the nite
37
38 CHAPTER 3. THE CAFE SOLUTION
element parametric coordinates tensor, k = 1, 2, 3 and
k
= [1 . . . +1]; u
ij
(
k
)
is the displacement tensor at point
k
of a nite element (HKS, 2001).
If heat transfer is not taken into account then material behaviour at an
integration (or Gaussian or material) point r is described by a constitutive
equation of the form

r
ij
= f(
r
ij
) (3.3)
where
r
ij
and
r
ij
are stress and strain rate tensors respectively at an integration
point r, r = 1 . . . R, R is the total number of integration points per nite
element.
The separation of a conventional material nite element into a structural
and a material units is shown schematically in Figure 3.1.
Material FE Structural FE Material
= +
D
ijkl
, N
p
(
k
), D
ijkl
, N
p
(
k
)
r
ij
= f(
r
ij
)

r
ij
= f(
r
ij
)
Figure 3.1: A nite element as a structural and a material unit.
In a CAFE model the role of material unit is given to an appropriate number
of arrays of cellular automata (CA).
A CA is a discrete time entity composed of a nite number of cells. The
space of cell states is also discrete. In a classical CA formulation (Von Neumann,
1966) the state of each cell
m
at time t
i+1
is completely dened by the state
of this and the neighbouring cells at time t
i
:
m
(t
i+1
) = (
m
(t
i
),
m
l
(t
i
)) (3.4)
where
m
l
(t
i
) is the state of cell l from the neighbourhood of cell m at t
i
, l =
1 . . . L and L is the number of cells in the neighbourhood of cell m, m = 1 . . . M,
M is the total number of cells in the CA; is the transition rule and i N.
3.1. A CAFE MODEL 39
Thus a CA is completely dened by the initial state of each cell, by the
transition rules for each cell and by the neighbourhood of each cell. Usually the
same transition rules and neighbourhood are applied to all cells in a CA.
One important property of the CA dened above is that in itself it is a
non-spatial entity. This means that cells do not need to have any size, shape
or location in physical space for the successful functioning of a CA. Moreover
whatever the spatial meaning given to cells there will be no eect on the CA
functioning. This property makes CA a very general tool suitable for numerous
applications in mathematics and engineering.
A clear spatial meaning has to be given to cells in the present CAFE model
since the purpose of a CA is a representation of material behavior where size
scale is important. We shall relate a CA cell to a damage cell.
A damage cell or a computational cell concept was described in section
2.3.1. It was shown there that the microstructurally signicant size scales related
to ductile and brittle fractures are dierent. Therefore a material has to be
modelled with damage cells of two distinctly dierent sizes. This can be done
easily if a CA is chosen to represent material behaviour.
Two independent CAs, called hereafter the brittle CA array and the ductile
CA array, are created. The cell size in the brittle CA array is related to the
brittle damage cell size. Accordingly the ductile CA array cell size is related to
the ductile damage cell size.
The cubic shape of CA cells is adopted by analogy with the square (in
two dimensions) or cubic (in three dimensions) damage cells routinely used in
pure nite element fracture modelling (Xia and Shih, 1996; Howard et al., 2000).
Cubic CA cells are also the easiest for visualisation which is a problem for three-
dimensional structures. Finally the CA cell neighbourhood can be dened very
easily for a cubic cell.
A 26-cell neighbourhood is adopted in the present model for each cell. If one
imagines a 333 = 27 cell cube then the 26 cells lying around the central one
are its neighbourhood. Six cells of this neighbourhood have a common side with
the central cell; 12 cells have a common edge and 8 a common corner. Such
a neighbourhood is a three-dimensional analogy of Moores two-dimensional 8-
cell neighbourhood (Hesselbarth and G obel, 1991; Das, 2002). The properties
of this 26-cell neighbourhood are given in Appendix A.
A CA must have self-closing boundary conditions if each cell is to have the
same 26-cell neighbourhood. Self-closing means that for a cell lying at the edge
of a CA the corresponding cells of the opposite edge are considered adjacent.
So a 26-cell neighbourhood of an edge cell consists of cells located at opposite
CA edges.
The 26-cell neighbourhood of a corner cell is shown in Figure 3.2.
Figure 3.2: A 26-cell neighbourhood of a corner cell in a CA with self-closing
boundary conditions.
A corner cell, x, is located at a corner of a three-dimensional cubic CA.
This corner is an intersection of three CA edges. The numbers of the neighbo-
uring cells are shown according to the convention given in Appendix A. In the
projection shown in Figure 3.2 the neighbouring cells 9 and

9 are located ex-
actly behind cell x and are therefore not visible. Cell

9 is situated immediately
behind cell x, and cell 9 occupies the corner of CA opposite to cell x.
The CA structure described above is a classical CA formulation (Von Neu-
mann, 1966). We shall now depart from the classical CA model by assigning a
set of properties to each cell of a CA. These are the properties which are set at
the beginning of the simulation and remain constant throughout the analysis.
3.1. A CAFE MODEL 41
We shall also characterise a cell by a set of time-dependent state variables.
By adding the cell properties and state variables to the right part of equation
(3.4) we get the following evolution equation:
m
(t
i+1
) = (
m
(t
i
),
m
l
(t
i
),
n
m
,
n
l
,
q
m
(t
i
)) (3.5)
where
n
m
is property n of cell m, n = 1 . . . N, N is the total number of properties
of each cell;
n
l
is property n of a neighbouring cell l and
q
m
(t
i
) is a state variable
q of cell m at time t
i
, q = 1 . . . Q, Q is the total number of state variables dened
at each CA cell.
For simplicity we shall require that all cells of a CA have the same N and
the same Q. If this requirement is met then all cells of a CA can be processed
according to a unied algorithm.
We shall use the cell properties,
n
m
, to store some intrinsic material infor-
mation throughout the analysis. On the other hand the state variables,
q
m
(t
i
),
would come from the solution of material constitutive equations at time t
i
.
The number of cell properties and state variables is theoretically unlimited.
Exactly which material properties and solution-dependent state variables are
being represented in a CA depends on the particular realisation of the above
CAFE generalisation.
Finally we have to address the fact that two CA arrays, the brittle and the
ductile, occupy the same physical space. Therefore the state of cell m of one
CA array will depend on the states of a group of S corresponding cells of the
other CA array. This is necessary to ensure that any loss of material integrity,
whether due to the ductile or the brittle failure mechanism, is accounted for in
both CA arrays.
The full transfer rules for both CA arrays thus have the following form:
m(D)
(t
i+1
) =
D
_
m(D)
(t
i
),
m
l(D)
(t
i
),
n
m(D)
,
n
l(D)
,
q
m(D)
(t
i
),
s(B)
(t
i+1
)
_
(3.6)
m(B)
(t
i+1
) =
B
_
m(B)
(t
i
),
m
l(B)
(t
i
),
n
m(B)
,
n
l(B)
,
q
m(B)
(t
i
),
s(D)
(t
i+1
)
_
(3.7)
where subscripts D and B refer to cells from the ductile and the brittle CA
arrays accordingly, cell s belongs to the group of S cells of one CA array, the
states of which will have an inuence on cell m of the other CA array, s = 1 . . . S,
1 S M.
The number of cells S of one CA array which will aect the state of cell m
of the other CA array is dicult to establish exactly. As will be shown later S
depends on the total number of cells in each array, M
D
and M
B
.
The brittle and the ductile CA arrays are totally independent of each other
as far as their construction is concerned. This means that the number and
the types of cells states,
m(D)
and
m(B)
, the total numbers of cells in these
arrays, M
D
and M
B
, the total numbers and types of cell properties, N
D
and
N
B
, and state variables, Q
D
and Q
B
, and nally the transfer functions,
D
and
B
can be chosen for each array independently. This gives us great freedom as
to how exactly material behaviour is represented through the two CA arrays.
However the cell states in the ductile CA array are aected by the states
of cells in the brittle CA array and vice versa. This property ensures that any
change in material integrity, no matter what fracture mechanism caused it, is
accounted for in both CA arrays.
Similarly to the way we introduced time-dependent state variables at each
CA cell,
q
m
(t
i
), to link the state of each cell with the solution of material con-
stitutive equations, we shall now introduce solution-dependent state variables
Y
r
a
at each nite element integration point r linked to the states of both CA
arrays:
Y
r
a
(t
i+1
) =
_
m(D)
(t
i+1
),
m(B)
(t
i+1
)
_
(3.8)
where Y
r
a
(t
i+1
) is state variable a at time t
i+1
and integration point r, a =
1 . . . A, A is the total number of state variables per integration point and
is the CA to FE transfer function. The same transfer function is used for all
material points.
Up until now we described the general principles of the ductile and the brittle
CA organisation and the link between them. The link between the FE and the
CA parts of a CAFE model depends on the exact technical realisation of a
CAFE generalisation. The rest of this chapter will deal only with the CAFE
model realised in the present work.
3.2. THE FULL MODEL 43
3.2 The full model
The CAFE model was realised via the user material subroutine VUMAT in the
Abaqus/Explicit nite element code. This program utilises explicit dynamic
integration of the equations of motion. Reduced integration 8-node nite el-
ements C3D8R (HKS, 2001) were used to mesh the anticipated damage zone.
These elements have only one integration point (R=1).
The explicit dynamic version of the Abaqus code was chosen because of the
element removal feature which is not available in the Abaqus/Standard. The
removal of dead nite elements from the mesh is necessary in large deformation
analysis. Otherwise the dead elements, which have the highest strains, might
turn inside out. The solution will terminate in this case.
The general structure of a CAFE model is shown in Figure 3.3.
FE Ductile CA
ij
(t
i+1
),
ij
(t
i
)
ij
(t
i+1
), Y
a
(t
i+1
)

m(D)
(t
i+1
)
m(B)
(t
i+1
)
Brittle CA
Figure 3.3: Flow of information between the three parts of the full CAFE model.
The ow of information between the three parts of the CAFE model is shown
schematically with the arrows. However, the scheme does not show how the data
is being processed within each part of the whole model.
3.2.1 The ductile CA array
The Rousselier continuous ductile damage model is used as a material constitu-
tive routine for the ductile CA array (section 2.1.6).
As we said in section 3.1 the ductile CA cells are related to the ductile
damage cells. Therefore the total number of cells per ductile CA, M
D
, has to
be chosen so that the linear size of an individual ductile CA cell is close to
the ductile damage cell size, L
D
. If we assume a cubic nite element of size
L
FE
L
FE
L
FE
then the following equation can be used to choose M
D
:
L
FE
3
M
D
= L
D
(3.9)
where
3
M
D
is the number of cells per dimension of a cubic ductile CA.
Each ductile cell m can take one of the two possible states,
m(D)
,: alive or
dead. In the beginning of the simulation
m(D)
(t
0
) = alive.
Each ductile cell m has only one cell property (N
D
= 1), this is the value of
initial void volume fraction,
1
m(D)
= f
m
0
. A random number generator is used
to generate f
m
0
for each cell at the beginning of the simulation.
Each ductile cell m carries only one state variable (Q
D
= 1), this is the
current value of the damage parameter,
1
m(D)
(t
i
) =
m
(t
i
).
The same critical value of the damage variable,
F
, is used for all ductile
cells. Therefore
F
is a model parameter rather than a cell property.
The state of each cell m is determined by the following criterion:
m(D)
(t
i+1
) =
_
_
_
alive if
m
(t
i+1
) <
F
dead otherwise
(3.10)
3.2.2 The brittle CA array
Similarly to section 3.2.1 the link between the brittle damage cell size, L
B
, and
the total number of brittle CA cells is as follows:
L
FE
3
M
B
= L
B
(3.11)
where
3
M
B
is the number of cells per dimension of a cubic brittle CA.
As was shown in section 2.3.1 a brittle damage cell is typically 10 20 times
larger than the mean (or median or mode) grain size.
Ideally, a brittle cell size has to be related to a grain size for a progressive
grain-to-grain fracture simulation. However for steels with small grain sizes this
will lead to extremely high numbers of brittle cells.
The following compromise approach is proposed in this work. The brittle
cell size is chosen with equation (3.11). However a randomly generated grain
size value is assigned to each brittle cell. Therefore computational eciency
can be achieved while some real metallurgical data is retained in the brittle CA
array.
According to the present understanding of brittle fracture initiation (section
2.2.1) a brittle crack typically initiates from a hard particle. Most usually
this is a grain boundary carbide or a large inclusion, e.g. MnS. We simplify
this idea for the purpose of the present modelling and formulate the following
necessary condition for brittle crack initiation at a particular cell. Only cells
with an adjacent grain boundary carbide can initiate brittle fracture. Thus large
inclusions are not taken into account at present. However, as will be shown in
Chapter 5 the inuence of large inclusions can be easily incorporated into the
model if only the information regarding the size, number and locations of such
particles is available.
Accordingly a special state of the brittle CA cell is created alive with a grain
boundary carbide or simply aliveC. Only brittle cells with
m(B)
(t
0
) = aliveC
can initiate a brittle crack.
We have to address the problem of synchronising both CA arrays. All ductile
failures must be reected into the brittle CA array. However a distinction must
be made between the brittle cells failed due to the brittle failure mode, and
those which were made dead articially to synchronise the integrity of both CA
arrays. A special state of brittle CA cell is created dead in the ductile CA
array or simply deadD.
Finally the state deadB is reserved for the brittle CA cells which fail when
the brittle failure criterion is satised.
Thus each brittle cell can take one of the four possible states,
m(B)
,: alive,
aliveC, deadB or deadD.
In the beginning of the simulation
u(B)
(t
0
) = aliveC and
v(B)
(t
0
) =
alive, u = 1 . . . U, v = 1 . . . V , U + V = M
B
. So the fraction of brittle cells
which have a grain boundary carbide is = U/M
B
. A random number generator
is used to assign the initial state to cells.
Each brittle cell m has two cell properties (N
B
= 2), these are the fracture
stress,
1
m(B)
=
m
F
, and the grain orientation angle,
2
m(B)
=
m
. The grain
orientation angle is obtained with a random number generator.
The use of only one grain orientation angle is, of course, a modelling simpli-
cation. In principle two angles are required to describe a crystal orientation
(Kelly and Groves, 1970). However, what really matters in modelling crack
propagation from one grain to another, is the grain misorientation angle, that
is the minimum of all angles formed by the pairs of the crystallographic planes,
where each pair contains one crystallographic plane of one grain and one crystal-
lographic plane of the other grain. The calculation of the grain misorientation
angle is very easy if each grain is described by only one orientation angle. How-
ever, it is a much more computationally expensive task if the orientation of each
grain is described by two angles.
Perhaps it would be more correct to call
m
a grain orientation angle class or
type. This would imply that
m
denotes a particular combination of two grain
orientation angles. Accordingly if l is a grain (brittle cell) adjacent to grain m
then |
m

l
| is a dierence between the orientation classes (types) of grains
m and l. This is taken as an analogue of the true grain misorientation angle.
The fracture stress of a cell is linked to the size of the grain that this cell
embodies (equation (2.26) of section 2.2.1). A random number generator is used
to generate a grain size, d
g
, for each cell. Then a fracture stress is assigned to
each cell based on the generated grain size.
Each brittle cell m carries only one state variable (Q
B
= 1), this is the
current value of the maximum principal stress,
1
m(B)
(t
i
) =
m
I
.
As was shown in section 2.2.3 a high-angle misorientation grain boundary
can inhibit or even arrest crack growth (Nohava et al., 2002; Bhattacharjee
and Davis, 2002; Bhattacharjee et al., 2003). Again we simplify this idea and
formulate the following necessary condition for crack propagation. A crack
will propagate from one brittle CA cell, m, to another, l, if the misorientation
angle for these two cells, dened as the absolute value of the dierence of their
orientation angles, |
m
l
|, is less than a misorientation threshold,
F
. It is
assumed that
F
is a material property.
Finally the following simple propagation criterion is used which must be
satised in all cases. A brittle cell m will become dead at time t
i+1
if the
maximum principal stress,
m
I
(t
i
), is greater than or equal to the fracture stress
of this cell,
m
F
. This criterion is identical to that of Folch (1997), section 2.3.3,
if we require the probability of failure, = 1, and
u
=
F
.
Thus the state of each cell m is determined by the following criterion:
m(B)
(t
i+1
) =
=
_
_
deadB if
m
I
(t
i
)
m
F

__
m(B)
(t
i
) = aliveC
_

__
m
l(B)
(t
i
) = deadB
m
l(B)
(t
i
) = deadD
_

|
m
l
|<
F
__
m(B)
(t
i
) otherwise
(3.12)
3.2.3 The FE part
Each material point has three solution-dependent variables (A = 3): state,
Y
1
(t
i
), integrity, Y
2
(t
i
), and the fraction of brittle phase, Y
3
(t
i
). As follows from
equation (3.8) the FE state variables depend on the states of the brittle and
ductile CA cells.
First, the number of brittle CA cells with
m(B)
(t
i
) = deadB,
X
B
(B)
(t
i
) =
MB
m=1
m m :
m(B)
(t
i
) = deadB, (3.13)
the number of brittle CA cells with
m(B)
(t
i
) = deadD,
X
D
(B)
(t
i
) =
MB
m=1
m m :
m(B)
(t
i
) = deadD, (3.14)
the total number of dead brittle CA cells,
X
(B)
(t
i
) = X
B
(B)
(t
i
) +X
D
(B)
(t
i
), (3.15)
and the total number of dead ductile CA cells,
X
(D)
(t
i
) =
MD
m=1
m m :
m(D)
(t
i
) = dead, (3.16)
are calculated.
Then the FE state variables are calculated according to the following three
equations:
Y
3
(t
i
) =
X
B
(B)
(t
i
)
X
(B)
(t
i
)
(3.17)
Y
2
(t
i
) = 1
X
(D)
(t
i
)
X
max
(D)
X
B
(B)
(t
i
)
X
max
(B)
(3.18)
Y
1
(t
i
) =
_
_
_
dead if Y
2
(t
i
) 0
alive otherwise
(3.19)
where X
max
(D)
and X
max
(B)
are the maximum numbers of dead cells allowed in the
ductile and the brittle CA arrays respectively. If the number of dead cells in any
array exceeds its maximum then a crack (or ductile void linkage) is assumed to
propagate across the whole of the FE. The load-bearing capacity of this FE is
then considered zero and the FE is removed from the mesh.
So Y
1
[alive, dead]; Y
2
[1 . . . 1] and Y
3
[0 . . . 1]. In the beginning of
the analysis Y
1
(t
0
) = alive, Y
2
(t
0
) = 1 and Y
3
(t
0
) = 0 for all nite elements
included in the CAFE model. When a FE fails Y
1
(t
f
) = dead and Y
2
(t
f
) = 0,
where t
f
is the time of a FE failure.
3.2.4 How the model works
As seen from sections 3.2.1 and 3.2.2 there is a fundamental dierence between
the roles of the ductile and the brittle CA arrays in the full model. While the
ductile CA array is used to calculate material constitutive response, the brittle
CA array is only used to assess the onset of brittle fracture at each cell. Before
the brittle CA cells can be processed to decide if any of them have died in this
time increment, the material response has to be calculated via the ductile array
(Figure 3.3).
Below is the sequence of operations performed at each time increment for
each FE and the corresponding two CA arrays of the CAFE model.
All tensors given to the VUMAT subroutine are in the local material orienta-
tion.
Step 1. Strain increment tensor at time t
i+1
,
ij
(t
i+1
), and the stress tensor
at time t
i
,
ij
(t
i
), at the FE integration point are given by the Abaqus solver
to the VUMAT subroutine.
Step 2. The maximum principal stress,
I
(t
i
), and its direction cosines, d
k
(t
i
),
are calculated from
ij
(t
i
).
Step 3. The strain increment tensor at each ductile CA cell m at time t
i+1
,
m
ij
(t
i+1
), is calculated. The following criteria are used:
m
m
ij(D)
(t
i+1
) =
ij
(t
i+1
) (3.20)
m :
m(D)
(t
i
) = dead:
if d
l
k
d
k
(t
i
) 1 then
l
ij(D)
(t
i+1
) = c
D

ij
(t
i+1
) (3.21)
where c
D
> 1 is the concentration factor for the ductile CA array and d
l
k
are
the direction cosines of the line connecting the centres of cells m and l (see
Appendix A).
Then the strains at all M
D
ductile CA cells are scaled so that the average
of the cell strains gives the FE strain:
1
M
D
MD
m=1
m
ij(D)
(t
i+1
) =
ij
(t
i+1
) (3.22)
The condition of equation (3.21) means that if there is a dead ductile cell
then all neighbouring cells which lie on or near the plane perpendicular to the
direction of the maximum principal FE stress will receive some strain concen-
tration. This condition reects the strain concentration in material surrounding
a void. The sign rather than = is used in the if part of equation (3.21)
because there are only 13 pairs of neighbouring cells with unique combinations
of d
l
k
(see Appendix A). So it is very unlikely that d
l
k
d
k
(t
i
) = 1.
Step 4. The stress at time t
i+1
at each ductile cell m,
m
ij
(t
i+1
), and the value
of the damage variable,
m
(t
i+1
), are obtained via the solution of the set of
equations of the Rousselier continuous ductile damage model (section 2.1.6 and
Appendix B).
Step 5. The state of each ductile cell m at time t
i+1
,
m(D)
(t
i+1
), is obtained
according to equation (3.10).
Step 6. All dead ductile CA cells are reected into the brittle CA array (section
3.2.2). A special mapping function, M
DB
, distributes the array of ductile
CA cell states,
m(D)
(t
i+1
), across the brittle CA array. The result is the
synchronisation array of the brittle CA cell states,
m(BD)
(t
i+1
).
m(BD)
(t
i+1
) = M
DB
_
m(D)
(t
i+1
)
_
(3.23)
The subscript BD means that each brittle cell m has the state of the ductile
cell occupying the same physical space. The subscript m instead of usual m is
used in the right part of equation (3.23) because M
B
= M
D
, so m = 1 . . . M
D
and m = 1 . . . M
B
. This change of notation is only used in equations (3.23) and
(3.29).
The space of states of
m(BD)
is the same as of
m(D)
, either dead or alive.
The
m(BD)
(t
i+1
) = dead means that there is a ductile void in the physical
space associated with brittle cell m. So the state of brittle cell m is changed to
deadD to acknowledge this fact. This is expressed by the following equation
m(B)
(t
i+1
) =
_
_
deadD if
m(BD)
(t
i+1
) = dead
_
m(B)
(t
i
) = alive
m(B)
(t
i
) = aliveC
_
m(B)
(t
i
) otherwise
(3.24)
The work of mapping function M
DB
is illustrated in Figures 3.4.a and
3.4.b.
1 2 3 4 5
1
2
3
4
5
1 2 3 4 5 6 7 8 9 10111213
1
2
3
4
5
6
7
8
9
10
11
12
13
1 2 3 4 5
1
2
3
4
5
m(D)
(t
i+1
)
m(BD)
(t
i+1
)
m(B)
(t
i+1
)
m(DB)
(t
i+1
)
a. Ductile (in) b. Brittle (out/in) c. Ductile (out)
Figure 3.4: Illustration of the mapping operations, M
DB
and M
BD
.
Figure 3.4.a shows a two-dimensional slice of
m(D)
(t
i+1
). Dead cells are
grey and the white ones are alive. Figure 3.4.b shows a two-dimensional slice
of
m(BD)
(t
i+1
). It is easy to see that the locations of groups of grey cells
in
m(BD)
(t
i+1
) are close to locations of grey cells in
m(D)
(t
i+1
). Because
of the discrete nature of the CA space, the dead cells in
m(BD)
(t
i+1
) will
occupy exactly the same physical space as the dead cells in
m(D)
(t
i+1
) only if
the number of cells per linear brittle CA dimension,
3
M
B
, is divisible by the
number of cells per linear ductile CA dimension,
3
M
D
:
mod
3
_
M
B
M
D
= 0 (3.25)
In the example shown in Figures 3.4.a and 3.4.b
3
M
B
= 13 and
3
M
D
= 5,
so the condition of equation (3.25) does not hold. Therefore the locations of
dead (grey) cells in physical space in
m(D)
(t
i+1
) and in
m(BD)
(t
i+1
) are only
close to each other, but not identical.
Step 7. This step is the brittle CA analogue for the ductile CA (Step 3).
The maximum principal stress in each brittle CA cell m at time t
i+1
,
m
I
(t
i
),
is calculated. The following criteria are used:
m
m
I
(t
i
) =
I
(t
i
) (3.26)
m :
m(B)
(t
i
) = deadB : if d
l
k
d
k
(t
i
) 1 then
l
I
(t
i
) = c
B

I
(t
i
) (3.27)
m :
m(B)
(t
i
) = deadD : if d
l
k
d
k
(t
i
) 1 then
l
I
(t
i
) = c
D

I
(t
i
) (3.28)
where c
B
> 1 is the concentration factor for the brittle CA array.
The meaning of equations (3.27) and (3.28) for the brittle CA array is similar
to that of equation (3.21) for the ductile CA array (Step 3).
The state of each brittle cell m at time t
i+1
,
m(B)
(t
i+1
), is obtained ac-
cording to equation (3.12).
All dead brittle CA cells are reected into the ductile CA array. A spe-
cial mapping function, M
BD
, distributes the array of brittle CA cell states,
m(B)
(t
i+1
), across the ductile CA array. The result is the synchronisation
array of the ductile CA cell states,
m(DB)
(t
i+1
).
m(DB)
(t
i+1
) = M
BD
_
m(B)
(t
i+1
)
_
(3.29)
The subscript DB means that each ductile cell m has the state of the brittle
cell occupying the same physical space. The subscript m instead of usual m is
used in the left part of equation (3.29) because M
B
= M
D
, so m = 1 . . . M
D
and m = 1 . . . M
B
. This change of notation is used only in equations (3.29) and
(3.23).
The space of states of
m(DB)
is the same as of
m(D)
, either dead or alive.
The
m(DB)
(t
i+1
) = dead means that there is a brittle crack in the physical
space associated with ductile cell m. So the state of ductile cell m is changed
to dead to acknowledge this fact. This is expressed by the following equation
m(D)
(t
i+1
) =
_
_
_
dead if
m(DB)
(t
i+1
) = dead
m(D)
(t
i
) = alive
m(D)
(t
i
) otherwise
(3.30)
In contrast with the brittle CA array, no special cell state is created in the
ductile CA to distinguish between the dead ductile cells due to the ductile failure
mode and those made dead articially for synchronisation (equations (3.24) and
(3.30)). So the percentage of brittle phase can only be calculated via the brittle
CA.
The operation of mapping function M
BD
is illustrated in Figures 3.4.b and
3.4.c.
It is important to note that although each of the two mapping functions,
M
DB
and M
BD
, reects the state of one CA array onto the state of another
CA array only approximately, mapping errors do not accumulate. As shown
in Figure 3.4 the sequential application of both mapping operations produces
a cell state array identical to the initial one. It is easy to see that Figures
3.4.a and 3.4.c are identical. This property of mapping functions can be written
symbolically as follows
m(D)
(t
i+1
) = M
BD
_
M
DB
_
m(D)
(t
i+1
)
__
(3.31)
Step 10. All dead ductile cells receive zero stress:
m :
m(D)
(t
i+1
) = dead
m
ij
(t
i+1
) = 0. (3.32)
Step 11. The FE stress at time t
i+1
,
ij
(t
i+1
) is calculated as
ij
(t
i+1
) =
1
M
D
m
ij
(t
i+1
). (3.33)
So the stress at the FE integration point is the average of the stresses of all
ductile CA cells, including the dead cells.
Step 12. The FE solution-dependent variables at time t
i+1
, Y
k
(t
i+1
), are cal-
culated according to equations (3.17), (3.18) and (3.19).
Step 13. The FE stress tensor and solution-dependent variables at time t
i+1
,
ij
(t
i+1
) and Y
k
(t
i+1
), are returned by the VUMAT subroutine to the Abaqus
solver.
Step 13 completes the cycle.
3.2.5 Problems
1. High computational costs.
The running times of the full CAFE model are two three orders of magni-
tude smaller comparing with the pure nite element model where each ductile
(or brittle) CA cell is substituted by a nite element of equal size. From this
point of view the CAFE model is very fast.
However, a typical simulation time for a full 3D Charpy impact test is days
rather than hours (see section 4.1.3). This is because the integration of the
Rousselier damage model (Step 4, page 49) has to be performed for each ductile
cell. The system of two nonlinear equations is solved in this step using Newtons
iterative method (Appendix B). This is a relatively time-consuming operation.
2. Loss of precision due to averaging.
As is shown in Appendix B the elastic strain tensor,
e
ij
, which is used in
equation (B.8) and is updated according to equation (B.56), has to be stored
from one time increment to another throughout the analysis. The elastic strains
are typically very small for steels, therefore the components of
e
ij
must be cal-
culated with high precision. A small change in
e
ij
will lead to a very signicant
change of the Rousselier model response (
ij
and ).
However, the maintenance of high precision of
e
ij
is dicult due to a simplis-
tic strain redistribution criterion (Step 3, page 49) and due to stress averaging
(Step 11, page 52).
Each ductile cell will have a unique strain history due to a unique, randomly
assigned, initial void volume fraction, f
m
0
. A ductile cell m should theoretically
receive the strain increment at time t
i+1
,
m
ij(D)
(t
i+1
), based on the stress
at this cell in the previous time increment, t
i
,
m
ij
(t
i
). However Steps 11 of
increment t
i
and Step 3 of the next time increment, t
i+1
, eectively result in
the value of
m
ij(D)
(t
i+1
) being based on the averaged (macro or nite element)
stress. This strain increment might be quite far from the true strain increment
for cell m.
m
ij(D)
(t
i+1
) directly aects
e
ij
for cell m, as expressed by equation
(B.56). As a result
e
ij
might not be calculated accurately.
The greater the dissimilarity between the deformation histories of two neigh-
bouring ductile cell (related to a dissimilarity of f
0
for these cells), the further
e
ij
will be from the true value. In its extreme this loss of precision results in
e
ij
values which are so far from the accurate ones that the solution of the system
of equations (B.1) (B.7) is meaningless (e.g. negative
eq
) or it cannot be
found at all. The analysis will terminate at this stage.
The probability of encountering this problem obviously increases with the
number of nite elements in the model. For example this problem was never
encountered during the CAFE simulation of a single-FE model (sections 4.1.1
and 4.1.2). However it happened on several occasions during the Charpy test
modelling, where the CAFE model included 900 nite elements (section 4.1.3).
A simplied CAFE model is proposed which aims to solve the above prob-
lems.
3.3 The simplied model
The general structure of the simplied CAFE model is shown in Figure 3.5.
The major dierence between the full and the simplied models is that the
Rousselier model integration in the latter is performed at the nite element level.
The damage variable, (t
i+1
), is given to the ductile CA array in the simplied
model instead of the strain increment tensor,
ij
(t
i+1
). Accordingly only the
solution-dependent state variables are returned to the FE from the ductile CA
array, as the new macro stress tensor,
ij
(t
i+1
), is calculated already at the FE
3.3. THE SIMPLIFIED MODEL 55
FE Ductile CA
(t
i+1
),
ij
(t
i
)
Y
a
(t
i+1
)

m(D)
(t
i+1
)
m(B)
(t
i+1
)
Brittle CA
Figure 3.5: Flow of information between the three parts of the simplied CAFE
model.
level (Figures 3.3 and 3.5). All other dierences between the two models are
the consequences of this major one.
The critical value of the damage variable is now a randomly assigned cell
property,
1
m(D)
=
m
F
. So
m
F
instead of
F
is used in equation (3.10) for the
simplied model. The other properties of the ductile CA array are as described
in section 3.2.1.
The brittle CA array and the FE part are used exactly as in the full model
(sections 3.2.2 and 3.2.3).
3.3.1 How the model works
Below is the sequence of operations performed at each time increment for each
FE and the corresponding two CA arrays of the CAFE model.
Step 1. The same as Step 4 of the full model but performed at the FE level.
The damage variable at time t
i
, (t
i+1
), at the FE integration point is given by
the Abaqus solver to the VUMAT subroutine.
Step 2. The same as for the full model.
Step 3. The same as for the full model but applied to the damage variable.
The damage variable at each ductile CA cell m at time t
i+1
,
m
(t
i+1
), is
calculated. The following criteria are used:
m
m
(t
i+1
) = (t
i+1
) (3.34)
m :
m(D)
(t
i
) = dead:
if d
l
k
d
k
(t
i
) 1 then
l
(t
i+1
) = c
D
(t
i+1
) (3.35)
No scaling is performed as opposed to the full model.
Step 4. Not present in the simplied model.
Step 5. Step 9. The same as for the full model.
Step 10. Step 11. Not present in the simplied model.
Step 12. The same as for the full model.
Step 13. The same as for the full model but only solution-dependent variables
at time t
i+1
, Y
a
(t
i+1
), are returned by the VUMAT subroutine to the Abaqus
solver.
Step 13 completes the cycle.
It is easy to see that in the simplied model the roles of the ductile and the
brittle CA arrays are very similar as opposed to the full model (section 3.2.4).
In the full model the ductile CA array is used for the calculation of material
constitutive response and for the simulation of the fracture propagation while
the brittle CA array is only used to simulate fracture propagation. In the
simplied model both the brittle and the ductile CA arrays are used only for
the simulation of the fracture propagation at each CA scale. The material
constitutive response is calculated at the FE level.
A signicant reduction of computational time is thus achieved. Also there
is no averaging in the simplied model since Step 11 is not present. So the
accuracy of
e
ij
is not reduced.
3.4. IMPORTANT PROPERTIES 57
3.4 Important properties
Both the full and the simplied CAFE models have the following two important
properties.
1. In this particular realisation of a CAFE generalisation nite elements with
a single integration point (C3D8R) have been used. As a consequence there are
no stress or strain gradients across nite elements. Therefore mapping proce-
dures such as triangulation (Das, 2002) or the Wigner-Seitz algorithm (Raabe
and Becker, 2000) need not be applied to decide which CA cells fall under the
inuence of a particular integration point.
This property makes the present CAFE models fast. However, care must
be taken to ensure that nite element sizes are appropriate for regions of high
strain gradients as the accuracy of a single integration point nite element is
inevitably not as good as that of a nite element with several integration points.
2. Neither rotation nor deformation of a nite element are transferred to the
corresponding CAs. This is a very important property of the present CAFE
models. As a consequence the fracture propagation path can only be visualised
in the initial nite element conguration. This was a conscious decision.
The present work aims to model and understand the behaviour of a macro-
scopic sample of material based on micromechanics of fracture. If, however,
the emphasis is shifted towards the modelling and understanding of the frac-
ture at the microscale, then all CA arrays have to be rotated and deformed
using the deformation gradient tensor at the corresponding integration point
(Das, 2002). In this case the exact locations of all cells in physical space will
be known throughout the analysis and the fracture propagation path can be
visualised on a deformed mesh.
3.5 The list of model parameters
To complete this chapter, the full list of the parameters for both CAFE models
is shown below.
1. M
D
the total number of cells in the ductile CA.
2. M
B
the total number of cells in the brittle CA.
3. a probability density function, f(f
0
), for the full model and f(
F
), for the
simplied model.
4. a fraction of the brittle CA cells to have a grain boundary carbide.
5. a probability density function f(d
g
).
6. a probability density function f().
7. X
max
(D)
the maximum number of dead cells allowed per ductile CA.
8. X
max
(B)
the maximum number of dead cells allowed per brittle CA.
9. c
D
the concentration factor for the ductile CA array.
10. c
B
the concentration factor for the brittle CA array.
Chapter 4
Results
All results in this Chapter were obtained with the following soft- and hard-
ware. The VUMAT subroutine was written in FORTRAN 95, and compiled and
linked with Compaq Visual Fortran compiler, version 6.1 (Compaq, 1999). The
Abaqus/Explicit version 6.2-1 code was used. The platform used was a Pen-
tium III, 1 GHz PC run under Windows 2000 operating system. The amount of
operating memory used was approximately 42 MB for the examples in sections
4.1.1 and 4.1.2 and 50 MB for the examples in sections 4.1.3 and 4.2.1.
All examples in this Chapter illustrate the ability of the present CAFE mod-
els, both the full and the simplied, to simulate fracture propagation typical of
that in a thermomechanically control rolled (TMCR) microalloyed steel. Sec-
tions 4.1.1 and 4.1.2 illustrate some important features of the full CAFE model
by simulating the fracture propagation in a typical TMCR steel. Fracture of a
particular TMCR steel is modelled in sections 4.1.3 and 4.2.1.
4.1 The full CAFE model
4.1.1 Single FE, tension compression
This simple example is used to demonstrate the behaviour of the full CAFE
model.
The model consists of a single cubic 1mm1mm1mm nite element. The
initial shape of the nite element (dashed lines) and the boundary conditions
59
60 CHAPTER 4. RESULTS
(arrows) are shown in Figure 4.1.a.
1
2
3
0 0.002 0.004 0.006 0.008 0.01
0.4
0.2
0
0.2
0.4
0.6
Time, sec
D
i
s
p
l
a
c
e
m
e
n
t
,

m
m
a. Initial and deformed meshes b. Applied displacement, u
2
(t)
Figure 4.1: A single nite element model under alternating uniaxial tension and
compression.
The vertical displacement of the four bottom nodes is zero and for the four
top nodes it is described by the function u
2
(t) shown in Figure 4.1.b.
The material plastic properties are described by the power hardening law of
the form:
Y
=
Y0
_
p
eq
3G
Y0
+ 1
_
n
(4.1)
where
Y0
is the rst yield stress and n is the hardening exponent.
The simulation was performed at T = 20
C. For this temperature the values
Y0
= 447 MPa and n = 0.0575 were chosen based on data shown in Figure
4.11. The Youngs modulus and the Poisson ratio are E = 2 10
5
MPa and
= 0.3 respectively. The hardening curve is shown in Figure 4.2.
A 555 cell ductile (M
D
= 125) and a 202020 cell brittle (M
B
= 8000)
CA array were created. Thus the ductile damage cell size is L
D
= 1/5 = 0.2 mm
and the brittle damage cell is L
B
= 1/20 = 0.05 mm.
A two-parameter Weibull distribution was used to simulate the distribution
of the initial void volume fraction, f
0
. The shape parameter was taken as
W
= 2 and the scale parameter was taken as W
= 2.82 10
4
. The resulting
histogram of f
0
is shown in Figure 4.3.
The fracture stress value,
F
, was assigned to each brittle CA cell based on
the normal distribution with
F
= 1.910
3
MPa and STD(
F
) = 3.610
2
Mpa.
The resulting histogram of
F
4.1. THE FULL CAFE MODEL 61
0 0.5 1 1.5 2 2.5 3
400
450
500
550
600
650
700
Equivalent plastic strain

Y
,

M
P
a
Figure 4.2: Yield stress,
Y
(
p
eq
).
0 2 4 6 8
x 10
4
0
1
2
3
4
5
6
7
Mean =2.43e004
Median =2.31e004
Mode =2.33e004
STD =1.27e004
Max =6.16e004
Min =2.69e005
f
0
N
u
m
b
e
r

o
f

c
e
l
l
s
Figure 4.3: The initial void volume fraction histogram.
As this example does not simulate the behaviour of any particular steel, but
rather illustrates the CAFE model behaviour, the actual numerical values of f
0
and
F
assigned to each ductile or brittle cell are not very important, providing
they lie within the reasonable ranges of these parameters for TMCR steels.
Davis (2003) reported that no large second phase particles were observed in the
fracture surfaces of the broken tensile and Charpy samples of a TMCR material
500 1000 1500 2000 2500 3000 3500
0
10
20
30
40
50
60
70
Mean =1.90e+003
Median =1.91e+003
Mode =1.92e+003
STD =3.59e+002
Max =3.19e+003
Min =6.19e+002

F
, MPa
N
u
m
b
e
r

o
f

c
e
l
l
s
Figure 4.4: The fracture stress histogram.
described in section 4.1.3. However, if f
0
= 0, then the Rousselier model cannot
simulate material softening, equations (B.1) (B.7). Even if f
0
= 0 but very
small, then the softening behaviour predicted by the model is too slow, and the
strains to fracture are unrealistically high. Therefore the author had to choose
a reasonable value of

f
0
based on his previous experience, and that of other
people, of modelling ductile fracture (Shterenlikht et al., 2003; Andrews et al.,
2002; Burstow, 1998)
To further simplify this example no grain size, d
g
, or grain orientation, ,
distribution is used. Also for reasons of simplicity = 1, so that brittle fracture
can initiate at any brittle CA cell.
The values for other model parameters used in this example (section 3.5)
are c
D
= 1.4, c
B
= 1.6, X
max
(D)
= M
2/3
D
and X
max
(B)
= M
2/3
B
.
There is little guidance as to how to choose the values for the above four
parameters. We can assume that c
B
> c
D
as the strain concentration ahead of
a crack is higher than at the surface of a void. On the other hand M
2/3
is the
number of cells in a square section of a cubic CA array. So X
max
= M
2/3
is
chosen on the assumption that a FE will fail when a planar crack, perpendicular
to one of the three basis directions, crosses either of the two CA arrays.
The parameters of the Rousselier damage model (section 2.1.6) are: D = 3,
1
= 500 MPa and
F
= 8. These parameters are tuned on the CAFE modelling
of the Charpy test at the upper shelf (section 4.1.3).
The computation of this example took less than 10 minutes.
The modelling results are shown in Figures 4.5 4.9. Of these Figures 4.5
4.7 present results on the macro (FE) scale. The micro (CA) scale outcome of
the simulation is shown in Figures 4.8 and 4.9.
Figure 4.5.a shows the vertical stress against time,
22
(t). It is easy to see
that the failure of the nite element occurred before the second change of the
direction of applied displacement (Figure 4.1.b). Accordingly the shape of the
nite element at the end of the simulation is that of a compressed element, as
shown in Figure 4.1.a (solid lines).
0 0.002 0.004 0.006 0.008 0.01
750
500
250
0
250
500
750
Time, sec

2
2
,

M
P
a
0 0.002 0.004 0.006 0.008 0.01
0
250
500
750
Time, sec

e
q
,

M
P
a
a. Axial stress,
22
(t) b. Equivalent stress,
eq
(t)
0 0.002 0.004 0.006 0.008 0.01
0
0.5
1
1.5
Time, sec
E
q
u
i
v
a
l
e
n
t

p
l
a
s
t
i
c

s
t
r
a
i
n
0 0.5 1 1.5
0
250
500
750

e
q
,

M
p
a
c. Equivalent plastic strain,
p
eq
(t) d. Equivalent stress,
eq
(
p
eq
)
compression.
The plot of von Mises equivalent stress against time,
eq
(t), is shown in
Figure 4.5.b. Two points worth noting are the time delay at t = 4 10
3
sec
associated with elastic unloading and loading in the opposite direction, and the
dierence in the maximum
eq
values in tension and compression. The maxi-
mum
eq
in tension and compression are 585 MPa and 604 MPa respectively.
This dierence is due to the Rousselier model depending on the sign of the
mean stress,
m
(equation (2.16), section 2.1.6). This property of the Rousse-
lier model is in a good agreement with various experimental observations that
the positive
m
values cause much greater damage than negative values (e.g.
Bridgman, 1952).
Elastic unloading and loading can be seen also in the plot of equivalent
plastic strain against time,
p
eq
(t), shown in Figure 4.5.c.
The equivalent plastic strain equivalent stress plot,
eq
(
p
eq
), is shown in
Figure 4.5.d. It is easy to see, by comparing this plot with that of
eq
(t) shown
in Figure 4.5.b, that damage indeed accumulates much slower (in terms of
p
eq
)
in compression than in tension. For
p
eq
= 0.6, from
p
eq
= 0.4 to
p
eq
1,
eq
remains virtually constant thus indicating that there is no loss of load bearing
capacity during this large strain increment.
The evolution of three nite element solution-dependent state variables,
Y
1
, Y
2
and Y
3
is shown in Figures 4.6 and 4.7.
0 0.5 1 1.5
0
1
F
E

s
t
a
t
e
0 0.5 1 1.5
0
0.2
0.4
0.6
0.8
1
F
E

i
n
t
e
g
r
i
t
y
a. FE state, Y
1
(
p
eq
) b. FE integrity, Y
2
(
p
eq
)
compression.
The state of the nite element, Y
1
, changed from alive (1) to dead (0) at
the moment when its integrity, Y
2
, reached zero (Figures 4.6.a and 4.6.b). Each
step in Figure 4.6.b is associated with the failure of one or more cells in either
of the CA arrays. These steps are also visible in Figures 4.5.a, 4.5.b and 4.5.d.
The evolution of the brittle phase per FE, Y
3
, is shown in Figure 4.7. It
is important to note that the brittle phase can decrease during the simulation
according to equation 3.17. This is what happened in the present example. The
brittle phase in the end of the analysis is only Y
3
= 6.25 10
4
or 0.06%. In
practise such fracture would be routinely called a 100% ductile.
0 0.2 0.4 0.6 0.8 1 1.2 1.4
0
0.002
0.004
0.006
0.008
0.01
0.012
0.014
0.016
B
r
i
t
t
l
e

p
h
a
s
e
Figure 4.7: Brittle phase per FE, Y
3
(
p
eq
).
The state of the brittle CA array at the end of the simulation is shown
in Figure 4.8. As was said in section 3.4 the visualisation of the fracture on
a CA scale in the present model is only possible in the initial (undeformed)
conguration. Therefore the bounding box in Figure 4.8 represents the initial
shape of the nite element.
All results on the CA scale in this and the following examples were generated
with the EnSight visualisation software (CEI, 2002).
All aliveC cells are transparent and only the deadB (black) and deadD
(grey) cells are shown. In fact there is only one deadB cell and all other dead
cells are deadD. All deadD cells are grouped in 444 = 64 cell cubes because
there are 64 brittle CA cells per each ductile cell (M
B
/M
D
= 8000/125 = 64).
There is a tendency for the dead ductile cells to form clusters in XZ planes,
which are the planes perpendicular to the direction of the maximum principal
stress. However this tendency is quite weak and can only be seem in the bottom
0.5
0.500006
0.5
0.750001
X-Axis
0.750007
1
Y-Axis
0.750001
1.25
1.00001
1.5
0.500006
1.25001
Z-Axis
1
0.750007
1.50001
0.5
0.5
Y-Axis
1.00001
0.750001
1.25
X-Axis
1
1.25001
0.750001
1.25
1.50001
1.5
1.5
1
Z-Axis
1.25
1.5
X
Y
Z
Figure 4.8: The state of the brittle CA array at the end of the simulation,
p
eq
= 1.38. Grey cells are
m(B)
= deadD and black cells are
m(B)
= deadB.
layer of the dead ductile cells. This suggests that the concentration factor for
the ductile CA array, c
D
= 1.4, used in this analysis is not high enough.
On the other hand there is experimental evidence that there might be many
voids in the vicinity of the main ductile crack which never coalesce with other
voids or with the main crack (Puttick, 1959). Taking these observations into
account the simulation result shown in Figure 4.8 seems reasonable.
The sequence of eight snapshots of the brittle CA array throughout the
simulation is shown in Figure 4.9.
The fracture process started at
p
eq
= 0.38 when one ductile CA cell became
dead. Mapping function M
DB
(Step 6, page 49) then translated the location
of this cell into the locations of the corresponding 64 brittle cells which became
deadD. These cells constitute the grey block in Figure 4.9.b. Figure 4.10 shows
another projection of the same state of the brittle CA array as Figure 4.9.b.
The locations of the deadD cells are more clear in Figure 4.10. There are
4 4 5 = 80 cells adjacent to the deadD cells which lie on the XZ planes.
These are the cells which satisfy the condition of equation (3.28), i.e. these
cells will be given the maximum principal stress higher than the value for the
0.5
0.500006
0.5 0.750001
X-Axis
0.750007
1
Y-Axis
0.750001
1.25
1.00001
1.5
0.500006
1.25001
Z-Axis
1
0.750007
1.50001
0.5 0.5
Y-Axis 1.00001
0.750001
1.25
X-Axis
1
1.25001
0.750001
1.25
1.50001
1.5
1.5
1
Z-Axis
1.25
1.5
X
Y
Z
0.5
0.500006
0.5 0.750001
X-Axis
0.750007
1
Y-Axis
0.750001
1.25
1.00001
1.5
0.500006
1.25001
Z-Axis
1
0.750007
1.50001
0.5 0.5
Y-Axis 1.00001
0.750001
1.25
X-Axis
1
1.25001
0.750001
1.25
1.50001
1.5
1.5
1
Z-Axis
1.25
1.5
X
Y
Z
a.
p
eq
= 0 b.
p
eq
= 0.38
0.5
0.500006
0.5 0.750001
X-Axis
0.750007
1
Y-Axis
0.750001
1.25
1.00001
1.5
0.500006
1.25001
Z-Axis
1
0.750007
1.50001
0.5 0.5
Y-Axis 1.00001
0.750001
1.25
X-Axis
1
1.25001
0.750001
1.25
1.50001
1.5
1.5
1
Z-Axis
1.25
1.5
X
Y
Z
0.5
0.500006
0.5 0.750001
X-Axis
0.750007
1
Y-Axis
0.750001
1.25
1.00001
1.5
0.500006
1.25001
Z-Axis
1
0.750007
1.50001
0.5 0.5
Y-Axis 1.00001
0.750001
1.25
X-Axis
1
1.25001
0.750001
1.25
1.50001
1.5
1.5
1
Z-Axis
1.25
1.5
X
Y
Z
c.
p
eq
= 0.66 d.
p
eq
= 0.95
0.5
0.500006
0.5 0.750001
X-Axis
0.750007
1
Y-Axis
0.750001
1.25
1.00001
1.5
0.500006
1.25001
Z-Axis
1
0.750007
1.50001
0.5 0.5
Y-Axis 1.00001
0.750001
1.25
X-Axis
1
1.25001
0.750001
1.25
1.50001
1.5
1.5
1
Z-Axis
1.25
1.5
X
Y
Z
0.5
0.500006
0.5 0.750001
X-Axis
0.750007
1
Y-Axis
0.750001
1.25
1.00001
1.5
0.500006
1.25001
Z-Axis
1
0.750007
1.50001
0.5 0.5
Y-Axis 1.00001
0.750001
1.25
X-Axis
1
1.25001
0.750001
1.25
1.50001
1.5
1.5
1
Z-Axis
1.25
1.5
X
Y
Z
e.
p
eq
= 1.01 f.
p
eq
= 1.15
0.5
0.500006
0.5 0.750001
X-Axis
0.750007
1
Y-Axis
0.750001
1.25
1.00001
1.5
0.500006
1.25001
Z-Axis
1
0.750007
1.50001
0.5 0.5
Y-Axis 1.00001
0.750001
1.25
X-Axis
1
1.25001
0.750001
1.25
1.50001
1.5
1.5
1
Z-Axis
1.25
1.5
X
Y
Z
0.5
0.500006
0.5 0.750001
X-Axis
0.750007
1
Y-Axis
0.750001
1.25
1.00001
1.5
0.500006
1.25001
Z-Axis
1
0.750007
1.50001
0.5 0.5
Y-Axis 1.00001
0.750001
1.25
X-Axis
1
1.25001
0.750001
1.25
1.50001
1.5
1.5
1
Z-Axis
1.25
1.5
X
Y
Z
g.
p
eq
= 1.25 h.
p
eq
= 1.38
Figure 4.9: Damage propagation across the FE.
nite element integration point,
l
I
(t
i
) = c
D

I
(t
i
). Because of the self-closing
boundary condition (Figure 3.2, page 40) the black cell in Figure 4.10 is one of
these cells.
X Z
Y
Y-Axis 1.00001
1.25001
0.750007
0.500006
1.50001
1.5
1.5 1.25
1.25
X-Axis
1
1
X-Axis
0.750001
0.750001 0.5
0.5 1.5
1.25 1.5
1.25
1.00001
0.500006
0.750001
1.50001
0.750007
1.25001
0.5
Z-Axis
0.5 0.750001 1
1
Z-Axis
Y-Axis
Figure 4.10: The state of the brittle CA array at
p
eq
= 0.38.
At
p
eq
= 0.38,
I
= 576 MPa, (in this example
I
=
22
, Figure 4.5.a)
therefore
l
I
= 1.4 576 = 806 MPa. The data from Figure 4.4 shows that
there only 5 brittle CA cells with
F
806 MPa. Accidentally one of these
ve happened to be the black cell in Figure 4.10. So this cell became deadB in
the same time increment as the rst ductile cell became dead. Accordingly the
brittle phase at this moment was 1/64 = 0.0156 (Figure 4.7).
The data for
F
(Figure 4.4) and
22
=
I
(Figure 4.5.a) help to explain
why there is only one deadB brittle CA cell at the end of the simulation.
The maximum
I
in the analysis is
I
= 604 (Figure 4.5.a). Therefore the
maximum
l
I
is
l
I
= 1.4 604 = 845 MPa. However Figure 4.4 shows that
there are only 7 brittle CA cells with
F
845 MPa. Thus the probability that
max(
m
I
)
F
for a brittle CA cell m is only 7/8000 = 8.75 10
4
.
The above discussion reveals the importance of the whole fracture stress
distribution rather than its single characteristics (mean, median or mode) in
transitional ductile brittle fractures. It shows that transitional behaviour is
only possible if there is a reasonable scatter of the fracture stress values. If there
is no scatter and all brittle CA cells receive the same fracture stress, then the
resulting fracture will be either 100% brittle or 100% ductile, because all brittle
CA cells will behave as one.
The importance of the fracture stress is one of the major issues of the present
work and it will be illustrated in greater detail in the next examples.
4.1.2 Single FE, forward tension simulation of scatter
In this example the capacity of the full CAFE model to simulate the scatter
at the transitional temperatures is illustrated. Similarly to the previous exam-
ple (section 4.1.1), material properties were chosen to qualitatively reproduce
fracture propagation in a TMCR steel.
The model consists of a single cubic 1mm 1mm 1mm nite element
(Figure 4.1, dashed lines) under uniaxial tension.
The power hardening law of equation (4.1) was used. The rst yield stress,
Y0
, and the hardening exponent, n, are temperature-dependent as shown in
Figure 4.11.
200 150 100 50 0 50
400
500
600
700
800
900
Temperature,
o
C

Y
,

M
p
a
200 150 100 50 0 50
0.01
0.02
0.03
0.04
0.05
0.06
Temperature,
o
C

n
a. First yield stress,
Y0
(T) b. Hardening exponent, n(T)
Figure 4.11: Temperature dependence of material properties. Based on the
experimental data provided by Davis (2003).
As in the previous example (section 4.1.1) E = 2 10
5
MPa and = 0.3.
A 555 cell ductile (M
D
= 125) and a 101010 cell brittle (M
B
= 1000)
CA array were created. Thus the ductile damage cell size is L
D
= 1/5 = 0.2 mm
and the brittle damage cell is L
B
= 1/10 = 0.1 mm.
As in the previous example a two-parameter Weibull distribution was used
to simulate the distribution of the initial void volume fraction, f
0
. The shape
parameter was taken as W
= 2 and the scale parameter was taken as W
=
2.82 10
4
. The resulting histogram of f
0
in each of the simulations of this
example was similar to that shown in Figure 4.3.
The fracture stress value,
F
, was assigned to each brittle CA cell based on
the normal distribution with
F
= 1.910
3
MPa and STD(
F
) = 510
2
MPa.
An example of the resulting histogram of
F
0 500 1000 1500 2000 2500 3000 3500
0
1
2
3
4
5
6
7
8
9
10
Mean =1.89e+003
Median =1.90e+003
Mode =2.03e+003
STD =4.77e+002
Max =3.38e+003
Min =3.69e+002

F
, MPa
N
u
m
b
e
r

o
f

c
e
l
l
s
Figure 4.12: The fracture stress histogram.
The values for other model parameters are = 1, c
D
= 1.4, c
B
= 1.7,
X
max
(D)
= M
2/3
D
and X
max
(B)
= M
2/3
B
. As in the previous example the parameters
of the Rousselier damage model are D = 3,
1
= 500 MPa and
F
= 8.
The simulations were performed at six temperatures, from T = 196
C to
T = 0
C, three runs at each temperature.

All simulations of this example took less than 4 minutes.
The modelling results are shown in Figures 4.13 4.17. Of these, Figures
4.13 4.16 present results on the macro (FE) scale. The micro (CA) scale
results of the simulations are shown in Figure 4.17.
Figure 4.13 shows the evolution of the equivalent plastic stress,
eq
(
p
eq
),
during the simulation for all 18 runs.
The scatter of
p
eq
to failure is very low at 196
C, higher at 175
C, the
highest at 150
C, lower at 100
C and very low again at 50
C and at
0
C. This Figure 4.13 demonstrates that the full CAFE model can generate the
scatter in terms of
p
eq
to failure. And the level of this scatter is temperature-
dependent.
Figures 4.14 and 4.15 show the evolution curves of two of the FE solution-
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
0
100
200
300
400
500
600
700
800
900

e
q
,

M
P
a
0
o
C
0
o
C
0
o
C
50
o
C
50
o
C
50
o
C
100
o
C
100
o
C
100
o
C
150
o
C
150
o
C
150
o
C
175
o
C
175
o
C
175
o
C
196
o
C
196
o
C
196
o
C
Figure 4.13: Equivalent von Mises stress,
eq
(
p
eq
).
dependent variables: integrity, Y
2
(
p
eq
), and the brittle phase, Y
3
(
p
eq
).
It is easy to see in Figure 4.14 that the diversity of the shapes of the three
Y
2
(
p
eq
) curves at each temperature is maximal at 150
C and decreasing to-

wards 196
C and towards 0
C. The three Y
2
(
p
eq
) curves at 196
C (dark blue)
lie virtually on top of each other. However the three curves at 0
C (black) can
be easily distinguished one from another. Thus the scatter in terms of Y
2
(
p
eq
)
is also temperature-dependent.
Figure 4.15 shows the evolution of the fraction of the brittle phase per FE
during the simulation, Y
3
(
p
eq
). The curves at 196
C and at 175
C (dark and
light blue) are not visible because in all six simulations Y
3
= 1 from the onset
of the fracture propagation until the failure of the FE (compare with Figures
4.13 and 4.14).
The scatter in terms of the dierence in curve shapes at each temperature is
temperature-dependent. The six curves, three at 196
C and three at 175
C,
lie on top of each other. However for all other temperatures the shapes the
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
I
n
t
e
g
r
i
t
y
0
o
C
0
o
C
0
o
C
50
o
C
50
o
C
50
o
C
100
o
C
100
o
C
100
o
C
150
o
C
150
o
C
150
o
C
175
o
C
175
o
C
175
o
C
196
o
C
196
o
C
196
o
C
Figure 4.14: FE integrity, Y
2
(
p
eq
).
curves are distinctly dierent.
The values of the brittle phase at the point of a FE failure are not very clear
from Figure 4.15. They are shown separately in Figure 4.16 which is a typical
illustration of the transitional ductile brittle fracture behaviour.
In this example the lower shelf temperatures are below 175
C, the upper
shelf is above 50
C and the transitional temperatures are from 150
C to
100
C. The level of scatter in the transitional region is higher than in the

lower or in the upper shelf. It is interesting to note that the transition from the
lower shelf is quite sharp, whereas the transition from the upper shelf is much
smoother.
As in the example of section 4.1.1 the results of this section are rather qual-
itative than quantitative. Although the model is very capable of simulating the
transitional behaviour, the transitional temperature range might be shifted from
that obtained experimentally. This is primarily due to diculties of nding the
proper values for c
D
, c
B
, X
max
(D)
, X
max
(B)
, W
and W
(for the distribution of f

0
)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
B
r
i
t
t
l
e

p
h
a
s
e
0
o
C
0
o
C
0
o
C
50
o
C
50
o
C
50
o
C
100
o
C
100
o
C
100
o
C
150
o
C
150
o
C
150
o
C
175
o
C
175
o
C
175
o
C
196
o
C
196
o
C
196
o
C
Figure 4.15: Brittle phase per FE, Y
3
(
p
eq
).
200 150 100 50 0
0
0.2
0.4
0.6
0.8
1
Temperature,
o
C
B
r
i
t
t
l
e

p
h
a
s
e
Figure 4.16: The fraction of the brittle phase at the end of the simulation.
and STD(
F
). This issue is discussed in greater detail in section 4.2.1. On the
other hand there is no necking in a single FE model. Therefore the stress triax-
iality in this and the previous examples is very low, unlike in any real fracture
test.
Figure 4.17 shows the state of the brittle CA array at the point of FE
failure at six temperatures. The results of two of the three simulations at each
temperature are shown. As in Figures 4.8 and 4.9 the black cells are deadB and
the grey ones are deadD. In this example each ductile cell occupies the space
of 2 2 2 = 8 brittle cells since M
B
/M
D
= 1000/125 = 8.
A major brittle fracture plane is obvious at 196
C (Figures 4.17.a and

4.17.a) and at 175
C (Figures 4.17.c and 4.17.d). However in all four cases

there are several deadB cells apart from the main crack. These cells represent
1.50001
0.5 0.5
Y-Axis
1.25001
1.00001
0.750001
0.750005
X-Axis
0.500005
0.750001
1
0.5 0.5
1.25
0.750001
Z-Axis
1
X-Axis
1.5
0.750001
1
1.50001
1.25001
1.25
Y-Axis 1.00001
1.25
0.750005
1
Z-Axis
1.5
0.500005
1.5
1.25
Y
X
1.5
Z
1.50001
0.5 0.5
1.25001
Y-Axis
0.750001
1.00001
X-Axis
0.750005
1
0.500005
0.750001
0.5
1.25
0.5 0.750001
1.5
1.50001
X-Axis
1
1.25001
Z-Axis
Y-Axis
1
1.00001
0.750001 1.25
0.750005
1.5
0.500005
1.25
1
Z-Axis
Y
1.5
1.25
X Z
1.5
a. 196
C b. 196
C
1.50001
0.5
Y-Axis
0.5
1.25001
0.750001
X-Axis
1.00001
1
0.750005
0.500005
1.25
0.5 0.5
1.5
0.750001
1.50001
X-Axis
1.25001
Y-Axis
1
1.00001
1.25
0.750001
0.750005
1.5
0.500005
0.750001
Z-Axis
1
Y
X
1
Z-Axis
1.25
Z
1.25
1.5
1.5
X-Axis
1.5
1.25
0.500005
1
0.750001
0.5
0.500005
0.5
0.750005
0.750005
Y-Axis
1.00001
0.750001
Y-Axis 1.00001
1.25001
Z-Axis
1.25001
1
1.50001
1.5
1.25
X-Axis
1
1.50001
0.750001
0.5
0.5
1.25
0.750001
1.5
1
Z-Axis
1.25
X
1.5
Y
Z
c. 175
C d. 175
C
1.50001
0.5
0.5
1.25001
0.750001
Y-Axis
X-Axis
0.750001
1.00001
1
0.750005
1.25
Z-Axis
1
0.500005
1.5
0.5
0.5
1.50001
1.25
0.750001
1.25001
0.750001
X-Axis
1
Y-Axis
1.00001
1.5
1.25
0.750005
1
Z-Axis
1.5
0.500005
1.25
Y
1.5
X
Z
Y-Axis
0.500005
0.750005
1.00001
0.5
1.25001
1.50001
0.5
0.5 0.5
0.750001
0.750001
X-Axis
X-Axis
1
1
0.750001
0.750001
1.25
1.25
1.5
Y-Axis
0.500005
0.750005
1.5
1.00001
1.25001
1.50001
Z-Axis
1
1
Z-Axis
1.25
1.25
Y
X
1.5
1.5
Z
e. 150
C f. 150
C
Figure 4.17: Damage propagation across the FE. Grey cells are
m(B)
= deadD
and black cells are
m(B)
= deadB. (Continued on the next page).
1.50001
1.5
1.25001
Y-Axis
1.25
1.00001
X-Axis
0.750005
1
0.500005
1.5
0.750001
1.25
0.5
1.50001
0.5
X-Axis
1
1.25001
Y-Axis 1.00001
0.750001
0.750001
0.750005
0.5
0.500005
Z-Axis
0.5
1
X
Y
0.750001
1.25
Z
1
Z-Axis
1.25
1.5
0.500005
1.5
0.750005
Y-Axis
1.25
1.00001
X-Axis
1.25001
1
1.50001
1.5
0.750001
1.25
0.5
0.500005
0.5
X-Axis
1
0.750005
1.00001 Y-Axis
0.750001
0.750001
1.25001
0.5
1.50001
Z-Axis
1
0.5
1.25
0.750001
X
1
Z-Axis
Y
1.25
Z
1.5
g. 100
C h. 100
C
0.500005
1.5
0.750005
1.25
Y-Axis
X-Axis
1.00001
1
1.25001
0.750001
1.50001
1.5
0.5
0.500005
0.5
1.25
0.750005
X-Axis
0.750001
1
Y-Axis 1.00001
0.750001
Z-Axis
1.25001
1
0.5
1.50001
0.5
1.25
0.750001
X
1
Z-Axis
Y
1.25
Z
1.5
1.5
0.500005
1.25
0.750005
X-Axis
Y-Axis
1
1.00001
0.750001
1.25001
0.5
1.50001
0.500005
0.5
1.5
0.750005
1.25
0.750001
Y-Axis 1.00001
X-Axis
1
Z-Axis
1
1.25001
0.750001
1.50001
0.5
1.25
0.5
0.750001
X
1
Z-Axis
1.25
Y
Z
1.5
i. 50
C j. 50
C
1.50001
0.5 0.5
1.25001
Y-Axis 1.00001
0.750001
0.750005
0.750001
X-Axis
0.500005
0.5
1
0.5
Z-Axis
1
0.750001
1.25
0.750001
X-Axis
1.25
1
1.5
1.50001
1
1.25001
Z-Axis
1.25
1.00001 Y-Axis
1.5
0.750005
1.25 1.5
0.500005
Y
1.5
X Z
Y-Axis
1.50001
1.25001
0.5
1.00001
0.5
0.750005
0.500005
0.5 0.5
0.750001
0.750001
0.750001
X-Axis
1
0.750001
X-Axis
1
1.25
Z-Axis
1
1.25
1
Z-Axis
1.5
1.50001
1.25001
Y-Axis 1.00001
1.5
0.750005
0.500005
1.25
1.25
1.5
1.5
Y
X
Z
k. 0
C l. 0
C
Figure 4.17: Continued.
the microcracks which were arrested due to a drop of applied stress or due to a
very high fracture stress in all neighbourhood grains (cells). These microcracks
have been observed in experiments (Lin et al., 1987; Nohava et al., 2002). The
overall similarity between these four brittle CA states is clear. That is what one
would expect after seeing the results in Figures 4.13 4.16.
The two brittle CA states at 150
C (Figures 4.17.e and 4.17.f) are quite

dierent. There is a large brittle crack in Figure 4.17.e; however it did not cross
the whole of the FE and the nal failure was due to ductile fracture. In contrast,
the brittle cleavage plane crossed virtually the whole of the FE in Figure 4.17.f
and only two ductile cells failed. Consequently the brittle phase for the state
shown in Figure 4.17.e is Y
3
= 0.4 and for the state shown in Figure 4.17.f it is
Y
3
= 0.85 (see Figure 4.16).
With temperatures increasing to 100
C (Figures 4.17.g and 4.17.h), 50
C
(Figures 4.17.i and 4.17.j), and 0
C (Figures 4.17.k and 4.17.l), one can see an

increasing number of deadD cells. The dierence between the states at each
temperature is insignicant.
Thus the results on the micro (CA) scale (Figure 4.17) are complementary
to those of the macro (FE) scale (Figures 4.13 4.16).
Finally a special transition temperature diagram can be constructed by com-
bining the data from Figures 4.11.a, 4.12, 4.13 and 4.16. This diagram is shown
in Figure 4.18.
200 180 160 140 120 100 80 60 40 20 0
0
500
1000
1500
2000
2500
3000
3500
Temperature,
o
C
S
t
r
e
s
s
,

M
p
a
Mean
F
=1.89e+003
Median
F
=1.90e+003
Mode
F
=2.03e+003
STD
F
=4.77e+002
Max
F
=3.38e+003
Min
F
=3.69e+002
Lower Shelf Transition Upper Shelf

F(LS)

F(US)

Y0

I
TT range
TT
F
region
Figure 4.18: Transition temperature diagram.
This diagram illustrates schematically the inuence of the fracture stress
distribution on the transition temperature range.
The fracture stress distribution, assumed to be independent of temperature,
is shown on the left. The maximum principal stress line,
I
(T), is drawn slightly
higher than c
B

Y 0
(T) because the data in Figure 4.13 shows that brittle
fracture can happen after the onset of plasticity, so
I
>
Y 0
. The transition
range obtained from Figure 4.16 is approximately from 160
C to 90
C. By
establishing the intersection points of the transition temperature limits with the
I
(T) curve, the transition temperature range of
F
can be obtained:
F(US)
and
F(LS)
.
The number of cells which have the fracture stress lower than the upper
shelf limit,
F(US)
, is so small that it is very unlikely that many brittle cells m
with
m
F
<
F(US)
will be located near the fracture propagation path. Although
some brittle cells with particularly low
F
will still fracture even at the upper
shelf temperatures (see e.g. Figures 4.17.i 4.17.l), the number of these cells is
so small that they do not aect the fracture process signicantly.
On the other hand the lower shelf limit,
F(LS)
, is so high that it is very
likely that many cells m with
m
F
<
F(LS)
will be located near the fracture
propagation path.
Thus the number of brittle cells which can fracture is inversely related to
temperature in the transitional region but saturation is reached as temperature
approaches either shelf.
At very low temperatures (the lower shelf) the number of brittle cells which
can fail is so high that all but a few brittle cells located at the fracture propaga-
tion path will fail. Further decrease of temperature results in saturation and no
signicant change in fracture behaviour occurs. Similarly, at very high tempe-
ratures (the upper shelf) this number is so small that very few brittle cells will
fail. For all practical purposes this is 100% ductile fracture. The saturation is
thus reached and further increase of temperature does not change the picture.
As both the upper and the lower shelves can be described as saturation, the
levels of scatter are small at these temperatures.
The fracture behaviour is dierent at transitional temperatures. The number
of brittle cells which can fail is now temperature-dependent and the fracture
propagation behaviour is aected by the locations of these cells as well as the
total number of them. So there might be a situation when no suitable brittle
cell is located at a point of highest macro (FE) stress. This will lead to further
increase of macroscopic stress until the brittle fracture criterion of equation
(3.12) is met at a point where a suitable brittle cell is located. However, in the
next simulation the locations of the brittle cells which could fail will be dierent
and thus the whole of the fracture propagation path will change. This is the
reason for the high scatter in the transitional temperature range.
The diagramin Figure 4.18 suggests that there are two important parameters
of the fracture stress distribution,
F(LS)
and
F(US)
. If we assume that these are
true material parameters and that they do not change with temperature, then
the diagram will indicate that the increase of the maximum principal stress,
whether due to higher yield stress or due to high stress triaxiality, shifts the
transition temperature range right, towards higher temperatures. There is some
experimental evidence to support this point.
Results presented by Kohout (2001) show that the transition temperatures
for V-notched Charpy samples are higher than for U-notched ones. This happens
because the stresses below the V-notch are higher than below the U-notch.
Hertzberg reproduced results (Hertzberg, 1996, page 389) which show that
the transition temperature range obtained in drop weight (DWTT) and in dy-
namic (DT) tear tests for A541 Class 6 steel is signicantly higher than that
recorded with the Charpy test. The reason is that DWTT and DT may be
considered to be oversized Charpy samples (Hertzberg, 1996, page 387). As a
consequence there is much higher stress triaxiality in DWTT and DT samples
compared with the Charpy specimen.
Materials which exhibit elevation of the rst yield stress with strain rate pro-
duce a transition temperature shift when tested at dierent strain rates. Results
reproduced by Hertzberg for impact and slow-bend Charpy tests (Hertzberg,
1996, pages 392 and 393) show that the transition temperature range for the
impact Charpy test is higher that for the slow-bend test.
The analysis of the diagramin Figure 4.18 has one very important conclusion.
This is that the whole of the grain size distribution is necessary for proper
simulation of the transitional fracture behaviour. If all grains have the same
grain size, e.g. the mean grain size, then no transitional behaviour is possible, as
the model will show that all cells fail in either ductile or brittle mode, depending
on the simulation temperature.
Indeed the diagram of Figure 4.18 indicates that the higher the STD(
F
),
the wider is the transition temperature range. If STD(
F
) is very small, then
the values
F(LS)
and
F(US)
are very close and the transition range is very
narrow. If STD(
F
)=0, which is the case of using a uniform grain size, then
F(LS)
=
F(US)
and the width of the transitional region is zero. Eectively the
model in that case can only simulate the upper or the lower shelf behaviour.
The grain size distribution has such an important role in this example be-
cause there are no explicit initiation sites of brittle fracture, or rather every
cell is a potential initiation site as = 1. However, if < 1, as in the next
two examples (sections 4.1.3 and 4.2.1), then the distribution of the initiation
sites (fracture stress distribution of aliveC brittle cells) is a key factor aecting
transitional behaviour.
4.1.3 The Charpy test
This example illustrates the simulation of the Charpy impact test with the use
of the full CAFE model.
Figures 4.19 and 4.20 show the meshes of all bodies included in the model.
These are the specimen, the anvils and the BS EN 10045-1 (1990) striker. Only
the rst 35 mm of the striker tup is modelled.
A 0.15 friction coecient is adopted for all contact surfaces.
The full CAFE model is only used in the nite elements located near the
anticipated fracture propagation path (damage zone). The damage zone consists
of 900 C3D8R elements (HKS, 2001) which are shown in Figure 4.21. Majority of
these are 1mm1mm1mm cubic elements. However slightly smaller elements
are used near the root of the notch.
5 5 5 cell ductile (M
D
= 125) and 10 10 10 cell brittle (M
B
= 1000)
CA arrays were created for each nite element in the damage zone. Thus the
ductile damage cell size is L
D
1/5 = 0.2 mm and the brittle damage cell is
L
B
1/10 = 0.1 mm. The full CAFE model will therefore have 112500 ductile
and 900000 brittle CA cells.
The material simulated in this example is a laboratory control rolled TMCR
steel. The data for this steel is obtained from published papers and by pri-
vate communication (Bhattacharjee and Davis, 2002; Bhattacharjee et al., 2003;
1
2
3
Figure 4.19: The nite element mesh of the Charpy specimen, the anvils and
the striker.
Davis, 2003).
Figure 4.22 shows two typical illustrations of the brittle and the ductile
fracture surfaces obtained on the present TMCR steel. The average spacing
between the larger voids is approximately 100 micron. The average cleavage
facet size is approximately 50 micron. Thus the damage cell sizes used in this
example are roughly two times larger than the microstructural features to which
they are usually related. This was a conscious decision aimed at cutting the
simulation time. The author believes that such discrepancy can be dealt with
appropriately by the present CAFE model, e.g. by assigning the fracture stress
to each brittle cell based on the grain size distribution (section 3.2.2).
The chemical composition of this TMCR steel is shown in Table 4.1.
The rolling process starts at 1120
C from a 145 mm thick slab and includes

1
2
3
Figure 4.20: A three-dimensional view of the Charpy test model.
1
2
3
Figure 4.21: The mesh of the damage zone of the Charpy specimen containing
900 nite elements.
26 passes to a nish roll temperature of 717
C. The nal slab thickness was

approximately 30 mm. Figure 4.23 shows the microstructure of this steel at the
a. Ductile fracture b. Brittle fracture
Figure 4.22: Typical ductile and brittle fracture surfaces of this TMCR steel.
Figure a. courtesy of Davis (2003), Figure b. is reproduced from Bhattacharjee
and Davis (2002).
Figure 4.23: Microstructure of the TMCR steel. From Bhattacharjee and Davis
(2002).
mid-thickness location, 15 mm below the surface. The grain size distribution at
the mid-thickness plane is shown in Figure 4.24.
Figures 4.23 and 4.24 suggest a duplex distribution. There is an obvious
drop in Figure 4.24 at d
g
1.1 micron which separates the two parts of the
C Si Mn P S Al Nb V
0.1 0.31 1.42 0.017 0.005 0.046 0.045 0.046
Table 4.1: Chemical composition of the TMCR steel used (weight %). From
Bhattacharjee and Davis (2002).
0 5 10 15 20 25 30
0
20
40
60
80
100
120
Mean =5.48e+000
Median =4.79e+000
Mode =5.71e001
STD =3.90e+000
Max =2.59e+001
Min =5.20e001
d
g
, micron
N
u
m
b
e
r

o
f

g
r
a
i
n
s
Figure 4.24: The 500-bin histogram of the grain size, d
g
. Data courtesy of Davis
(2003).
histogram. The two parts of the full histogram split at d
g
= 1.1 micron are
shown in Figures 4.25 and 4.26.
It must be said that Dr C Davis, who provided the raw grain size data
used in the histogram shown in Figure 4.24, suggested that high number of
grains with d
g
< 1 micron might be a by-product of the measuring technique.
However, it is not obvious what is the minimum reliable grain size (Davis,
2003). For this reason the grain size data was used in this example exactly as
provided without any ltering. Moreover, as will be shown below, grains with
d
g
1.1 micron do not contribute to the fracture propagation due to very high
fracture stresses. Nevertheless, this example demonstrates a useful technique
for simulating duplex grain size distributions.
We shall approximate each part of the full histogram with a separate Weibull
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2
0
20
40
60
80
100
120
Mean =7.04e001
Median =6.55e001
Mode =5.32e001
STD =1.71e001
Max =1.10e+000
Min =5.20e001
d
g
, micron
N
u
m
b
e
r

o
f

g
r
a
i
n
s
Figure 4.25: The 15-bin histogram of d
g
1.1 micron.
0 5 10 15 20 25 30
0
5
10
15
20
25
30
35
Mean =6.09e+000
Median =5.29e+000
Mode =4.20e+000
STD =3.73e+000
Max =2.59e+001
Min =1.13e+000
d
g
, micron
N
u
m
b
e
r

o
f

g
r
a
i
n
s
Figure 4.26: The 500-bin histogram of d
g
> 1.1 micron.
distribution based on the mean and the standard deviation values of this part.
The d
g
1.1 micron part of the histogram (Figure 4.25) is simulated with a
Weibull distribution with parameters W
= 1.065, W
= 0.189 and W
= 0.52
(W
l
distribution). The d
g
> 1.1 micron part of the histogram (Figure 4.26) is
simulated with a Weibull distribution with parameters W
= 1.298, W
= 5.401
and W
= 1.1 (W
r
distribution).
There are 4087 data points in the full histogram (Figure 4.24), 464 data
points in the d
g
1.1 micron histogram (Figure 4.25) and 3623 data points in
the d
g
> 1.1 micron histogram (Figure 4.26). Therefore the fraction of cells for
which the grain size will be generated with the W
l
distribution is 464/3623 =
0.128. For the rest of the brittle cells the W
r
distribution will be used.
The histogram of the brittle CA cell grain sizes generated with the W
l
and
the W
r
distributions is shown in Figure 4.27.
0 5 10 15 20 25 30 35 40 45
0
0.5
1
1.5
2
2.5
x 10
4
Mean =5.40e+000
Median =4.53e+000
Mode =5.61e001
STD =4.04e+000
Max =4.18e+001
Min =5.20e001
d
g
, micron
N
u
m
b
e
r

o
f

c
e
l
l
s
Figure 4.27: The 1000-bin histogram of the generated grain size.
It is much smoother that the experimentally obtained histogram shown in
Figure 4.24. This is partly due to a much larger number of cells (900000) than
the number of grains for which size was obtained experimentally (4087). On the
other hand the Weibull distribution is only an approximation of the very com-
plex experimental data. However, statistical characteristics of the histograms
shown in Figures 4.24 and 4.27 are quite similar.
Wu and Davis (2003) reported an eective fracture surface energy of 52 J/m
2
for a TMCR steel very similar to that used in this work. This value agrees well
with typical surface energy values (Lin et al., 1987).
Taking this value as the eective surface energy of a ferrite ferrite interface,
ff
, we can generate the fracture stress values using the grain size histogram
(Figure 4.27) according to equation (2.26), page 28. The histogram of the
generated brittle CA cell fracture stress is shown in Figure 4.28.
0 1000 2000 3000 4000 5000 6000 7000 8000 9000
0
500
1000
1500
2000
2500
3000
Mean =3.39e+003
Median =2.81e+003
Mode =2.22e+003
STD =1.74e+003
Max =8.31e+003
Min =9.70e+002

F
, MPa
N
u
m
b
e
r

o
f

c
e
l
l
s
Figure 4.28: The 1000-bin histogram of the generated fracture stress,
F
.
As in the example of section 4.1.1 a two-parameter Weibull probability den-
sity function was used to simulate the distribution of the initial void volume
fraction, f
0
. The shape parameter was taken as W
= 2 and the scale parame-

ter was taken as W
= 2.82 10
4
. The resulting histogram of f
0
is shown in
Figure 4.29.
The value of the misorientation threshold,
F
, is chosen based on the experi-
mental data reported by Bhattacharjee and Davis (2002). The authors observed
misorientation angles as high as 12
within a single cleavage facet. On the other

hand the angles between the facets were reported in the range 17
45
. Ac-
cordingly
F
= 30
is used in this example.

The maximum possible misorientation angle for this steel is
max
60
(Bhattacharjee and Davis, 2002; Bhattacharjee et al., 2003). Therefore the ori-
entation angle, , assigned to each brittle CA cell, is generated using a uniform
distribution, [0 . . .
max
].
0 0.2 0.4 0.6 0.8 1 1.2
x 10
3
0
500
1000
1500
2000
2500
3000
Mean =2.50e004
Median =2.35e004
Mode =1.92e004
STD =1.31e004
Max =1.02e003
Min =2.38e007
f
0
N
u
m
b
e
r

o
f

c
e
l
l
s
Figure 4.29: The initial void volume fraction histogram.
The value = 510
3
is used in this example. So only 510
3
900000 =
4500 randomly chosen brittle CA cells can initiate brittle fracture.
The values for other model parameters are c
D
= 1.4, c
B
= 3, X
max
(D)
= M
2/3
D
and X
max
(B)
= M
2/3
B
.
As in the previous example the power hardening law of the form expressed by
equation (4.1) was used. The temperature dependence of the rst yield stress,
Y0
, and the hardening exponent, n, for this TMCR steel is shown in Figure
4.11. Both curves are tted over the data extracted from the tensile test results
provided by Davis (2003).
One useful feature of the present CAFE models (both the full and the sim-
plied) is that each fracture mechanism can be switched on and o according
to users wish. This feature is used in this example to tune the parameters of
the Rousselier ductile damage model. By switching o the brittle failure we can
ensure that all fractures will occur in the ductile CA array.
If both the brittle and the ductile fracture mechanisms are switched on then
the running time of this example is about one week. If the brittle fracture is
switched o then the running time can be cut roughly by a factor of three, to
two three days.
However, even two days is a long time for tuning as tens of runs might be
necessary to nd the best combination of modelling parameters. Accordingly
a simplied tuning procedure was adopted in this example to cut the number
of simulations (section 2.4). All three Rousselier model parameters, D,
1
and
F
, were tuned together. The tuning was aimed at achieving the total energy
absorbed similar to that obtained in the Charpy test at the upper shelf tempe-
ratures.
A more thorough tuning could have been performed if experimental force
time Charpy data were available. This sort of data is usually obtained with the
instrumented Charpy test (Shterenlikht et al., 2003). Unfortunately no such
data were available for this steel.
The Charpy impact data from Bhattacharjee et al. (2003) shown in Figure
4.39 was used for tuning. The experimentally obtained upper shelf temperatures
for this TMCR steel are T 20
C and the corresponding Charpy energy values

are C
v
180J (Figure 4.39).
The best tted values for the Rousselier model parameters are D = 3,
1
=
500 MPa and
F
= 8.
After the tuning the simulation of the Charpy test at T = 50
C was
performed. The simulation results on the macro (FE) scale are shown in Figures
4.30 and 4.31. Accordingly Figures 4.32 4.35 demonstrate results on the micro
(CA) scale.
Figure 4.30.a shows the modelling force curve. As no instrumented Charpy
test data were available for this TMCR steel, Figure 4.30.b shows a typical force
curve obtained during the instrumented Charpy test on a TMCR steel dierent
from the one used in the present work. This experimental data is courtesy of
Davis (2003) and is given here only for qualitative comparison of force drop in
the model and in the experiment. It is very likely that the experimental force
curve obtained on the TMCR steel used in this example would exhibit a similar
drop to that shown in Figure 4.30.b.
It is easy to see that brittle fracture propagates slower in the model than in
the experiment. Brittle fracture propagation can be easily identied in the ex-
perimental force curve (Figure 4.30.b) as an instantaneous drop of force from the
maximum value to almost zero. The model, however, exhibits a more gradual
0 2 4 6
0
5
10
15
20
Time, msec
F
o
r
c
e
,

k
N
0 2 4 6
0
5
10
15
20
Time, msec
F
o
r
c
e
,

k
N
For qualitative
comparison only
a. Simulation b. A typical experiment
Figure 4.30: Modelling (T = 50
C) and experimental (T = 60
) Charpy
force curves. Experimental data obtained on a dierent TMCR steel is courtesy
of Davis (2003).
reduction in force.
The present CAFE model cannot simulate a sharp drop of force because the
fracture cannot propagate from one nite element into another. This is a limita-
tion imposed by the organisation of the Abaqus code (section 5.1). Thus brittle
fracture must reinitiate in each nite element (each brittle CA) in the fracture
propagation path which results in a more gradual brittle fracture propagation
process.
The total energy absorbed predicted by the model is 89J. The energies ob-
tained in the experiment (Figure 4.39) at this temperature range from 117J to
142J. Therefore the full CAFE model probably underpredicts the total energy
absorbed at T = 50
C. However, at least three simulations have to be per-

formed before some indication of the scatter of simulated energy values can be
obtained. Because of the long running times it was not feasible to do this with
the full CAFE model.
A deformed mesh at the end of the simulation is shown in Figure 4.31. The
failed FEs are removed from the mesh.
The visualisation of results on the micro (CA) scale is not straightforward
because all CAs are cubic but some FEs are not. To avoid the overlapping of the
CA blocks during visualisation they are shown shrunk. However the position
and the size of each CA is related to the position and size of the corresponding
FE. This technique is illustrated in Figure 4.32. The brittle CA arrays are used
1
2
3
Figure 4.31: Deformed FE mesh of the Charpy test simulation at T = 50
C.
for visualisation in this example.
Figure 4.33 shows all 4500 AliveC brittle cells generated for this example.
Only these cells can initiate the brittle fracture.
Fracture propagation on the CA scale is shown in Figure 4.34. Only DeadB
(black) and DeadD (green) cells are shown. So the brittle fracture areas are
black and green areas represent the ductile fracture. Fracture propagates on the
XZ plane along direction X.
Fracture starts in the ductile mode. However several small brittle cracks
are also visible in Figures 4.34.a and 4.34.b. It is easy to see that the fracture
1
2
3
a. FE mesh of the damage zone
Y
X Z
0.245664
-4.51959 -4.51959
-2.13697
5.01092
2.62829
0.245664
-2.13697
2.62829
5.01092
b. Corresponding brittle CA blocks
Figure 4.32: FE mesh of the damage zone and the corresponding brittle CAs.
propagation front is not a straight line but rather an arc (Figures 4.34.c
4.34.g). The brittle fracture starts after some ductile crack extension (Figure
4.34.d) and stops when the remaining ligament is very small and the dominant
-4.51959
0.500054
0.5
-2.13697
2.90621
Y-Axis
X-Axis
0.245664
5.31237
2.90406
2.62829
7.71852
5.01092
10.1247
-4.51959
Z-Axis
2.90621
5.30812
-2.13697
5.31237
Y-Axis 0.245664
2.90406 7.71852
2.62829
7.71217
10.1247
5.01092
5.30812
10.1162
Y
X
7.71217
Z
10.1162
Figure 4.33: Locations of AliveC brittle cells.
fracture mode is ductile plastic collapse (Figures 4.34.g and 4.34.h).
This is a very realistic sequence of fracture propagation events.
The fracture surface at the end of the simulation is shown in Figure 4.35.
Photographs of three broken Charpy samples of this TMCR steel at T = 50
are shown in Figures 4.43.c 4.43.e. The simulated percentage of the brittle
phase, 54%, is within the range of the experimental values which are from 45%
to 75%.
The regions of initial ductile fracture, of brittle fracture and of the nal
ductile fracture, associated with the plastic collapse in the remaining ligament,
which are found in the experimental fracture surfaces at this temperature (Fig-
ures 4.43.c 4.43.e.) are adequately reproduced in the simulated fracture surface
(Figure 4.35). Moreover the shapes of these regions are not too dissimilar espe-
cially if we remember that the simulated fracture surface is drawn on the initial
(undeformed) mesh geometry.
One point of dierence between the simulated and the experimental fracture
surfaces is that the remaining ligament in the model is unbroken. The nite
elements in the last raw of the fracture propagation plane are still alive at the
end of the simulation. Figure 4.31 shows that these FEs are highly distorted
0.5
-4.51959
0.500054
2.90406
-2.13697
2.90621
Y-Axis
Z-Axis
5.30812
0.245664
X-Axis
5.31237
7.71217
10.1162
7.71852
Z
Y
X
0.5
-4.51959
0.500054
2.90406
-2.13697
2.90621
Y-Axis
Z-Axis
5.30812
0.245664
X-Axis
5.31237
7.71217
10.1162
7.71852
Z
Y
X
a. 3.6 10
4
sec b. 4.8 10
4
sec
0.5
-4.51959
0.500054
2.90406
-2.13697
2.90621
Y-Axis
Z-Axis
5.30812
0.245664
X-Axis
5.31237
7.71217
10.1162
7.71852
Z
Y
X
0.5
-4.51959
0.500054
2.90406
-2.13697
2.90621
Y-Axis
Z-Axis
5.30812
0.245664
X-Axis
5.31237
7.71217
10.1162
7.71852
Z
Y
X
c. 6 10
4
sec d. 7.2 10
4
sec
0.5
-4.51959
0.500054
2.90406
-2.13697
2.90621
Y-Axis
Z-Axis
5.30812
0.245664
X-Axis
5.31237
7.71217
10.1162
7.71852
Z
Y
X
0.5
-4.51959
0.500054
2.90406
-2.13697
2.90621
Y-Axis
Z-Axis
5.30812
0.245664
X-Axis
5.31237
7.71217
10.1162
7.71852
Z
Y
X
e. 8.4 10
4
sec f. 9.6 10
4
sec
0.5
-4.51959
0.500054
2.90406
-2.13697
2.90621
Y-Axis
Z-Axis
5.30812
0.245664
X-Axis
5.31237
7.71217
10.1162
7.71852
Z
Y
X
0.5
-4.51959
0.500054
2.90406
-2.13697
2.90621
Y-Axis
Z-Axis
5.30812
0.245664
X-Axis
5.31237
7.71217
10.1162
7.71852
Z
Y
X
g. 1.32 10
3
sec h. 6 10
3
sec
Figure 4.34: Fracture propagation on the CA scale at T = 70
C. Only DeadB
(black) and DeadD (green) cells are shown.
suggesting high shear and very low mean stresses in these elements. Accordingly
these elements still carry some load bearing capacity as the Rousselier damage
Z
X
Z-Axis
7.71217
10.1162
0.5
2.90406
5.30812
Y
10.1162
5.30812
7.71217
2.90406
0.5
Z-Axis
Figure 4.35: The simulated Charpy fracture surface at T = 50
C. The brittle
phase is 54%.
model can only account for a volumetric void growth, but not for a shape change
(section 2.1.6).
An unbroken remaining ligament of 0.1 0.25 mm is sometimes observed
in the tested Charpy samples of this TMCR steel though usually at higher
temperatures.
The results of this example show that the full CAFE model is capable of
simulating transitional ductile-brittle fracture in a Charpy specimen. On the
micro (CA) scale the model can predict the percentage of the brittle phase, the
locations and shapes of the ductile and the brittle fracture regions (Figure 4.35)
and the crack front throughout the simulation (Figure 4.34). The force time
data, Figure 4.30.a, and the total energy absorbed are obtained on the macro
(FE) scale.
As was said earlier the major problem of the present CAFE model is the
inability to simulate brittle cracks running across the nite element boundaries.
Such is the limitation of the present realisation of a CAFE approach via the
VUMAT user subroutine and the Abaqus code. Accordingly the brittle fracture
4.2. THE SIMPLIFIED CAFE MODEL 95
has to reinitialise in each nite element in the fracture propagation path. Thus
the brittle zone in the simulated Charpy fracture surface consists of many small
cracks rather than of a single large crack (Figure 4.35). On the macro-scale
such model behaviour results in a gradual reduction of force as opposed to an
instantaneous drop obtained experimentally (Figure 4.30).
Slow brittle fracture propagation can be also caused by a high misorientation
threshold,
F
. The value used in this example,
F
= 30
, is only a rough guess

based on very limited experimental data (Bhattacharjee and Davis, 2002).
This example also shows how the fracture stress distribution (Figure 4.28)
is obtained from the grain size data (Figure 4.24). Although in this example
the separation of the full grain size histogram into two parts probably was not
necessary, because the cells with
F
> 5000 MPa did not take part in the fracture
process, this is a useful technique that can be explored in future (section 5.1).
4.2 The simplied CAFE model
4.2.1 The Charpy test
In this example the simplied CAFE model is used to simulate the scatter
usually obtained experimentally in the Charpy test at transitional temperatures.
The ability of the full CAFE model to predict transitional behaviour and
scatter in terms of the percentage of the brittle phase was demonstrated for a
single-FE model in section 4.1.2 (Figure 4.16). In this example the Charpy test
model described in section 4.1.3 with the following parameters was used.
The initial void volume fraction for each ductile CA cell was taken a f
0
=
1 10
4
.
The critical value of the damage variable,
F
, was distributed across the duc-
tile CA arrays using the normal distribution with

F
= 8 and STD(
F
) = 1.2.
The

F
value is based on the best-tted critical damage variable obtained in
section 4.1.3. A typical
F
histogram is shown in Figure 4.36.
The other two Rousselier model parameters were tuned following the pro-
cedure described in section 4.1.3 and the best-tted values were found to be
1
= 517 MPa and D = 2.8.
The best-tted values for this parameters in the full CAFE model were
2 4 6 8 10 12 14
0
50
100
150
200
250
300
350
400
450
Mean =8.00e+000
Median =7.99e+000
Mode =8.04e+000
STD =1.20e+000
Max =1.33e+001
Min =2.92e+000

F
N
u
m
b
e
r

o
f

c
e
l
l
s
Figure 4.36: A typical histogram of
F
.
1
= 500 MPa and D = 3. Thus the dierence between the full and the
simplied models is small in this respect.
A temperature-dependent misorientation threshold,
F
, is used in this ex-
ample. The following temperature dependence is used:
F
=
_
max
if T 80
C
a T +b if 80
C < T < 20
C
0 if T 20
C
(4.2)
where a and b are chosen to ensure continuity of
F
(T), as shown in Figure 4.37.
The value of the concentration factor for the brittle CA array was increased
to c
B
= 11 in this example in an attempt to promote faster brittle fracture prop-
agation. The value c
B
= 3 used in the previous example (section 4.1.3) probably
was not high enough and as a consequence the brittle fracture propagation was
too slow (Figure 4.30.a).
Also the fraction of brittle cells capable of initiating brittle fracture was
increased two times compared to the full CAFE model, to = 0.01. So the
number of brittle AliveC cells was 9000 in the simplied CAFE model.
All other model parameters, which are not mentioned here, have exactly the
80 70 60 50 40 30 20
0
20
40
60
T,
o
C

F
,

o
Figure 4.37: The misorientation threshold,
F
(T).
same values as in section 4.1.3.
In all simulations of this example variable mass scaling was applied to par-
ticularly small nite elements which form contact surfaces for the specimen
anvil interactions (Figure 4.19). A stable time increment thus was increased
approximately ten times whereas the change in total energy of the whole model
was less than 0.1%. Accordingly the total simulation times were cut by a fac-
tor of ten. No mass scaling was applied to the nite elements in the damage
zone (Figure 4.21). The running time of each simulation in this example was
approximately ten hours.
The simulations were performed at several temperatures from T = 80
C
to T = 0
C, three runs at each temperature.

The results are shown in Figures 4.38 4.43. Results on the macro (FE)
scale are shown in Figures 4.38 4.40. Accordingly Figures 4.41 4.43 show
results on the micro (CA) scale.
The value of 50% impact transition temperature (ITT) shown in Figure 4.38
is based on the fracture surface appearance. Data presented by Bhattacharjee
et al. (2003) shows that for this TMCR steel 50% brittle phase is obtained in
the Charpy tests at T 50
C.
The simulation results at T = 0
C and at T 75
C agree very well with

the experimental data. However, for other temperatures the model predicts a
100 80 60 40 20 0 20
0
10
20
30
40
50
60
70
80
90
100
T,
o
C
%

B
r
i
t
t
l
e

p
h
a
s
e
CAFE model
Experiement
50% ITT
Figure 4.38: Charpy percentage of brittle phase. 50% ITT is taken from Bhat-
tacharjee et al. (2003). Experimental data was provided by Corus UK Ltd.
more sharp transition compared with that obtained experimentally. The 50%
ITT on the simulated data is around T = 30
C.
The scatter of the simulated brittle phase values is higher in the transitional
region than at the upper and the lower shelf. At the upper shelf, at T 20
C,
the three values of the brittle phase obtained at each temperature are virtually
identical. At the lower shelf, at T 60
C, the maximum dierence between

the three values is approximately 3% at T = 75
C. However at transitional
temperatures the simulated scatter is much higher. The maximum dierence
between the three simulated brittle phase values is about 21% at T = 30
C
and T = 35
C. The experimental scatter can only be assessed at T = 50
C
where the maximum and the minimum values dier by 67%.
The experimental transition temperature range is from T = 75
C to T =
0
C, however, the upper bound value is arguable because there is no clear upper
shelf. The simulated data suggests that the modelling transitional temperature
range is less wide, approximately from T = 60
C to T = 20
C.
Finally the simulated transition in Figure 4.38 towards the lower shelf is
smoother than towards the upper shelf. This result is in contrast to that ob-
tained in section 4.1.2 for the single-FE model, Figure 4.16. Moreover the lower
shelf behaviour in Figure 4.16 is similar to the upper shelf behaviour in Figure
4.38 in that the scatter is very low in both cases. The limit values of the brittle
phase (100% for the lower shelf in Figure 4.16 and 0% for the upper shelf in
Figure 4.38) were achieved in both cases.
So the single-FE model predicts 100% brittle phase at the lower shelf tem-
peratures, but does not predict 0% brittle phase at the upper shelf, Figure 4.16.
In contrast the simplied model of the Charpy test does not predict 100% brit-
tle phase at the lower shelf, but predicts 0% brittle phase at the upper shelf,
Figure 4.38. The brittle phase prediction of the simplied model of the Charpy
test is probably more physically based and it does agree very well with the
experimental data at the lower and the upper shelf.
Such a dierence of the CAFE model performances in these two examples is
probably caused by dierent values of the CAFE model parameters and by the
fact that fracture cannot cross the nite element boundary due to the limitations
of the Abaqus code (section 5.1).
Figure 4.39 shows the experimental and the simulated values of the total
energy absorbed in the Charpy test.
The model can predict the upper shelf energy values reasonably well, however
the lower shelf and the transition temperature range are quite dierent from
the experimental data. The simulated lower shelf is achieved at T 60
and the lower shelf energies are about 60J. There is no clear lower shelf in
the experimental data. However, it is reasonable to assume that it starts at
T 80
and the energies there are 20J 40J.

The simulated transition temperature range according to the energy data
is the same as the range obtained from the brittle phase results, from T =
60
C to T = 20
C. However, the experimental transition temperature range

obtained from Figure 4.39 is from T = 80
C to T = 20
C, which is slightly
lower than the range that can be perceived from Figure 4.38.
The simulated scatter in energy values at the transitional temperatures is
higher than at the lower or at the upper shelf. The maximum dierence between
100 80 60 40 20 0 20
0
20
40
60
80
100
120
140
160
180
200
T,
o
C
E
n
e
r
g
y
,

J
CAFE model
Experiment A
Experiment B
Figure 4.39: Charpy transition energy data. Experimental data A is reproduced
from Bhattacharjee et al. (2003), and data B is courtesy of Davis (2003).
the three simulated energy values at the lower shelf is 14% at T = 70
C. At
the upper shelf the maximum dierence is only 1% at T = 20
C. However, at
T = 35
C the maximum dierence is 60%.

The experimental data shows very high scatter at T = 60
C (40% vari-
ation) and at T = 75
C (more than 5 times). However, these experimental

results must be compared with some caution as they were obtained at dierent
times on dierent machines by dierent people, so there is a possibility that
dierent experimental practices could have contributed to such a large scatter.
The energy values predicted by the simplied CAFE model at T = 50
C,
66J 82J (Figure 4.39) are similar to that generated by the full CAFE model,
89J, section 4.1.3.
The lateral expansion measured on the deformed meshes and on the broken
Charpy samples is presented in Figure 4.40.
On the whole the simulated lateral expansion data resembles the simulated
energy, Figure 4.39. The transition temperature range in Figure 4.40 is from
100 80 60 40 20 0 20
0
0.5
1
1.5
2
2.5
3
T,
o
C
L
a
t
e
r
a
l

E
x
p
a
n
s
i
o
n
,

m
m
CAFE model
Experiement
Figure 4.40: Simulated and experimental Charpy lateral expansion data.
T = 60
C to T = 20
C. The scatter at the upper shelf is very low, less than

1%. The scatter at the lower shelf is higher, 19% at T = 80
C. However, the
scatter is much higher at the transitional range, from T = 60
C to T = 35
C.
The maximum dierence between the three simulated lateral expansion values
is 51% at T = 35
C.
The simulated lateral expansion values agree well with those obtained from
the broken Charpy samples using a digital calliper. However, the scatter in
the experimental values is smaller than that predicted by the simplied CAFE
model at T = 50
C.
Figure 4.41 shows four simulated fracture surfaces obtained at the lower
shelf temperatures. Because of the very small level of scatter only one of three
simulated surfaces is shown for each temperature. As in section 4.1.3 deadD
cells are shown green and black cells are deadB.
The four fracture surfaces shown in Figure 4.41 dier very slightly. This is
an additional evidence in support of the point that the lower shelf behaviour
is achieved at T 60
C. All four fracture surfaces show very little areas of

a. 80
C, 96% brittle b. 75
C, 96% brittle
c. 70
C, 94% brittle d. 60
C, 94% brittle
Figure 4.41: The simulated Charpy fracture surfaces at the lower shelf tempe-
ratures.
ductile fracture, most of them are located immediately below the root of the
notch. However, some isolated islands of ductile fracture can be seen far from
the notch root in all four surfaces.
Experimental Charpy fracture surfaces obtained at T = 70
C (Figure
4.43.a) and at T = 60
C (Figure 4.43.b) show lower values of the brittle

phase, 80% and 60% accordingly. This is primarily due to signicantly large
ductile shear regions. As was said earlier the present CAFE model cannot sim-
ulate shear fracture because the Rousselier damage model can only account for
volumetric void growth (section 5.1). Nevertheless there is a qualitative agree-
ment between the simulated and the experimental fracture surfaces at the lower
shelf.
Two of the simulated fracture surfaces at the upper shelf are shown in Figure
4.42. Both show virtually 100% ductile fracture although very few isolated
islands of brittle fracture can be seen.
a. 20
C, 0% brittle b. 0
C, 0% brittle
Figure 4.42: The simulated Charpy fracture surfaces at the upper shelf tempe-
ratures.
Experimental Charpy fracture surfaces obtained at T = 20
C (Figure
4.43.g) and at T = 0
C (Figure 4.43.h) have 10% and 0% brittle phase accord-

ingly. These experimental fracture surfaces exhibit substantial delaminations
which cannot be modelled at present.
Finally Figure 4.43 shows the simulated fracture surfaces at transitional
temperatures. Three surfaces are shown for each temperature to illustrate sig-
nicant variation from one simulation to another.
At T = 50
C and T = 40
C the dominant fracture mode predicted by

the simplied CAFE model is still brittle, Figures 4.43.a 4.43.f. However,
there are more islands of ductile fracture compared with the surfaces shown in
Figure 4.41. At the same time the locations and sizes of the ductile areas vary
much more at T = 50
C and T = 40
C than at the lower shelf temperatures,

Figure 4.41.
The simplied model predicted higher fraction of brittle phase at T = 50
C
(Figures 4.43.a 4.43.c) that that obtained with the full CAFE model at the
same temperature, Figure 4.35. This is most probably due to dierent values of
a. 50
C, 92% brittle
c. 50
C, 88% brittle
Figure 4.43: The simulated Charpy fracture surfaces at transitional temperatu-
res. (Continued on pages 105 106).
modelling parameters, e.g. c
B
, used in the full and the simplied CAFE models.
The simulated fracture surfaces at T = 50
C and T = 40
C are quite far

from the experimental ones shown in Figures 4.43.c 4.43.f. This is because the
simulated percentage of brittle phase is substantially higher than that obtained
experimentally.
At temperatures T = 35
C and T = 30
C, Figures 4.43.g 4.43.l, the

variation from one simulation to another is greater than at T = 50
C and
T = 40
C.
For example at T = 35
C Figure 4.43.g shows more islands of ductile

fracture towards the back side of the Charpy sample, whereas in Figure 4.43.h
e. 40
C, 89% brittle f. 40
C, 80% brittle
g. 35
C, 74% brittle h. 35
C, 72% brittle
i. 35
C, 87% brittle j. 30
C, 57% brittle
k. 30
C, 69% brittle l. 30
C, 60% brittle
m. 25
C, 25% brittle n. 25
C, 28% brittle
o. 25
C, 22% brittle
a. 70
C, 60% brittle
c. 50
C, 70% brittle
Figure 4.43: Fracture surfaces of the broken Charpy samples. The brittle phase
values are courtesy of Corus UK Ltd. (Continued on the next page).
ductile fracture areas are located closer to the notch of the specimen. However,
the nal fractions of brittle phase for these two fracture surfaces are very similar,
74% and 72% respectively. The third fracture surface shown in Figure 4.43.i has
signicantly higher brittle phase, 87%, and accordingly the ductile fracture zones
are mostly found next to the notch, although the is one ductile fracture island
located in the last (counting from the notch) row of the nite elements.
Similar variations can be found among simulated fracture surfaces at T =
e. 50
C, 45% brittle f. 40
C, 25% brittle
g. 20
C, 10% brittle h. 0
C, 0% brittle
30
C, Figures 4.43.j 4.43.l. In Figure 4.43.l there is a large brittle region in

the centre of the fracture surface surrounded by the areas of ductile fracture.
This fracture surface is quite close qualitatively to that shown in Figure 4.35.
Finally at T = 25
C, Figures 4.43.m 4.43.o, the variation from one

simulation to another is smaller than at T = 35
C and T = 30
C. There is
no pattern as to where the brittle fracture areas are located. No brittle crack
running across the whole of a nite element can be found in any of the three
fracture surfaces. This is caused by a relatively small value of the misorientation
threshold,
F
= 5
at T = 25
C, Figure 4.37.
Thus low
F
at the upper-transition and the upper shelf temperatures in-
hibits or stops brittle fracture propagation which otherwise would run very fast
across virtually any brittle cell due to very high value of the brittle concentration
factor, c
B
= 11. In contrast, at the lower shelf the brittle fracture can propagate
with very few deviations because
F
is so high that virtually all neighbouring
cells m and l will have |
m
l
| <
F
.
Thus it is probably true to say that in this example the simplied CAFE
model was able to simulate transitional behaviour largely due to the chosen
value of c
B
and the chosen temperature dependence of
F
, equation (4.2). If
these two model parameters had other values then the whole transition curve
would have looked dierently, or maybe there would have been no transitional
behaviour at all.
Chapter 5
Discussion
The results of the four examples of Chapter 4 show that both the full and the
simplied CAFE models are suitable for the modelling of transitional ductile-
brittle fracture in steels. The performance of the two models diers very slightly,
however, the simplied model is much faster.
The ability of the full model to correctly address the strain reversal eects
was demonstrated in section 4.1.1 (Figure 4.5). The simulation of strain reversal
is a delicate task because many technical details must be taken into account (e.g.
elastic unloading, anisotropy of the Rousselier damage model, accurate strain
rate decomposition, etc.). Thus among other results this example demonstrates
that the Rousselier damage model was coded correctly (Appendix B).
Both the full and the simplied models can predict realistic transitional
behaviour, including the levels of scatter, in terms of the percentage of brittle
phase (Figures 4.16, 4.34, 4.35, 4.38, 4.41, 4.42 and 4.43) or in terms of the
total energy absorbed in the Charpy test (simulation results of section 4.1.3 and
Figure 4.39).
It is interesting to note that some results indicate that transition towards the
upper shelf is smoother than that towards the lower shelf, e.g. Figures 4.16, 4.38
(experimental data) and Figure 4.39 (both experimental and simulated data).
Such a behaviour can be explained with the use of the diagram of Figure 4.18.
Dierent slopes of the transition curve towards the upper and the lower
shelves, either in terms of the energy absorbed or in terms of the brittle phase,
111
112 CHAPTER 5. DISCUSSION
could be caused by the fact that temperature dependence of the rst yield stress,
Y 0
(T), is such that
Y 0
T
2
> 0 (5.1)
within the range of temperatures for which experimental data was available
(Figure 4.18). Consequently
|
Y 0
(T + T)
Y 0
(T)| < |
Y 0
(T T)
Y 0
(T)|. (5.2)
Accordingly the fracture propagation behaviour will change more with decreas-
ing than with increasing temperatures.
However, the simulated brittle phase data in Figure 4.38 exhibits the op-
posite trend, that the transition towards the lower shelf is smoother than that
towards the upper shelf. Probably a temperature dependence of
Y 0
is only one
of the factors resulting in transition behaviour, and perhaps not the major one.
The transitional behaviour was achieved in the single-FE CAFE model in
section 4.1.2, Figure 4.16, solely due to the temperature dependence of
Y 0
.
This, of course, is only a qualitative prediction because the single-FE model
has the uniaxial stress state until the nal failure. The model does not account
for necking or the strain concentration associated with it. Such a deformation
history (the uniaxial stress until failure) is virtually impossible to reproduce
experimentally, hence the simulation results obtained on the single-FE model
are only suitable for qualitative analysis of the performance of the full CAFE
model.
However, the temperature dependence of
Y 0
was not enough to simulate a
realistic transitional behaviour in the Charpy test. The dierence between the
rst yield stresses at the upper and lower shelf temperatures is relatively small:
Y 0
(80
C) = 505 MPa and

Y 0
(0
C) = 453 MPa (Figure 4.11.a) and the

hardening exponent is virtually constant within this temperature range, Figure
4.11.b. Thus the temperature dependence of the misorientation threshold,
F
,
was introduced in equation (4.2), Figure 4.37. With this addition results ob-
tained with the simplied CAFE model were much closer to the experimental
data, Figures 4.38, 4.39 and 4.40.
There is little experimental evidence that
F
is temperature-sensitive. One
of the possible explanations for this is that usually misorientation analysis is
113
performed on fracture surfaces obtained at the lower shelf temperatures (Bhat-
tacharjee and Davis, 2002; Bhattacharjee et al., 2003).
Other authors also had to introduce additional (apart from the temperature
dependence of
Y 0
) temperature-dependent modelling parameters to make their
models predict transitional behaviour. Of these, the most popular is temper-
ature dependence of the scale parameter of Weibull distribution,
u
, equation
(2.28), page 29. Burstow (1998) and Burstow et al. (2003) reported the tuned
values of
u
from 1722.8 MPa at T = 45
C to 2699.9 MPa at T = 20
C. The
possibility of the temperature dependence of the shape parameter of the Weibull
distribution, m, (equation (2.28), page 29) is discussed by Burstow (2003).
Probably a similar eect could have been achieved if the concentration factor
for the brittle CA array, c
B
, was made temperature-dependent. However such
modication has no physical basis and, if introduced, will be just a modelling
trick.
Useful information about the fracture progression and change of the per-
centage of brittle phase during crack propagation can be obtained from Figures
4.6.b, 4.7, 4.14 and 4.15. These outputs, if required, can be easily requested
from any nite element in the damage zone of the Charpy model.
With the use of an appropriate visualisation technique many useful details of
fracture propagation can be revealed by plotting the states of either the ductile
or the brittle CA arrays. In this particular realisation of the CAFE model the
brittle CA arrays were designed for visualisation (section 3.2.2). Figures 4.8,
4.9, 4.10, 4.17, 4.34, 4.35, 4.41, 4.42 and 4.43 provide information about the
fracture propagation on the micro (CA) scale. The shape of the crack front, the
speed of fracture propagation, locations and shapes of the areas of brittle and
ductile fracture can be obtained from these Figures.
The ways in which material properties can be taken into account by both
models were demonstrated in Chapter 4. Particular attention should be drawn
towards the grain size distribution.
It was shown in section 4.1.3 how a duplex grain size distribution was simula-
ted by the full CAFE model. The strategy of applying dierent random number
generators to simulate dierent parts of the grain size histogram can be used
to simulate a grain distribution data of a much greater complexity. Moreover,
additional data (e.g. the histogram of the second phase particle sizes), if avail-
able, can be easily included into the model through the fracture stress histogram
following the technique used in section 4.1.3. The method used in section 4.1.3
is very simple. However, proper statistical methods are available if one wants
to use them (Devroye, 1987; McLachlan and Basford, 1988; Scott, 1992). These
methods, however, are signicantly more complex and require some knowledge
of probability and statistics.
TMCR steels were chosen to verify the CAFE model in Chapter 4 only
because the author had an access to the data generated by a thorough study of
one particular TMCR steel. However, the CAFE model described in Chapter 3
is designed to be able to simulate other widely used types of steels. Of particular
importance for the fracture mechanics community are line pipe steels, ship steels,
pressure vessel steels, nuclear reactor steels etc. Moreover the model can be
applied to other materials, e.g. aluminium, provided that cell properties and
state variables for the brittle and ductile CA arrays,
m(D)
,
m(B)
,
m(D)
and
m(B)
, are chosen according to the knowledge of the micromechanics of fracture
of this material.
As was said in section 3 the present CAFE model was built around C3D8R,
a nite element with a single integration point. This is the only 8-noded nite
element allowed in the Abaqus\Explicit program. However, if nite elements
with several integration points are used (e.g. in another FEA program) then
many additional modelling possibilities are open because the macro strain gra-
dients can be supplied to the CA arrays. The algorithms for strain redistribution
across CA cells (Step 3, page 49, Step 7, page 51) and for calculating the FE
stress (Step 11, page 52) can be modied to take this into account. This will
allow for better simulation of local strain and stress concentration ahead of a
crack.
The present CAFE models were primarily designed to model fracture propa-
gation. Accordingly the size of the smallest modelling entity (CA cell) is chosen
equal to that of the fracture propagation step. However, the model formalism
allows for much smaller entities to be used. This can be utilised e.g. for explicit
simulation of crack initiation. In this case the smallest modelling entity has to
be chosen on the basis of a crack initiation site size, typical of the simulated
5.1. UNRESOLVED PROBLEMS AND FUTURE WORK 115
material. For steels these are most usually carbides. Accordingly a much more
thorough modelling of the fracture initiation can be achieved if the CAFE model
has cells of a typical carbide size. This can be technically done in two ways.
The rst way is to create all brittle CA cells with size equal to that of a
typical carbide in a simulated steel. This approach was explored by Das et al.
(2001); Das (2002); Das et al. (2003). The results presented by Das (2002) show
very good correlation between the prediction of oxide cracking during hot rolling
and the experimental observations of the quality of the slab surface after rolling.
The second approach is to create a third layer of CA only for aliveC brittle
cells. This additional layer can be used for detailed modelling of stress elds at
a grain boundary carbide. The second approach seems more favourable because
the third level CA arrays are created only around the grain boundary carbides.
So the highest level of detalisation in the CAFE model is used only where it
is really required. Moreover, this approach allows for three CA arrays with
independent cell sizes.
A detailed simulation of fracture initialisation according to the second ap-
proach described above is recommended for future work.
5.1 Unresolved problems and future work
1. Fracture cannot propagate across the nite element boundary at present.
The information given by the Abaqus solver to the VUMAT subroutine is lim-
ited. For instance no nite element number, either external, given by the user, or
internal, used in the solver stiness matrix, is given to VUMAT. This means that it
is not possible to establish which nite elements are being processed in this call
to the subroutine. Consequently it is not possible to nd the neighbouring nite
elements from VUMAT. Thus fracture must reinitiate in each nite element in the
fracture propagation path. Fracture initiation requires more energy than prop-
agation. Therefore this limitation is likely to result in over-estimated energy
values.
One way of obtaining the numbers of neighbouring elements is to use material
point coordinates which are given to VUMAT. However, this procedure is non-
trivial and computationally expensive as the full deformation history must be
traced for each nite element in the damage zone. Work in this direction is
recommended for the future.
2. Shear ductile fracture cannot be modelled at present. The Rousselier contin-
uous ductile damage model, used in this work, can only account for volumetric
void growth.
An additional criterion has to be applied to estimate the onset of shear
instability of a ductile damage cell. There are a number of works addressing
this issue (Rice, 1977; Rousselier and Barbier, 1997; Barbier et al., 1998). The
possibility of including an appropriate shear localisation model into the present
CAFE structure was not explored due to time constraints. However this is
possible and is recommended for future work.
3. There are many model parameters which require proper tuning before good
correlation with the experiment is achieved.
Parameters such as X
max
(D)
, X
max
(B)
, c
D
, c
B
strongly aect model performance.
Yet, in this work the values for these parameters were based on a rough guess
and on some data tting. A more detailed understanding of the fracture process
at the micro-scale might help to nd metallurgically meaningful values.
It would be interesting to investigate in greater detail the inuence of the
CAFE model parameters on various simulation results, e.g. the eect of
F
on
the transition temperature range, the link between c
D
and c
B
and the upper
and the lower Charpy energy values, the inuence of the grain size distribution
on the transition temperature range. It might be better to use a dierent idea
of how X
max
(D)
and X
max
(B)
should be chosen.
This analysis was not performed as part of the present work because of
signicant demands for the computing power and time constraints.
5.2 Overall conclusions
The CAFE model proposed in this work is designed for fast three-dimensional
multi-scale analysis. One typical application of such a model is a simulation
of transitional ductile-brittle fracture in steels. This engineering problem is of
high practical importance. However, so far it could not be solved successfully
5.2. OVERALL CONCLUSIONS 117
by the pure nite element methods. The results presented in this work show
how this problem can be solved with the CAFE model.
The author believes that the model can be useful in other areas of engi-
neering. It might need to be modied, additional parameters might need to be
included, but the basic structure can remain as described in Chapter 3.
Appendix A
The CA cell neighbourhood
Figure A.1 illustrates the 26-cell neighbourhood used in both the full and the
simplied CAFE models. The neighbourhood is shown as three sections of a
3 3 3 cell cube along direction k (or 3).
Each neighbourhood cell has three characteristics. These are:
1. Direction cosines, i.e. the cosines of angles that a line, connecting the
centres of each cell with the central cell, makes with the basis axes (i, j,
k),
2. Cell coordinates relative to the coordinates of the central cell,
3. The cell number.
These characteristics are written in three lines inside each neighbourhood
cell in Figure A.1.
Because of the central symmetry only 13 out of 26 neighbouring cells have
unique combination of direction cosines. These cells have unique numbers. Their
symmetrical partner cells have the same number but with a bar on top, e.g. cell
6 is a symmtrical partner of cell 6.

Coordinates of each neighbourhood cell are, of course, unique.
119
120 APPENDIX A. THE CA CELL NEIGHBOURHOOD
j, 2
1
3
;
1
3
;
1
3
0;
1
2
;
1
2
1
3
;
1
3
;
1
3
i-1; j+1; k-1 i; j+1; k-1 i+1; j+1; k-1
3

10 7
1
2
; 0;
1
2
0; 0; 1
1
2
; 0;
1
2
i-1; j; k-1 i; j; k-1 i+1; j; k-1
6

11 4
1
3
;
1
3
;
1
3
0;
1
2
;
1
2

1
3
;
1
3
;
1
3
i-1; j-1; k-1 i; j-1; k-1 i+1; j-1; k-1
9

13 1 i, 1
2
;
1
2
; 0 0; 1; 0
1
2
;
1
2
; 0
i-1; j+1; k i; j+1; k i+1; j+1; k
2 12 8
1; 0; 0 1; 0; 0
i-1; j; k i; j; k i+1; j; k
5 central cell 5
1
2
;
1
2
; 0 0; 1; 0
1
2
;
1
2
; 0
i-1; j-1; k i; j-1; k i+1; j-1; k
8

12 2
3
;
1
3
;
1
3
0;
1
2
;
1
2
1
3
;
1
3
;
1
3
i-1; j+1; k+1 i; j+1; k+1 i+1; j+1; k+1
1 13 9
1
2
; 0;
1
2
0; 0; 1
1
2
; 0;
1
2
i-1; j; k+1 i; j; k+1 i+1; j; k+1
4 11 6
1
3
;
1
3
;
1
3
0;
1
2
;
1
2
1
3
;
1
3
;
1
3
i-1; j-1; k+1 i; j-1; k+1 i+1; j-1; k+1
7 10 3
k, 3
Figure A.1: Neighbouring cell numbering convention
Appendix B
Rousselier model
integration
Integration of the Rousselier continuous material damage model (Rousselier
et al. (1989), section 2.1.6) for a single integration point is shown below.
Where explicit time is not given, t
i+1
is assumed.
If we substitute dierentials by nite dierences the complete set of equia-
tions will have the form (Aravas, 1987; Rousselier et al., 1989; HKS, 2001):
p
m
p
eq
B()
3
1
Dexp
_
1
_
= 0 (B.1)
eq
H
_
p
eq
_
+B() Dexp
_
1
_
= 0 (B.2)
m
=
e
m
3K
p
m
(B.3)
eq
=
e
eq
3G
p
eq
(B.4)
=
p
eq
Dexp
_
1
_
(B.5)
() =
1
1 f
0
+f
0
exp
(B.6)
B() =

1
f
0
exp
1 f
0
+f
0
exp
(B.7)
where
e
m
=
1
3
e
ii
;
e
eq
=
_
3
2
S
e
ij
S
e
ij
; S
e
ij
=
e
ij

e
m
ij
;
e
ij
= E
ijkl

kl
and

kl
=
e
ij
(t
i
) +
ij
(B.8)
121
122 APPENDIX B. ROUSSELIER MODEL INTEGRATION
p
eq
=
p
eq
(t
i
) +
p
eq
(B.9)
= (t
i
) + (B.10)
For the rest of the variables refer to the Nomenclature, page 10, and to
section 2.1.6, page 23.
The equations (B.1) (B.7) are solved by the Newtons method. The strain
increments
p
m
and
p
eq
are primarily unknowns. We nd them by solving
the equations (B.1) and (B.2). If we write the equations (B.1) and (B.2) as:
_
_
_
f
_
p
m
,
p
eq
,
_
= 0
g
_
p
m
,
p
eq
,
_
= 0
(B.11)
then the solution is found by the iterative process. Each cycle the following
matrix equation is solved:
J c = A (B.12)
where
J =
_
_
f
p
m
f
p
eq
g
p
m
g
p
eq
_
_
, c =
_
_
c
m
c
eq
_
_
, A =
_
_
f
g
_
_
.
Then the strain increments are updated:
p
m

p
m
+c
m
(B.13)
p
eq

p
eq
+c
eq
(B.14)
The steps (B.12), (B.13) and (B.14) are repeated until the correction values c
m
amd c
eq
are less that the specied tolerance.
Components of Jacobian
f
p
m
= 1
D
p
eq
3
1
p
m
_
B() exp
_
1
__
(B.15)
f
p
eq
=
B()
3
1
Dexp
_
1
_
D
p
eq
3
1
p
eq
_
B() exp
_
1
__
(B.16)
g
p
m
=

p
m
_
eq
_
+D

p
m
_
B() exp
_
1
__
(B.17)
g
p
eq
=

p
eq
_
eq
H
_
p
eq
_
p
eq
+D

p
eq
_
B() exp
_
1
__
(B.18)
APPENDIX B. ROUSSELIER MODEL INTEGRATION 123
further:
p
m
_
B() exp
_
1
__
=
B()
p
m
exp
_
1
_
+
B() exp
_
1
_
1
1
_

m
p
m
1
+
m
p
m
_
(B.19)
p
eq
_
B() exp
_
1
__
=
B()
p
eq
exp
_
1
_
+
B() exp
_
1
_
1
1
_

m
epsilon
p
eq
1
+
m
p
eq
_
(B.20)
p
m
_
eq
_
=

eq
p
m
1
+
eq
p
m
(B.21)
p
eq
_
eq
_
=

eq
epsilon
p
eq
1
+
eq
p
eq
(B.22)
From (B.7):
B()
=
1
f
0
exp(1 f
0
+f
0
exp) f
0
(exp)
2
(1 f
0
+f
0
exp)
2
=
1
f
0
exp(1 f
0
)
(1 f
0
+f
0
exp)
2
(B.23)
From (B.6):
= f
0
exp (B.24)
From (B.3):
p
m
= 3K (B.25)
p
eq
= 0 (B.26)
From (B.4):
eq
p
m
= 0 (B.27)
eq
p
eq
= 3G (B.28)
Calculations of partial derivatives

p
m
and

p
eq
is slightly more complicated
because dependence between the variables ,
p
m
and
p
eq
is described by an
implicit function dened by equation (B.5).
If we regroup equation (B.5) to the following form
h
_
,
p
eq
,
p
m
_
=
p
eq
Dexp
_
1
_
= 0 (B.29)
we can get the sought partial derivatives using the formula for partial derivatives
of the implicit function (e.g. Brand, 1955):
p
m
=
h
p
m
h
(B.30)
p
eq
=
h
p
eq
h
(B.31)
In the above we used the fact that =
t
+ that leads to

p
m
=

p
m
and
p
eq
=

p
eq
.
Thus we can obtain from (B.29):
h
p
m
=
p
eq
Dexp
_
1
_
1
p
m
(B.32)
h
p
eq
= Dexp
_
1
_
p
eq
Dexp
_
1
_
1
p
eq
(B.33)
h
= 1
p
eq
Dexp
_
1
_

m
(B.34)
By sibstituting (B.32), (B.33) and (B.34) into (B.30) and (B.31), and taking
into account that
m
p
eq
= 0 and

1
=

1
, we get:
p
m
=
p
eq
Dexp
_
m
1
_
1
1
m
p
m
1
p
eq
Dexp
_
m
1
_
m
1
(B.35)
p
eq
=
Dexp
_
m
1
_
1
p
eq
Dexp
_
m
1
_
m
1
(B.36)
By simlifying these equations we get:
p
m
=
p
eq
Dexp
_
m
1
_
1
1
m
p
m
1
p
eq
Dexp
_
m
1
_
m
1
(B.37)
p
eq
=
Dexp
_
m
1
_
1
p
eq
Dexp
_
m
1
_
m
1
(B.38)
By substituting (B.24)-(B.28) into (B.19)-( B.22), (B.37) and (B.38) we get:
p
m
_
B() exp
_
1
__
= exp
_
1
__
B()
p
m
+
B()
1
_
3K
+
m
f
0
exp

p
m
__
(B.39)
APPENDIX B. ROUSSELIER MODEL INTEGRATION 125
p
eq
_
B() exp
_
1
__
= exp
_
1
__
B()
p
eq
+
B()
m
f
0
exp

p
eq
_
(B.40)
p
m
_
eq
_
=
eq
f
0
exp

p
m
(B.41)
p
eq
_
eq
_
=
3G
+
eq
f
0
exp

p
eq
(B.42)
p
m
=
3K
p
eq
Dexp
_
m
1
_
p
eq
Dexp
_
m
1
_
m
f
0
exp
_ (B.43)
p
eq
=
Dexp
_
m
1
_
1
p
eq
Dexp
_
m
1
_
m
1
f
0
exp
(B.44)
The components of Jacobian can be found as:
f
p
m
= 1
D
p
eq
3
1
p
m
_
B() exp
_
1
__
(B.45)
f
p
eq
=
D
3
1
_
B() exp
_
1
_
+
p
eq
p
eq
_
B() exp
_
1
___
(B.46)
g
p
m
=
eq
f
0
exp

p
m
+D

p
m
_
B() exp
_
1
__
(B.47)
g
p
eq
=
3G
+
eq
f
0
exp

p
eq
H
_
p
eq
_
p
eq
+
D

p
eq
_
B() exp
_
1
__
(B.48)
where
p
m
_
B() exp
_
1
__
= exp
_
1
__
B()
p
m
+
B()
1
_
3K
+
m
f
0
exp

p
m
__
(B.49)
p
eq
_
B() exp
_
1
__
= exp
_
1
__
B()
p
eq
+
B()
m
f
0
exp

p
eq
_
(B.50)
p
m
=
3K
p
eq
Dexp
_
m
1
_
p
eq
Dexp
_
m
1
_
m
f
0
exp
_ (B.51)
p
eq
=
Dexp
_
m
1
_
1
p
eq
Dexp
_
m
1
_
m
1
f
0
exp
(B.52)
B()
=

1
f
0
exp(1 f
0
)
(1 f
0
+f
0
exp)
2
(B.53)
Equations (B.1) (B.14) and (B.45) (B.53) are required to nd
p
m
and
p
eq
.
When
p
m
and
p
eq
are known then
m
and
eq
are found from equations
(B.3) and (B.4), from (B.10) and (B.5) and
ij
=

eq
e
eq
S
e
ij
+
m
ij
(B.54)
p
ij
=
3
2
S
e
ij
e
eq
p
eq
+
p
m
ij
, (B.55)
e
ij
=
e
ij
(t
i
) +
ij
p
ij
(B.56)
As seen from equations (B.8) (B.10) the elastic strain tensor,
e
ij
, equivalent
plastic strain,
p
eq
, and the damage variable, , have to be stored from one time
increment to another throughout the analysis.
Bibliography
Andersson, H. (1977), Analysis of a model for void growth and coalescence
ahead of a moving crack tip, Journal of the Mechanics and Physics of Solids
25, 217233.
Andrews, R. M., I. C. Howard, A. Shterenlikht and J. R. Yates (2002), The
Eective Resistance of Pipeline Steels to Running Ductile Fractures; Model-
ling of Laboratory Test Data, in A.Neimitz, I. V.Rokach, D.Koca nda and
K.Go los, eds, ECF14, Fracture Mechanics Beyond 2000, EMAS Publica-
tions, Sheeld, UK, pp. 6572.
Aravas, N. (1987), On the numerical integration of a class of pressure-dependent
plasticity models, International journal for numerical methods in engineering
24, 1395416.
Argon, A. S. and J. Im (1975), Separation of second phase particles in
spheroidized 1045 steel, Cu - 0.6% Cr alloy, and maraging steel in plastic
straining, Metallurgical Transactions A 6, 83951.
Argon, A. S., J. Im and R. Safoglu (1975), Cavity fomation from inclusions in
ductile fracture, Metallurgical Transactions A 6, 825837.
Ashby, M. F., F. J. Blunt and M. Bannister (1989), Flow characteristics of
highly constrained metal wires, Acta Metallurgica 37(7), 18471857.
Averbach, B. L., D. K. Felbeck, G. T. Hahn and D. A. Thomas, eds (1959),
Fracture. Proceedings of an International Conference on the Atomic Mecha-
nisms of Fracture held in Swampscott, Massachusetts, April 1216, 1959, The
127
128 BIBLIOGRAPHY
Technology Press of Massachusetts Institute of Technology and John Wiley
& Sons, New York.
Barbier, G., G. Rousselier, R. Lebrun and Y. Sun (1998), Analytical predic-
tion of intrinsic length scales in strain localisation, Journal de Physique IV
8(P8), 1927.
Bates, R. C. (1984), Modelling of ductile fracture by microvoid coalescence for
the prediction of fracture toughness, in J. M.Wells and J. D.Landes, eds,
Fracture: Interactions of Microstructure, Mechanisms and Mechanics, The
Metallurgical Society of AIME, pp. 117155.
Batisse, R., M. Bethmont, G. Devesa and G. Rousselier (1987), Ductile fracture
of a 508 Cl 3 steel in relation with inclusion content: the benet of the local
approach to fracture and continuum damage mechanics, Nuclear Engineering
and Design 105(1), 113120.
Beachem, C. D. (1963), An electron fractography study of the inuence of plas-
tic strain conditions upon ductile rupture processes in metals, Transactions
of the ASM 56, 318326.
Beachem, C. D. (1975), The eects of crack tip plastic ow directions upon
microscopic dimple shapes, Metallurgical Transactions A 6, 377383.
Belytschko, T., W. K. Liu and B. Moran (2000), Nonlinear Finite Elements for
Continua and Structures, John Wiley & Sons.
Beremin, F. M. (1983), A local criterion for cleavage fracture of a nuclear
pressure vessel steel, Metallurgical Transactions A 14, 22772287.
Berg, C. A. (1970), Plastic dilation and void interaction, in M. F.Kanninen, ed.,
Inelastic Behaviour of Solids, McGraw-Hill, pp. 171209.
Beynon, J. H., S. Das, I. C. Howard, E. J. Palmiere and A. Chterenlikht
(2002), The combination of Cellular Automata and Finite Elements for the
Study of Fracture; the CAFE Model of Fracture, in A.Neimitz, I. V.Rokach,
D.Koca nda and K.Go los, eds, ECF14, Fracture Mechanics Beyond 2000,
EMAS Publications, Sheeld, UK, pp. 241248.
BIBLIOGRAPHY 129
Bhattacharjee, D. and C. L. Davis (2002), Inuence of processing history on
mesotexture and microstructure-toughness relationship in control-rolled and
normalised steels, Scripta Materialia 47(12), 825831.
Bhattacharjee, D., C. L. Davis and J. F. Knott (2003), Predictability of charpy
impact toughness in thermomechanically control rolled (TMCR) microalloyed
steels, Ironmaking and Steelmaking 30(3), 249255.
Bilby, B. A., I. C. Howard, A. M. Othman, D. P. G. Lidbury and A. H. Sherry
(1994), Prediction of spinning cylinder tests 2 and 3 using continuum dam-
age mechanics, in H. S.Mehta, ed., Fracture Mechanics Applications, PVP-
Vol.287/MD-Vol.47, The American Society of Mechanical Engineeers, pp. 3
10.
Bilby, B. A., I. C. Howard and Z. H. Li (1994), Mesh independent cell models
for continuum damage theory, Fatigue and Fracture of Engineering Materials
and Structures 17(10), 12211233.
Bluhm, J. I. and R. J. Morrissey (1964), Preliminary investigation of the frature
mechanism in a tensile specimen, US Army Materials Research Agency, Wa-
tertown, Massachusetts, USA.
Bluhm, J. I. and R. J. Morrissey (1966), Fracture in a tensile specimen, in
Proceedings of First International Conference on Fracture, Vol. 3, Japanese
Society for Strength and Fracture of Materials, Sendai, Japan, pp. 17391780.
Brand, L. (1955), Advanced calculus, John Wiley & Sons, New York.
Bridgman, P. W. (1952), Studies in Large Plastic Flow and Fracture, Metallurgy
and Metallurgical Engineering Series, 1 edn, McGraw-Hill.
Brown, L. M. and J. D. Embury (1973), The initiation and growth of voids at
second-phase particles, in Proceedings of the 3rd International Conference
on Strength of Metals and Alloys, Institute of Metals, London, pp. 164169.
BS EN 10045-1 (1990), British Standard. Metallic materials Charpy impact
test Part 1: Test method, British Standards Institution.
130 BIBLIOGRAPHY
Burstow, M. C. (1998), Local approach analysis of fracture in the transition
region, Technical report, Department of Mechanical Engineering, The Uni-
versity of Sheeld.
Burstow, M. C. (2003), A re-assessment of parameters tuning for the beremin
model using the toughness scaling technique, to be submitted.
Burstow, M. C., D. W. Beardsmore, I. C. Howard and D. P. G. Lidbury (2003),
The prediction of constraint-dependent R6 failure assessment lines for a pres-
sure vessel steel via micromechanical modelling of fracture, to be submitted.
CEI (2002), EnSight User Manual, Version 7.4, Computational Engineerig In-
ternational, Inc., Apex, North Carolina, USA.
Compaq (1999), Visual Fortran Version 6.1, Digital Equipment Corporation,
USA.
Cottrell, A. H. (1959), Theoretical aspects of fracture, in Averbach et al. (1959),
pp. 2053.
Cottrell, A. H. (1967), The nature of metals, in D.Flannagan, F.Bello and
P.Morrison, eds, Meterials, W. H. Freeman and Company, San Francisco,
pp. 3955.
Cox, T. B. and J. R. Low, Jr (1974), An investigation of the plastic fracture of
AISI 4340 and 18 Nickel 200 grade maraging steels, Metallurgical Transac-
tions 5, 14571470.
Curry, D. A. and J. F. Knott (1978), Eects of microstructure on cleavage
fracture stress in steel, Metal Science 12, 511514.
Czochralski, J. (1924), Moderne Metallkunde in Theorie und Praxis (Modern
metal science in theory and practice), Springer-Verlag, Berlin.
Das, S. (2002), The eect of boundary conditions and material data repre-
sentation on the simulation of deformation during hot rolling, PhD thesis,
Department of Engineering Materials, The University of Sheeld, UK.
Das, S., E. J. Palmiere and I. C. Howard (2001), CAFE: A new approach to
the modelling of multipass hot rolling, in Proceedings of Modelling of Metal
BIBLIOGRAPHY 131
Rolling Processes Symposium 11 - Through Process Modelling, The Institute
of Materials, London, pp. 3340.
Das, S., E. J. Palmire and I. C. Howard (2003), CAFE: a tool for modelling
thermomechanical processes, in E. J.Palmire, M.Mahfouf and C.Pinna, eds,
Thermomechanical Processing: Mechanics, Microstructure and Control, Pro-
ceeding of the International Conference, Sheeld, 23 - 26 June, 2002, IMM-
PETUS, BBR Solutions, Chestereld, UK, pp. 296301.
Davis, C. L. (2003), Private communication. Metallurgy and Materials, The
University of Birmingham.
Devroye, L. (1987), A Course in Density Estimation, Vol. 14 of Progress in
Probability and Statistics, Birkh auser, Boston, USA.
Drucker, D. C. (1951), A more fundamental approach to stress-strain relations,
in Proceedings of the First U.S. National Congress of Applied Mechanics,
pp. 487491.
Drucker, D. C. (1959), A denition of stable inelastic material, Journal of
Applied Mechanics 26, 101106.
Eripret, C., D. P. G. Lidbury, A. Sherry and I. C. Howard (1996), Prediction
of fracture in the transition regime: application to an A533B pressure vessel
steel, Journal de Physique IV 6(C6), 315323.
Faleskog, J. and C. F. Shih (1997), Micromechanics of coalescenceI. Synergetic
eects of elasticity, plastic yelding and multi-size-scale voids, Journal of the
Mechanics and Physics of Solids 45(1), 2150.
Faleskog, J., X. Gao and C. F. Shih (1998), Cell model for nonlinear fracture
analysis I. micromechanics calibration, International Journal of Fracture
89(4), 355373.
Folch, L. C. A. (1997), Application of the local damage mechanics approach to
transition temperature behaviour in steels, PhD thesis, University of Manch-
ester Institute of Science and Technology, Department of Civil and Structural
Engineering.
132 BIBLIOGRAPHY
Folch, L. C. A. and F. M. Burdekin (1999), Application of coupled brittle-
ductile model to study correlation between charpy energy and fracture tough-
ness values, Engineering Fracture Mechanics 63(1), 5780.
Franklin, A. G. (1969), Comparison between a quantitative microscopic and
chemical methods for assessment of non-metallic inclusions, Journal of the
Iron and Steel Institute 207, 181186.
Gandin, C. A., J. L. Desbiolles, M. Rappaz and P. Thevoz (1999), A three-
dimensional cellular automaton nite element model for the prediction of
solidication grain structures, Metallurgical and Materials Transactions A
30(12), 31533165.
Gao, X., J. Faleskog and C. F. Shih (1998), Cell model for nonlinear frac-
ture analysis II. fracture process calibration and verication, International
Journal of Fracture 89(4), 375398.
Gao, X., J. Faleskog and C. F. Shih (1999), Analysis of ductile to cleavage
transition in part-through cracks using a cell model incorporating statistics,
Fatigue and Fracture of Engineering Materials and Structures 22(3), 239250.
Germain, P., Q. S. Nguyen and P. Suquet (1983), Continuum thermodynamics,
Journal of Applied Mechanics 50(12), 10101020.
Gilman, J. J. (1969), Micromechanics of Flow in Solids, Materials Science and
Engineering, McGraw-Hill, New York.
Gladman, T. (1997), The Physical Metallurgy of Microalloyed Steels, Iron and
Steel Institute, London.
Gladman, T., B. Holmes and F. B. Pickering (1970), Work hardening of low-
carbon steels, Journal of the Iron and Steel Institute 208, 172183.
Gladman, T., B. Holmes and I. D. McIvor (1971), Eects of Second Phase
Particles on the Mechanical Properties of Steels, Iron and Steel Institute,
London.
Gologanu, M., J. B. Leblond and J. Devaux (1993), Approximate models for
ductile metals containing non-spherical voids case of axisymmetric pro-
BIBLIOGRAPHY 133
late ellipsiodal cavities, Journal of the Mechanics and Physics of Solids
41(11), 17231754.
Gologanu, M., J. B. Leblond and J. Devaux (1994), Approximate models
for ductile metals containing non-spherical voids case of axisymmetric
oblate ellipsiodal cavities, Journal of Engineering Materials and Technology
116(3), 290297.
Goods, S. H. and L. M. Brown (1979), The nucleation of cavities by plastic
deformation, Acta Metallurgica 27, 115.
Grith, A. A. (1921), The phenomena of fracture and ow in solids, Philo-
sophical Transactions of the Royal Society of London, Series A 221, 16398.
Grith, A. A. (1924), The theory of rupture, in Proceedings of the rst inter-
national congress for applied mechanics, Delft, pp. 5563.
G ucer, D. E. and J. Gurland (1962), Comparison of the statistics of two fracture
modes, Journal of the Mechanics and Physics of Solids 10(4), 365373.
Gurland, J. (1972), Observations of the fracture of cementite particles in a
spheroidized 1.05%C steel deformed at room temperature, Acta Metallurgica
20, 735741.
Gurland, J. and J. Plateau (1963), The mechanism of ductile rupture of metals
containing inclusions, Transactions of the ASM 56, 443454.
Gurson, A. L. (1977a), Continuum theory of ductile rupture by void nucleation
and growth: Part I - yield criteia and ow rules for porous ductile media,
Journal of Engineering Materials and Technology 99, 215.
Gurson, A. L. (1977b), Porous rigid-plastic materials containing rigid inclusions
- Yield function, plastic potential, and void nucleation, in D. M. R.Taplin, ed.,
Proceedings of the international conference on fracture, Vol. 2A, Pergamon
Press, pp. 357364.
Hahn, G. T. (1984), The inuence of microstructure on brittle fracture tough-
ness, Metallurgical Transactions A 15(6), 947959.
134 BIBLIOGRAPHY
Hahn, G. T., B. L. Averbach, W. S. Owen and M. Cohen (1959), Initiation
of cleavage microcracks in polycrystalline iron and steel, in Averbach et al.
(1959), pp. 91116.
Hayden, H. W. and S. Floreen (1969), Observations of localised deformation
during ductile fracture, Acta Metallurgica 17, 213214.
Hertzberg, R. W. (1996), Deformation and Fracture Mechanics of Engineering
Materials, 4 edn, John Wiley & Sons.
Hesselbarth, H. W. and I. R. G obel (1991), Simulation of recrystallization by
cellular automata, Acta Metallurgica et Materialia 39(9), 21352143.
Hill, R. (1983), The Mathematical Theory of Plasticity, The Oxford Engineering
Science Series, Oxford: Clarendon Press.
HKS (2001), ABAQUS Theory Manual, Version 6.2.
Howard, I. C., Z. H. Li and M. A. Sheikh (2000), Modeling the ductile to cleavage
transition in steels and structures, in P. C.Paris and K. L.Jerina, eds, Fatigue
and Fracture Mechanics: 30th Volume, ASTM STP 1360, American Society
for Testing and Materials, West Conshohocken, PA, USA, pp. 152168.
Howard, I. C., Z. H. Li, M. A. Sheikh, D. P. G. Lidbury and A. H. Sherry (1996),
The simulation of the fourth spinning cylinder test using damage mechanics,
in K.Yoon, ed., Fatigue and Fracture, vol. 2, Proc. ASME PVP Conference,
Montreal, July 21-26., The American Society of Mechanical Engineeers, New
York, pp. 257265.
Huang, Y., J. W. Hutchnison and V. Tvergaard (1991), Cavitation instabili-
ties in elastic-plastic solids, Journal of the Mechanics and Physics of Solids
39(2), 223241.
Jang, J. I., B. W. Lee, J. B. Ju, D. Kwon and W. S. Kim (2003), Experimental
analysis of the practical LBZ eects on the brittle fracture performance of
cryogenic steel HAZs with respect to crack arrest toughness near fusion line,
Engineering Fracture Mechanics 70(10), 12451257.
Johnson, C. (1987), Numerical Solution of Partial Dierential Equations by the
Finite Element Method, Cambridge University Press.
BIBLIOGRAPHY 135
Kachanov, L. M. (1971), Foundations of the Theory of Plasticity, North-Holland
Publishing Company, Amsterdam, Netherlands.
Kelly, A. and G. W. Groves (1970), Crystallography and Crystal Defects, Long-
man, London, UK.
Khalili, A. and K. Kromp (1991), Statistical properties of Weibull estimators,
Journal of Materials Science 26(24), 67416752.
Khan, A. S. and S. Huang (1995), Continuum Theory of Plasticity, John Wiley
& Sons, New York.
Knott, J. F. (1973), Fundamentals of Fracture Mechanics, Butterworth, London.
Kohout, J. (2001), Various regressional functions tting transition curves, in
D.Fran cois and A.Pineau, eds, Proceedings of the Charpy Centenary Confer-
ence, Poitiers, 2 5 October 2001., ESIS, Societe Fran caise de Metallurgie
et de Materiaux, France, pp. 8188.
Koppenhoefer, C. and R. H. Dodds, Jr (1998), Ductile crack growth in pre-
cracked CVN specimens: numerical studies, Nuclear Engineering and Design
180(3), 185289.
Kroon, M. and J. Faleskog (2002), A probabilistic model for cleavage fracture
with a length scale-inuence of material parameters and constraints, Inter-
national Journal of Fracture 118(2), 99118.
Lemaitre, J. (1985), A continuum damage mechanics model for ductile fracture,
Journal of Engineering Materials and Technology 107(1), 8389.
Lemaitre, J. (1996), A Course on Damage Mechanics, 2 edn, Springer.
Lidbury, D. P. G., A. H. Sherry, B. A. Bilby, I. C. Howard, Z. H. Li and C.
Eripret (1994), Prediction of the rst spinning cylinder test using continuum
damage mechanics, Nuclear Engineering and Design 152(1-3), 110.
Lin, T., A. G. Evans and R. O. Ritchie (1987), Stochastic modeling of the
independent roles of particle size and grain size in transgranular cleavage
fracture, Metallurgical Transactions A 18(4), 641651.
136 BIBLIOGRAPHY
Liu, C. T. and J. Gurland (1968), The fracture behaviour of spheroidized carbon
steels, Transactions of the ASM 61, 156167.
Malik, L., L.N. Pussegoda, B.A. Graville and W.R. Tyson (1996), Crack ar-
rest toughness of a heat-aected zone containing local brittle zones, Journal
of Oshore Mechanics and Arctic Engineering, Transactions of The ASME
118(4), 292299.
McClintock, F. A. (1968), A criterion for ductile fracture by the growth of
voids, Journal of Applied Mechanics 35, 363371.
McDowell, D. L., ed. (1997), Applications of Continuum Damage Mechanics to
Fatigue and Fracture, American Society for Testing and Materials, STP 1315.
McLachlan, G. J. and K. E. Basford (1988), Mixture Models: Inference and
Applications to Clustering, Vol. 84 of Statistics: Textbooks ad Monographs,
Marcel Dekker, Inc., 270 Madison Avenue, New York, USA.
McMahon, Jr, C. J. and M. Cohen (1965), Initiation of cleavage in polycrys-
talline iron, Acta Metallurgica 13, 591604.
Nadai, A. (1950), Theory of Flow and Fracture of Solids, Engineering Societies
Monographs, 2 edn, McGraw-Hill, New York.
Needleman, A. (1972), Void growth in an elastic-plastic medium, Journal of
Applied Mechanics 39, 964970.
Nohava, J., P. Hausild, M. Karlk and P. Bompard (2002), Electron backscat-
tering diraction analysis of secondary cleavage cracks in a reactor pressure
vessel steel, Materials Characterization 49(3), 211217.
Pardoen, T. and J. W. Hutchinson (2000), An extended model for void
growth and coalescence, Journal of the Mechanics and Physics of Solids
48(12), 24672512.
Prager, W. (1959), An Introduction to Plasticity, The Addison-Wesley series in
the engineering sciences, Addison-Wesley, Reading, Massachusetts, USA.
Puttick, K. E. (1959), Ductile fracture in metals, Philosophical Magazine, Se-
ries 8 4(44), 964969.
BIBLIOGRAPHY 137
Raabe, D. and R. C. Becker (2000), Coupling of a crystal plasticity nite-
element model with a probabilistic cellular automaton for simulating primary
static recrystallization in aluminium, Modelling and Simulation in Materials
Science and Engineering 8(4), 445462.
Rabotnov, Yu. N. (1969), Creep Problems in Structural Members. Translated
from Russian, originally entitled: Polzuchest Elementov Konstruktsii, Nauka,
Moscow, 1966, North-Holland.
Rhines, W. J. (1961), Ductile fracture by the growth of pores, Masters thesis,
MIT, Department of Mechanical Engineering.
Rice, J. R. (1977), The localization of plastic deformation, in W. T.Koiter, ed.,
Theoretical and applied mechanics, North-Holland, Amsterdam, pp. 207
220.
Rice, J. R. and D. M. Tracey (1969), On the ductile enlargement of voids in
triaxial stress elds, Journal of the Mechanics and Physics of Solids 17, 201
217.
Ritchie, R. O., J. F. Knott and J. R. Rice (1973), On the relationship between
critical tensile stress and fracture toughness in mild steel, Journal of the
Mechanics and Physics of Solids 21, 395410.
Rogers, H. C. (1960), Tensile fracture of ductile metals, Transactions of the
Metallurgical Society of AIME 218, 498506.
Rousselier, G. (1981), Finite deformation constitutive relations including ductile
fracture damage, in S.Nemat-Nasser, ed., Three Dimensional Constitutive
Relations and Ductile Fracture, North-Holland, pp. 331355.
Rousselier, G. (1987), Ductile fracture models and their potential in local ap-
proach to fracture, Nuclear Engineering and Design 105(1), 97111.
Rousselier, G. and G. Barbier (1997), Analytical modeling of the shear mode
and opening mode of ductile fracture, in A. B.Geltmacher, P.Matic and
K.Sadananda, eds, Recent Advances in Fracture. Proceedings of the TMS
meeting, 1013 February 1997, Orlando, Florida, The Minerals, Metals &
Materials Society, USA, pp. 195203.
138 BIBLIOGRAPHY
Rousselier, G., J.-C. Devaux, G. Mottel and G. Devesa (1989), A methodology
of ductile fracture analysis based on damage mechanics: an illustration of a
local approach of fracture, in J. D.Landes, A.Saxena and J. G.Merkle, eds,
Nonlinear fracture mechanics: Volume II - Elastic-Plastic Fracture, ASTM
STP 995, Philadelphia, pp. 332354.
Ruggieri, C. (1998), Probabilistic treatment of fracture using two failure mod-
els, Probabilistic Engineering Mechanics 13(4), 309319.
Ruggieri, C., F. Minami and M. Toyoda (1995), A statistical approach for
fracture of brittle materials based on the chain-of-bundles model, Journal of
Applied Mechanics 62(2), 320328.
Scott, D. W. (1992), Multivariate Density Estimation: Theory, Practice, and
Visualization, Wiley Series in Probability and Mathematical Statistics, John
Wiley & Sons.
Sherry, A. H., D. W. Beardsmore, D. P. G. Lidbury, M. A. Sheikh and I. C.
Howard (1998), Remanent life assessment using the local approach a predic-
tion of the outcome of the NESC experiment, in Remanent Life Prediction
IMechE Seminar, London, 26 November 1997, Mechanical Engineering Pub-
lishing, Bury St Edmunds, UK, pp. 87103.
Shterenlikht, A., S. H. Hashemi, J. R. Yates, I. C. Howard and R. M. Andrews
(2003), Assessment of an instrumented charpy impact machine, Submitted
to the International Journal of Fracture, August.
Smith, E. (1966a), The formation of a cleavage crack in a crystalline solid I,
Acta Metallurgica 14, 985989.
Smith, E. (1966b), The formation of a cleavage crack in a crystalline solid II,
Acta Metallurgica 14, 991996.
Smith, G. D. (1965), Numerical Solution of Partial Dierential Equations, Ox-
ford Mathematical Handbooks, Oxford University Press.
Szczepi nski, W. (1982), On the mechanisms of ductile fracture of metals, in
G. C.Sih and H.Zorski, eds, Defects and fracture. Proceedings of First Inter-
BIBLIOGRAPHY 139
national Symposium on Defects and Fracture, held at Tuczno, Poland, October
13-17, 1980, Martinus Nijho Publishers, The Hague, pp. 155163.
Thomason, P. F. (1968), A theory for ductile fracture by internal necking of
cavities, Journal of the Institute of Metals 96, 360365.
Thomason, P. F. (1981), Ductile fracture and the stability of incompressible
plasticity in the presence of microvoids, Acta Metallurgica 29(5), 763777.
Thomason, P. F. (1982), An assessment of plastic-stability models of ductile
fracture, Acta Metallurgica 30(1), 279284.
Thomason, P. F. (1985a), A three-dimentional model for ductile fracture by the
growth and coalescence of microvoids, Acta Metallurgica 33(6), 10871095.
Thomason, P. F. (1985b), Three-dimentional models for the plastic limit-loads
at incipient failure of the intervoid matrix in ductile porous solids, Acta
Metallurgica 33(6), 10791085.
Thomason, P. F. (1990), Ductile Fracture of Metals, Pergamon Press.
Thomason, P. F. (1993), Ductile fracture by the growth and coalescence of
microvoids of non-uniform size and spacing, Acta Metallurgica et Materialia
41(7), 21272134.
Thomason, P. F. (1998), A view on ductile-fracture modelling, Fatigue and
Fracture of Engineering Materials and Structures 21(9), 11051122.
Thompson, A. W. and J. F. Knott (1993), Micromechanics of brittle fracture,
Metallurgical Transactions A 24(3), 523534.
Tipper, C. F. (1949), The fracture of metals, Metallurgia 39, 133137.
Tvergaard, V. (1981), Inuence of voids on shear band instabilities under plain
strain conditions, International Journal of Fracture 17(4), 389407.
Tvergaard, V. (1982a), Ductile fracture by cavity nucleation between larger
voids, Journal of the Mechanics and Physics of Solids 30(4), 265286.
Tvergaard, V. (1982b), On localization in ductile materials containing spherical
voids, International Journal of Fracture 18(4), 237252.
140 BIBLIOGRAPHY
Tvergaard, V. (1990), Material failure by void growth to coalescence, Advances
in Applied Mechanics 27, 83151.
Tvergaard, V. and A. Needleman (1984), Analysis of the cup-cone fracture in
a round tensile bar, Acta Metallurgica 32(1), 157169.
Vandyousse, M. and A. L. Greer (2002), Application of cellular automaton-
nite element model to the grain renement of directionally solidied Al-4.15
wt% Mg alloys, Acta Materialia 50(7), 16931705.
Von Neumann, J. (1966), Theory of Self-Reproducing Automata, Edited and
completed by A. W. Burks, University of Illinois Press, Urbana, Illinois, USA.
Wallin, K. (1991), Statistical modelling of fracture in the ductile-to-brittle tran-
sition region, in J. G.Blauel and K. H.Schwalbe, eds, Defect Assessment in
Components Fundamentals and Applications, ESIS/EGF9, Mechanical En-
gineering Publications, London, pp. 415445.
Weibull, W. (1951), A statistical distribution function of wide applicability,
Journal of Applied Mechanics 18, 293297.
Wu, S. J. and C. L. Davis (2003), Eect of duplex ferrite grain size distribution
on local fracture stresses of Nb-microalloyed steels, Presented at 13th Inter-
national Conference on Strength of Materials, Hungary, 25 30th August. To
be published in Journal of Materials Science and Engineering A.
Xia, L. and F. S. Shih (1996), Ductile crack growth III. transition to cleavage
fracture incorporating statistics, Journal of the Mechanics and Physics of
Solids 44(4), 603639.
Xia, L. and L. Cheng (1997), Transition from ductile tering to cleavage fracture:
A cell-model approach, International Journal of Fracture 87(3), 289306.
Yamamoto, H. (1978), Conditions for shear localization in the ductile fracture
of void-containing materials, International Journal of Fracture 14(4), 347
365.
Zhang, Z. L., C. Thaulow and J. degard (2000), A complete Gurson model
approach for ductile fracture, Engineering Fracture Mechanics 67(2), 155
168.
BIBLIOGRAPHY 141
Ziegler, H. (1977), An Introduction to Thermomechanics, Vol. 21 of North-
Holland series in Applied Mathematics and Mechanics, North-Holland.
Zikry, M. A. and M. Kao (1996), Inelastic microstructural failure mechanisms
in crystalline materials with high angle grain boundaries, Journal of the Me-
chanics and Physics of Solids 44(11), 17651798.

3D CAFE Modelling of Transitional Ductile - Brittle Fracture in Steels

Uploaded by

Copyright:

Available Formats

3D CAFE Modelling of Transitional Ductile - Brittle Fracture in Steels

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

3D CAFE Modelling of Transitional Ductile - Brittle Fracture in Steels

Uploaded by

Copyright:

Available Formats

Department of Mechanical Engineering

The University of Sheeld

shape parameter of Weibull dustribution

location parameter of Weibull dustribution

scale parameter of Weibull dustribution

- is the void nucleation intensity and

(f) to model the rapid loss of

(f) was chosen as:

C. For this temperature the values

= 2 and the scale parameter was taken as W

= 2 and the scale parameter was taken as W

C, three runs at each temperature.

C and very low again at 50

C and decreasing to-

C and three at 175

C and the transitional temperatures are from 150

C. The level of scatter in the transitional region is higher than in the

(for the distribution of f

C (Figures 4.17.a and

C (Figures 4.17.c and 4.17.d). However in all four cases

C (Figures 4.17.e and 4.17.f) are quite

C (Figures 4.17.g and 4.17.h), 50

C (Figures 4.17.k and 4.17.l), one can see an

C from a 145 mm thick slab and includes

C. The nal slab thickness was

= 2 and the scale parame-

within a single cleavage facet. On the other

is used in this example.

C and the corresponding Charpy energy values

C. However, at least three simulations have to be per-

, is only a rough guess

C, three runs at each temperature.

C agree very well with

C, the maximum dierence between

C. The experimental scatter can only be assessed at T = 50

and the energies there are 20J 40J.

C. However, the experimental transition temperature range

C the maximum dierence is 60%.

C (more than 5 times). However, these experimental

C. The scatter at the upper shelf is very low, less than

C. All four fracture surfaces show very little areas of

C (Figure 4.43.b) show lower values of the brittle

C (Figure 4.43.h) have 10% and 0% brittle phase accord-

C the dominant fracture mode predicted by

C than at the lower shelf temperatures,

C are quite far

C, Figures 4.43.g 4.43.l, the

C Figure 4.43.g shows more islands of ductile

C, Figures 4.43.j 4.43.l. In Figure 4.43.l there is a large brittle region in

C, Figures 4.43.m 4.43.o, the variation from one

C) = 505 MPa and

C) = 453 MPa (Figure 4.11.a) and the

6 is a symmtrical partner of cell 6.

You might also like