Petrinic (Impact)

VII International Conference on Computational Plasticity
COMPLAS VII
E. Oñate and D. R. J. Owen (Eds)
 CIMNE, Barcelona, 2003
AN APPROACH TO PREDICTIVE MODELLING OF FAILURE IN

METALS SUBJECTED TO IMPACT LOADING
Nikica Petrinic, Li Wang and Benjamin C.F. Elliott
Solid Mechanics Group
Department of Engineering Science
University of Oxford
Parks Road, Oxford OX1 3PJ, UK
e-mail: Nik.Petrinic@eng.ox.ac.uk, web page: http://www.eng.ox.ac.uk/~ftgwww
Key words: Impact engineering, high strain rate testing, constitutive modelling, finite
element method, inverse modelling.
Abstract: This paper introduces an approach to predictive modelling of failure in metals

subjected to rapidly applied loading based upon an integrated experimental-numerical
methodology. To illustrate the adopted approach several aspects of uniaxial tensile testing at
high rates of strain are presented and the use of such experimentally obtained data in
improving numerical modelling capabilities is described. Characterisation of an advanced
metallic alloy for aerospace applications is outlined together with the modifications of an
existing constitutive model that couples the temperature and path dependent elasto-visco-
plasticity and damage mechanics. The model has been adapted to incorporate the
regularisation of numerical solutions affected by the occurrence of material instability
preceding cracking and fracture. The non-measurable parameters of the constitutive model
needed to describe accurately the observed bahaviour have been identified by means of
inverse modelling.
1
1 INTRODUCTION
The successful application of metals in engineering designs follows from the ability of
metallic alloys to be fabricated to the desired shape/size and to satisfy the anticipated service
requirements. This ability is determined by materials’ chemical and physical properties. The
focus of this paper is on the identification of mechanical properties of metallic alloys which
describe the relationship between forces (surface tractions, body forces, inertial forces) acting
on the material and the material’s resistance to thus imposed deformation through to failure.
In particular, the response to impact loading, duration of which is less than the natural period
of vibration of the component or structure under consideration, is investigated.
Failure is a change in geometric and/or material properties of components or structures that
is sufficient for their function to be lost or impaired. In the case of impact loading the
inability of materials to further dissipate energy is most often considered a failure. The ability
of ductile metals to dissipate energy through inelastic deformation ceases with cracking.
Cracking to the extent that a component is separated into two or more pieces is termed
fracture. Further energy dissipation following fracture may occur through sliding contact but
this is largely dependent on component/structure’s geometry and loading conditions rather
than just upon the material under consideration.
Designers of sea/land/air/space-craft are under constant pressure to decrease the weight of
their components and structures in order to enable reductions in fuel consumption and
associated emissions12,20. In addition, it is desired that the minimum number of large-scale
and full-scale prototype tests required during the design and development not be exceeded.
Achieving these reductions whilst improving performance and safety is a major challenge
facing the automotive and aerospace industries today. The ability to predict material’s
response to given loading is key to accomplishing such objectives.
Unlike in the majority of manufacturing processes and in service, in which it is important
to prevent the materials from being loaded beyond the onset of material instability22, in the
case of impact loading resulting from collisions or explosions the materials are virtually
without exception loaded beyond that limit. It is therefore vital that the adopted predictive
numerical modelling methodologies used in impact engineering comprise constitutive models
capable of accurately representing the behaviour of materials subjected to rapidly applied
loading through to their failure characterised by cracking and fracture.
Most advanced ductile metallic alloys for automotive and aerospace applications exhibit
some form of rate dependent behaviour21,31,42. The origins of such behaviour lie in the
materials’ microstructure32. However, detailed analysis of underlying microscopic processes
requires discretisation into such small representative volume elements1,25,26 that its use
becomes prohibitive in large scale numerical simulations of real components and structures.
Consequently, a large number of physically based phenomenological constitutive models27,32
have been developed which rely upon homogenisation of relevant microscopic material
processes within representative volume elements determined by adopted discretisation
techniques. Such an approach to predictive modelling of material response to loading relies
on experimental observations and measurements for identification of parameters in the
mathematical equations devised to describe physical processes of deformation and failure at
2
macroscopic scales.
In this paper an iterative integrated experimental-numerical approach to predictive
modelling of the response of advanced metallic alloys to impact loading is presented. In this
approach (Figure 1), experimental observations and measurements at different rates of strain
are conducted first, in order to observe and classify the material behaviour. Subsequently, a
theoretical framework is chosen in order to enable numerical algorithms to be devised and
implemented within the existing computational framework(s) based upon the discrete/finite
element method3,6,7,9,11,34. In addition, the inability to identify experimentally all required
parameters is compensated by incorporation of inverse modelling techniques. This approach
enables not only the evaluation of non-measurable constitutive parameters but also the
validation of individual components of the utilised computational framework35.
Experimental observations Modify experiments
Experimental data
Theoretical abstraction Modify theory
Numerical algorithms Modify algorithms
Data for numerical analysis
Simulation ABAQUS
LS-DYNA
PAM-CRASH
DEST
Optimisation
Tolerance (modify data)
Figure 1: Integrated experimental-numerical approach to predictive modelling of deformation and failure
3
The adopted approach is presented here by describing one of the several concurrent routes
commonly jointly pursued in evaluating a given material sample with the objective of
developing the necessary predictive modelling capability for dealing with that material. An
example of a titanium alloy for aerospace applications is used with the aim of illustrating
several aspects of experimental characterisation based upon uniaxial tensile testing and to
describe the subsequently performed modification of an existing constitutive model that
enabled regularisation of the numerical simulation of observed softening and localization
preceding specimen fracture. In addition, the currently employed methodology for
identification of non-measurable material parameters by means of inverse modelling is
outlined.
2 EXPERIMENTAL METHODOLOGY
The currently employed experimental methodology for mechanical characterisation of
materials is a result of series of carefully conducted modifications of specialised loading rigs
(following the feedback from numerical simulations of actually performed experiments) as
well as the application of ever improving state-of-the-art sensors and data acquisition
equipment. This methodology incorporates a number of different experimental setups thus
enabling generation of a wide range of loading conditions and consequent material states
within the specimen gauge sections under observation24. It is central to the adopted integrated
experimental-numerical approach to predictive modelling of material response to impact
loading to provide detailed information on the conducted experiments in a format readily
employable within the adopted computational framework. This way the selected (existing or
newly developed) numerical algorithms can be validated and modified in order to provide
better agreement with results of conducted experiments. It is expected that the numerical
algorithms capable of accurately simulating the selected range of relevant laboratory
experiments employed in characterisation of the given material can correctly estimate the
behaviour of that material in real components and structures subjected to impact loading.
In this paper a set of uniaxial tensile tests is described as one of the several concurrent
routes (in addition to uniaxial compression, pure shear, biaxial tension, tension-torsion, etc.)
normally undertaken in the adopted experimental approach to characterisation of deformation
and failure of ductile materials. Uniaxial tensile tests are an integral part of the currently
employed experimental methodology and are normally used first in characterisation of
mechanical behaviour of ductile materials16. In the case of characterisation of materials for
application in components and structures subjected to rapidly applied loading the experiments
are performed at several distinct rates of strain in order to enable determination of the effects
(if any) of loading rate on material response. In the presented approach to characterisation of
advanced metallic alloys for aerospace applications by tensile testing a single specimen
geometry (Figure 2) is used at all employed loading rates in order to remove effects of
specimen size and shape/geometry from the obtained results. In addition, the specimen design
was influenced by the desire to reduce the effect of radial inertia on the measured resisting
stress14 and to allow sufficient resolution in calculating the effective stress and strain in the
neck from images recorded by means of high speed photography10,13.
4
10.97 mm 8 mm
M6
R = 1.5 mm
4.6 mm
1 mm 7 mm 2.97 mm d = 3 mm
32.6 mm
Figure 2: Tensile specimen for tests at different rates of strain
The presently utilised experimental technique for uniaxial tensile testing employs the
following loading systems
a) screw-driven loading machine (Hounsfield) for quasi-static loading (controlled cross-
head velocity, vMIN = 0.01 mm/min) (Figure 3a);
b) hydraulic loading machine (in-house) for medium rate loading (controlled piston
velocity, vMAX = 5 m/s) (Figure 3b);
c) Split-Hopkinson-Bar apparatus (in-house) for high rate loading (controlled impactor
velocity, vMAX = 35 m/s) (Figure 4).
v(t) ∆VF v(t)
↓
F(t)
u ∆Vu u
l0 ↓ l0
u(t) A0
A0 a) quasi-static testing b) medium rate testing
Figure 3: Loading rigs for tensile testing at quasi-static and medium rates of strain
Following a straightforward procedure based upon the analysis of electric signals obtained
from a set of calibrated load cells and extensometers used in experiments at quasi-static and
medium loading rates and by ignoring the inertia forces the required mechanical quantities
(extension, extension rate, resisting force) are obtained which are easily comparable against
the results of simulations of conducted experiments. The data acquired in tests at high rates of
strain by means of the Split-Hopkinson-Bar apparatus requires further analysis in order to
obtain the equivalent mechanical quantities as the inertia forces and stress propagation speed
must be considered.
5
bearing input bar

strain gauge specimen output bar
v(t) p0
t t t
arrester loading bar projectile chamber A valve C
bearing input bar
strain gauge specimen output bar
velocity seal chamber B

measuring F(t)
device
Strain gauge Computer
bridge box Amplifiers vINP Transient (IBM PC)
recorders FOUT
c) high rate testing
u(t)
vOUT
t t
Figure 4: Split-Hopkinson-Bar for tensile testing at high rates of strain
The interpretation of data generated by a tensile Split-Hopkinson-Bar relies upon 1D

elastic stress wave analysis18,19. The method has been further enhanced recently in order to
account for stress wave dispersion in loading bars and to enable handling of very ductile
metallic alloys. The method is presented here in detail since it represents a vital part of the
adopted integrated experimental-numerical approach to predictive modelling of material
response to impact loading.
2.1 Analysis of longitudinal elastic waves in thin rods

Let the displacements of particles at position x in a thin rod at time t loaded axially, as
illustrated in Figure 5, be expressed using the following scalar function
u = u ( x, t ) (1)
When the load is applied rapidly the dynamic equilibrium (ignoring the body forces) for a
1D rod follows from the translatorial equilibrium along the rod axis
6
∂σ ∂ 2u
=ρ 2 (2)
∂x ∂t
where σ is the longitudinal (axial) stress and ρ is the density.
dx
u u+du
A, E, ρ σ x
u(x, t) ∂σ
σ+ dx
∂ 2u ∂x
ρ Ad x
∂t 2
Figure 5: Axially loaded rod
For materials which exhibit linear-elastic behaviour when loaded by a force below certain
magnitude the strain is normally small and can be expressed as a function of particle
displacements as follows (compatibility equation)
∂u
ε= (3)
∂x
while the stress can be expressed as a function of strain by (constitutive equation)
σ = Eε (4)
where E is the Young’s modulus thus enabling the transformation of equation (2) as follows
∂ 2u ∂ 2u
c2 = (5)
∂x 2 ∂t 2
where
E
c= (6)
ρ
is the elastic stress wave speed.
The following substitution
ξ = x − ct (7)
η = x + ct (8)
into
u = u (ξ ,η ) (9)
results in
7
∂ 2u ∂ 2u ∂ 2u ∂ 2u 
= + 2 + 
∂x 2
∂ξ 2
∂ξ∂η ∂η 2

 (10)
∂ u 2 ∂ u ∂ u ∂ u 
2 2 2 2

=c  2 −2 + 
∂t 2  ∂ξ ∂ξ∂η ∂η 2 
and after the substitution in the equilibrium equation (5) it follows that
∂ 2u
=0 (11)
∂ξ∂η
thus
u (ξ ,η ) = f (ξ ) + g (η ) (12)
Hence the following D’Alambert’s solution of the 1D wave equation (2)

u ( x, t ) = f ( x − ct ) + g ( x + ct ) (13)
Functions f(x-0) and g(x+0) are obtained from initial conditions as follows
L
1 1
f ( x − 0) = u0 ( x) − ∫ u&0 ( x) d x (14)
2 2c 0
L
1 1
g ( x + 0) = u0 ( x) + ∫ u&0 ( x) d x (15)
2 2c 0
where L is the length of the rod.
The D’Alambert’s solution is illustrated graphically in Figure 6, Figure 7 and Figure 8.
u u(x,0) = f(x-0) + g(x+0)
f(x-0) = g(x+0)
x
g(x+c∆t) f(x-c∆t)
∆t
2∆t
t
Figure 6: Waves propagation in elastic thin rod
8
In Figure 6 the function

dx E
f(x-ct) with speed c = = (16)
dt ρ
denotes the forward travelling wave while the function
dx E
g(x+ct) with speed c = − =− (17)
dt ρ
denotes the backward travelling wave.
For a forward travelling wave (Figure 7) the value of stress is
∂u
σ = Eε = E = E [ f ′( x − ct ) + g ′( x + ct )] (18)
∂x
while the particle velocity is
∂u
v= = −c[ f ′( x − ct ) − g ′( x + ct )] (19)
∂t
where
∂f ( x − ct )
f ′( x − ct ) = (20)
∂ ( x − ct )
and
∂g ( x + ct )
g ′( x + ct ) = (21)
∂ ( x + ct )
while
g ( x + ct ) = 0 ⇒ g ′( x + ct ) = 0 (22)
u(0,t)
x
u
u(x,t0) = f(x-ct0)
t=t0
u x
c(t1- t0)
u(x,t1) = f(x-ct1)
t=t1
x
Figure 7: Forward travelling wave
9
g(x,t)
x
t
t x
f(x,t)
t
t
x
u(x,t)=f(x-ct)+g(x+ct)
t t
Figure 8: Waves superposition
By eliminating f ′( x − ct ) from the above expressions for stress and velocity it follows that
σ v dx
+ = 0 for =c (23)
E c dt
Equivalently, for a backward travelling wave it can be shown that
σ v dx
− = 0 for = −c (24)
E c dt
Recalling that
10
E
c= (25)
ρ
it follows that
dx
σ = ± ρcv for = mc (26)
dt
where ρc is a constant known as the acoustic impedance or characteristic impedance of the
material. In addition, the longitudinal elastic strain in the thin rod can be expressed as follows
σ v dx
ε= =± for = mc (27)
E c dt
2.2 Method of characteristics

D’Alambert’s solution of 1D wave equation (2)
∂ 2u ∂ 2u
c2 = (28)
∂x 2 ∂t 2
obtained in the following form
u ( x, t ) = f ( x − ct ) + g ( x + ct ) (29)
results in the following expression for stress

∂u
σ = Eε = E = E [ f ′( x − ct ) + g ′( x + ct )] (30)
∂x
and particle velocity
∂u
v= = m c[ f ′( x − ct ) − g ′( x + ct )] (31)
∂t
For waves travelling in only one direction it follows that (by eliminating f ′( x − ct ) or
g ′( x + ct ) from (30) and (31))
σ v dx
± = 0 for = ±c (32)
E c dt
For waves travelling in both directions subtracting (31) from (30) gives
σ v
− = 2 f ′( x − ct ) (33)
E c
while summing (30) and (31) gives
11
σ v
+ = 2 g ′( x + ct ) (34)
E c
By integrating
dx
= ±c (35)
dt
it follows that
x − ct = const = A (36)
and
x + ct = const = B (37)
Substituting into (33) and (34) it follows that

α
f ′( x − ct ) = f ′( A) = const = (38)
E
and
β
g ′( x + ct ) = g ′( B ) = const = (39)
E
thus resulting in
σ v 2α dx
− = for =c (40)
E c E dt
and
σ v 2β dx
+ = for = −c (41)
E c E dt
Adding and subtracting the above equations and using E = ρc 2 it follows that
σ = β +α
 (42)
ρcv = β − α 
where α and β are the constants which have units of stress. The constants α and β denote the
dx
magnitude of stress along the same characteristic lines = ±c (hence the name “method of
dt
characteristics”) as illustrated in the example presented in Figure 9.
If at one end of the thin rod shown in Figure 9 the following step load (prescribed velocity)
12
0 t≤0
v ( x, t ) x = 0 =  (43)
v0 t>0
is applied the material particles will gain the velocity v0 and the magnitude of thus initiated
forward travelling stress wave (α wave, β = 0) can be obtained from
ρcv0 = β − α = 0 − α = −α (44)
as follows
σ = β + α = 0 + α = − ρcv0 (45)
where negative sign denotes compression. The α wave reaches the other (free) end of the rod
l
after t = . The conditions at the free end require that
c
σ = β +α = 0 (46)
which determines the magnitude of backward travelling stress wave (β wave)

β = −α (47)
and the velocity of particles

ρcv = β − α = ρcv0 − (− ρcv0 ) = 2 ρcv0 (48)
as follows
v = 2v0 (49)
l
v(t)
x
1 α
dx
=c
dt
β
1
dx
= −c
dt
t
Figure 9: α and β waves
13
2.3 Split-Hopkinson-Bar data analysis

The term “rod” used in the description of 1D wave analysis assumed the uniformity of
material state over the area of any cross section within the domain under consideration. In
reality, the bars used in laboratory experiments are of finite diameter thus embodying the 3D
behaviour of real solid bodies in which longitudinal elastic stress waves coexist with
transverse and surface waves generated by rapidly applied loading at one of the bar ends.
Such interaction results in attenuation and dispersion of the longitudinal waves which are of
primary interest in uniaxial compression and tension experiments. However, in well designed
and manufactured Split-Hopkinson-Bar systems these effects are relatively small compared to
the magnitude of the main longitudinal stress waves commonly generated in experiments.
If precise strain gauges of short gauge length are attached appropriately onto the bars’
surfaces in particular cross sections (as in the example presented in Figure 10) the magnitude
of axial (longitudinal) elastic stress in those cross sections can be calculated directly from the
measured strains
σ ( xG , t ) = α ( xG , t ) + β ( xG , t ) = Eε ( xG , t ) (50)
where E is the Young’s modulus of the bar material and G denotes the position of the strain
gauges station. The velocity of particles at the strain gauges station cross section can be
obtained from
β ( xG , t ) − α ( xG , t ) ε (xG , t )
v( xG , t ) = = (51)
ρc c
where ρ is the density and c is the elastic stress wave speed for the bar material.
l
v(t)
x
α = -ρcv0 σ=0
v=0
σ = -ρcv0
v = v0
β = -ρcv0 σ=0
v = 2v0
0 t≤0
v ( x, t ) x = 0 =
v0 t >0
t
Figure 10: Magnitudes of stress at discrete points of the bar
Using the method of characteristics it is possible to calculate the magnitude of forward and
backward travelling elastic stress waves as functions of position and time α(x,t) and β(x,t)
thus enabling the calculation of total stress in any selected cross section including the
14
downstream end of the bar by measuring the strain in instrumented cross sections.
Consequently, if a tensile specimen is fixed between the two bars of a Split-Hopkinson-Bar
apparatus as illustrated in Figure 11 (where INP and OUT denote “input bar” and “output bar”
respectively) the velocity of particles and the stress in cross sections at the interfaces between
the bars and the specimen can be evaluated. The presented experimental analysis procedure
based upon two sets of strain gauge stations on the input bar allows for attenuation and
dispersion of incident stress waves to be taken into account.
l1
l2 l3 l4
v(t)
x
αINP
β INP
αOUT
t
β OUT
Figure 11: Split-Hopkinson-Bar experimental setup
The data obtained at strain gauge stations (instrumented cross sections on input and output
loading bars) can be used to calculate the relationship between specimen’s extension
∆ u (t ) = u xMAX
INP
( ) (
, t − u xMIN
OUT
,t )
t t
= ∫ v(x ,τ )d τ −∫ v(x ,τ )dτ
INP OUT
MAX MIN
0 0 (52)
t
β (x ,τ ) − α (x ,τ )
INP
β (x
INP t OUT
) (
,τ − α xMIN
OUT
)
,τ
= ∫
0
MAX
ρc
dτ −MAX
∫
0
MIN
ρc
dτ
and its corresponding resistance force as the following functions of time
F (t ) =
1
2
[ ( ) (
INP
F xMAX )]
, t + F xMIN
OUT
,t
(53)
=
1
2
{ [( ) (
AINP α xMAX
INP
)]
, t + β xMAX
INP
[(
, t + AOUT α xMIN
OUT
) (
, t + β xMIN
OUT
,t )]}
where MAX denotes the rightmost cross section on the input bar and MIN denotes the
leftmost cross section output bar (the interfaces between the bars and the specimen), F are the
15
total internal forces in the bars at the interfaces with the specimen, while AINP and AOUT are
the cross section areas of the input and output bar respectively. The magnitudes of forward
and backward travelling waves at the interfaces between the bars and the specimen are
calculated using the algorithm presented in Figure 12. The algorithm is the result of recent
enhancement that allows for multiple reflections from various interfaces within the Split-
Hopkinson-Bar system. This way it became possible to quantify the response of very ductile
materials.
l1 l3 l4
l2
v(t) EINP, ρINP, cINP EOUT, ρOUT, cOUT
G1 G2 G3 x
t1 α INP
t2
t3
t5
βINP t6
α INP α OUT
t4
βOUT
βINP t7
α OUT
t t8
(
α xMAX
t < t4 
INP
) [
, t = σ xGINP , t − (t3 − t1 ) ] α OUT
( ) [ ] [ ]
1
 β xMAX , t = σ xG 2 , t + (t3 − t2 ) − σ xG1 , t + (t3 − t 2 ) − (t2 − t1 )

INP INP INP
α (x INP
MAX , t ) = σ [x
INP
G1 , t − (t − t )]− β [x
3 1
INP
MAX ]
, t − (t 4 − t3 )
≤ t
, t ) = σ [x , t + (t − t )]− {σ [x ] [ ]}
t4
 β (x , t + (t3 − t 2 ) − (t 2 − t1 ) − β xMAX , t − 2(t 2 − t1 )
INP INP INP INP
MAX G2 3 2 G1
( ) [
α x OUT , t = σ xGOUT
t < t7  MIN
, t + (t6 5 − t )]
( )
3
 β xMIN , t = 0
OUT
( ) [
t7 ≤ t < t8  MIN
] [
, t + (t6 − t5 ) + α xMIN
OUT
, t − (t7 − t5 ) ]
( )
3
 β xMIN , t = 0
OUT
( ) [
t8 ≤ t  MIN
, t + (t6 − t5 ) + α xMIN ] [
OUT
, t − (t7 − t5 ) ]
( ) [ ]
3
 β xMIN , t = −α xMIN , t − (t8 − t5 )

OUT OUT
Figure 12: Enhanced Split-Hopkinson-Bar analysis
16
A new software (ELAN2) for analysis of data generated by Split-Hopkinson-Bar systems

has been developed that incorporates the presented analysis procedure thus providing insight
into the mechanical characteristics of tested materials immediately following the triggered
high-speed data acquisition. The software has also been designed to provide data on resulting
specimen boundary loading conditions in the format readily applicable in numerical
simulation of the conducted experiments. This way a straightforward direct comparison
procedure is enabled that has been exploited within the newly developed inverse modelling
framework designed for identification of non-measurable parameters of constitutive models.
A typical basic set of results of uniaxial tensile tests at high rates of strain (Figure 13)
comprises detailed information on specimen boundary and loading conditions ready for use in
numerical modelling of conducted experiments.
30 3500
28 L09
26 L10
L11 3000
24 T09
22 T10
Extension velocity (m/s)
Nominal strain rate (1/s)

2500
20 T11
18
2000
16
14 Specimen fracture
1500 L09
12
Specimen fracture L10
10
L11
8 1000
T09
6 T10
4 500 T11
2
0 0
1.0E-04 1.5E-04 2.0E-04 2.5E-04 3.0E-04 3.5E-04 4.0E-04 0 5 10 15 20 25 30
Time (s) Nominal strain (%)
5.0E-03 10
L09 L09
4.5E-03 L10 9
L10
L11 L11
4.0E-03 8
T09
T09
3.5E-03 T10 7
T10
T11
Extension (m)
3.0E-03 6 T11
Force (kN)
2.5E-03 5
2.0E-03 4
1.5E-03 3
Specimen fracture
1.0E-03 2
5.0E-04 1
0.0E+00 0
1.0E-04 1.5E-04 2.0E-04 2.5E-04 3.0E-04 3.5E-04 4.0E-04 0.00E+00 2.50E-04 5.00E-04 7.50E-04 1.00E-03 1.25E-03 1.50E-03 1.75E-03
Time (s) Extension (m)
10 2000
L09
9 1800
L10
8 L11 1600
T09 1400
7
True stress (MPa)
T10 L09
1200 L09
6 L10
Force (kN)
T11 Onset of specimen fracture L10

5 1000 L11
Onset of material instability L11
800 T09
4 T09
T10
3 600 T10
T11
T11
2 400
1 200
Specimen fracture
0 0
1.0E-04 1.5E-04 2.0E-04 2.5E-04 3.0E-04 3.5E-04 4.0E-04 0 5 10 15 20 25 30 35 40 45 50 55 60
Time (s) True strain (%)
Figure 13: Typical basic set of results of uniaxial tensile tests at high rates of strain
17
A typical set of results of uniaxial tensile tests at three distinct rates of strain comprises the
applied prescribed velocity and the measured resistance force boundary conditions. An
example of such results set is given in Figure 14 together with the corresponding
recommended numerical model.
18 10
16 9
HR
14 8
7
12 HR
6
Force (kN)
10
5
8
4
6
3
4
2
2 1
0 0
0.0E+00 2.0E-05 4.0E-05 6.0E-05 8.0E-05 1.0E-04 1.2E-04 1.4E-04 0.0E+00 2.0E-05 4.0E-05 6.0E-05 8.0E-05 1.0E-04 1.2E-04 1.4E-04
Time (s) Time (s)
0.7 9
MR MR
8
0.6
7
0.5
6
Force (kN)
0.4 5
0.3 4
3
0.2
2
0.1
1
0 0
0.0E+00 1.0E-03 2.0E-03 3.0E-03 4.0E-03 5.0E-03 6.0E-03 7.0E-03 8.0E-03 0.0E+00 1.0E-03 2.0E-03 3.0E-03 4.0E-03 5.0E-03 6.0E-03
Time (s) Time (s)
7.0E-06 8
QS QS
6.0E-06 7
6
5.0E-06
5
Force (kN)
4.0E-06
4
3.0E-06
3
2.0E-06
2
1.0E-06 1
0.0E+00 0
0 50 100 150 200 250 300 350 0 50 100 150 200 250 300 350
Time (s) Time (s)
10
L09
9 L10
L11
8 T09
T10
7 T11
v(t) FR(t) 6
L07
Force (kN)
L08
L12
5
T07
4 T08
T12
3 L04
L05
2 L06
T04
1 T05
T06
0
0.00E+00 2.50E-04 5.00E-04 7.50E-04 1.00E-03 1.25E-03 1.50E-03 1.75E-03 2.00E-03
Extension (m)
Figure 14: Typical set of idealised results of uniaxial tensile tests at three different rates of strain
18
The specimen boundary conditions can be used to estimate the material’s response in the
form of constitutive relation. Such calculation must take into account the geometry of the
specimen and the dynamic effects resulting in radial inertia as well as the triaxiality
conditions following the onset of necking.
The use of high frequency laser extensometers and high speed photographic equipment
enables further quantification of the induced material behaviour40. As a result, continuous
measurements of the reduction in cross section area and the recording of full field surface
displacement and temperature can be performed thus providing additional information on the
material behaviour for comparison against the results of numerical simulations of conducted
experiments. In addition, all specimens can be examined after testing using the X-ray
diffractometry as well as the TEM/EBSD microscopy (Figure 15) in order to relate
macroscopic observation and measurements to observable/measurable effects at microscopic
scale. Equivalent data can be gathered across a range of different types of experiments which
is instrumental to the development of new improved physically based predictive modelling
tools.
Figure 15: Microscopic observation of cracking in a ductile metallic alloy
19
3 COMPUTATIONAL FRAMEWORK
At microscopic scale matter consists of atoms that are held together by bonds resulting
from the interaction of electromagnetic fields and are grouped into molecules thus forming
chains, planes, crystals, fibres, grains. Engineers, however, mostly have to deal with pieces of
matter at a much larger scale. When observed at macroscopic scale, matter appears (to the
naked eye) continuous, hence the name continuum. In such context, the world appears to be
composed of dense matter or solids (Earth’s crust, man-made structures) and rare matter or
fluids (liquids, gases) that are continually and continuously interacting. One mathematical
description of such world regards matter as indefinitely divisible and is known as continuum
mechanics. The basis of this theory is defined at a material point or particle - an infinitely
small part of the continuum under consideration still retaining all characteristics of the
continuum as a whole. However, real materials are not continuous and the quantities (e.g.
strain, strain rate, temperature, stress) chosen to describe the observed behaviour of continua
actually represent averages over a representative volume element (RVE) at a macroscopic
scale. The actual size of RVE employed in predictive modelling of real phenomena depends
upon the geometric characteristics of domains under consideration (e.g. sharp notches, tight
radii, etc.) and the type of observed material behaviour (e.g. inviscid flow, localisation,
cracking, etc.). As a result, the required size of RVE can be reduced to a dimension of the
same order as the size of material’s constituents (e.g. grains, fibres, etc.) and the assumptions
of continuity and homogeneity upon which the continuum mechanics theory is based are no
longer applicable. Similarly, at a larger scale mechanical systems consist of numerous
components forming assemblies of discrete interacting bodies. Hence – the Discrete Element
Method.
3.1 The Discrete Element Method

The Discrete Element Method (DEM) is a family of related numerical techniques designed
to solve those problems in applied mechanics that exhibit large discontinuous geometrical and
material behaviour15. In order to predict the behaviour of such mechanical systems the
employed numerical methodology should incorporate
- combined discontinuum/continuum (multi-scale) representations of matter in which the
discontinuities can be characterised fully, as embedded or as smeared at all scales,
- efficient geometric intersection search algorithms that can support the discontinuous
(discrete) representation of matter thus involving extremely large numbers of interacting
entities,
- representations of interaction across discontinuities for all types of matter (solid-solid,
solid-fluid, fluid-fluid),
- representations of the transition from continuum to discontinuum at a given scale.
Extensive research has been conducted over the past two decades8 that has provided
fundamental contributions for inclusion of these techniques into newly developed software for
numerical analysis of mechanical systems. Convalescing from this approach, in the
implementation presented in this paper the geometry of discrete elements is described
(Figure 16) by either
20
a) single rigid body primitives (e.g. spheres, ellipsoids, superquadrics, regular and
arbitrary polyhedrons, etc.), or
b) single deformable bodies (which are further discretised using finite element method or
smoothed particle hydrodynamics technique), or
c) recursive use of either or both of the first two types of approximation while forming
assemblies with a specific interaction laws.
c)
b)
a)
Figure 16: Discrete element discretisation of solid bodies
This computational framework is ideally suited for research into predictive modelling of
deformation and failure of solids and solid assemblies as it contains all the necessary
ingredients required to simulate successfully the transition from continuum to discontinuum.
It has as such been adopted to accompany the presented experimental methodology thus
enabling an integrated approach to predictive modelling of deformation and failure in metallic
alloys subjected to impact loading. In the adopted approach presented in this paper large
mechanical systems (e.g. automotive structures in crash simulations, aeroengine fan structure
in birdstrike events, etc.) are treated as assemblies of continuous deformable solids.
The current implementation (DEST - Discrete Element Simulation Tools36) of this
methodology relies upon the updated Lagrangian description of motion and employs the latest
computational technology (object-oriented, function call overload optimised and parallelised
code) in the attempt to offer an accurate as well as efficient solution to the following
governing equations
- kinematic compatibility
∂x ∂v
 1
d = 2 l + l(T
)
F= ⇒ l = F& F −1
= ⇒  (54)
∂X ∂x 1
w = l − l T
 2
( )
- conservation of mass
21
ρ& + ρ div[u& ] = 0 (55)
- conservation of momentum
div[σ EXT − σ INT ] + ρ b = ρ u
&&
 ∀p ∈ (Ω ) (56)
σ ij = σ ji 
- conservation of energy
ρe& = σ : d + r − div[q ] (57)
- constitutive relation (rate form)
σ (t ) = Σ [σ (τ ), ε (τ ), d (τ ), w (τ ), θ (τ ), {V }], − ∞ < τ ≤ t
* (58)
- boundary conditions
σEXT n = t ∀p ∈ (∂ t Ω )
 (59)
u = ub ∀p ∈ (∂ u Ω ) 
- contact interface impenetrability condition

∆ vN = v A ⋅ nC + v B ⋅ nC = v NA − v NB ≤ 0 (60)
- initial conditions
u(0) = u0 
 (61)
v (0) = v 0 
where p denotes a material particle within the domain Ω or on the domain boundary ∂Ω, F is
the deformation gradient, X and x are the reference and current position of p respectively, l is
the velocity gradient, v is the velocity of p, d is the rate of deformation tensor, w is the spin
tensor, ρ is the material density, u is the displacement of p, σ is the Cauchy (true) stress
(where EXT denotes external and INT denotes internal stress), e is the internal energy density,
r is the internal heat production density (rate of radiant heat supply), q is the surface heat flux,
σ* is the objective rate of Cauchy stress tensor, Σ is the functional of deformation history, ε is
the true strain tensor, θ is the temperature, {V} is the set of internal variables, n is the surface
normal vector on the domain’s boundary, t is the surface traction vector along ∂t Ω portion of
the boundary (Neuman’s boundary condition), ub is the prescribed displacement along ∂u Ω
portion of the boundary (Dirichlet’s boundary condition), ∆ v N is the relative velocity of
contacting entities in the direction normal to the common contact surface defined by its
normal nC thus defining the initial value problem of discrete solid mechanics.
The algorithms presented in this paper rely upon the finite element sub-discretisation of
discrete elements (distinct continuous solids) within the domain under consideration. This
22
enables utilisation of computational methodology that has proved itself in dealing with the
bulk of materials’ response to loading preceding the onset of material instability38 leading to
strain localisation, cracking and component fracture (and potentially loss of structural
integrity). The adopted approach allows for other numerical techniques (Mesh Free Galerkin
- MFG, Smoothed Particle Hydrodynamics - SPH, Finite Volume Method - FVM, Boundary
Element Method - BEM) to be utilised within the same code. In addition, due to the short
lasting nature of impact events under consideration the solution to the governing partial
differential equations in the time domain is sought by means of the leap-frog central
difference explicit time integration scheme.
3.2 Constitutive modelling of rate dependent material behaviour and strain localisation
The accuracy of numerical predictions of material failure by cracking and fracture is
mainly dependent on the performance of employed constitutive models. In the case of impact
loading it is important that constitutive models are capable of correctly simulating material’s
rate dependent behaviour if any such behaviour is observed in experiments. A number of
models have been proposed that are capable of simulating rate dependent behaviour by using
simple empirical or more advanced physically based constitutive equations27,31. Most of such
models can successfully represent the observed response up to the onset of material
instability. In order to enable simulation of the observed response beyond the onset of
material instability it is important that the employed constitutive models can exhibit strain
softening6. However, the consequent numerical difficulties posed by the change in nature of
underlying differential equilibrium equations from hyperbolic to elliptic results in the inability
of the adopted solution method to represent accurately the fracture as solutions become
hugely dependent on the size of the RVE defined by spatial discretisation. Consequently, the
energy dissipation within the localisation zone vanishes as the size of elements within the
zone approaches zero value4,5,37. This obviously cannot represent fracture as fracture always
involves a certain amount of energy dissipation. As a result, a number of regularisation
techniques have been proposed23,33 in the attempt to resolve these problems.
In this paper, modifications and further regularisation of a constitutive model which
couples the temperature dependent isotropic elasto-visco-plasticity and scalar damage
mechanics2 are presented. The constitutive model has been adapted in order to remove the
dependency of its predictions on the size of the RVE determined by spatial discretisation.
This has been achieved by introducing a length scale into the equations defining the evolution
of the scalar damage variable within the model. In this model the material response to
deformation is formally expressed in the following rate form
[ ]
σ (t ) = Σ σ (τ ), ε (τ ), d (τ ), w (τ ), θ (τ ), α (τ ), κ (τ ), {D, D& } , − ∞ < τ ≤ t
* (62)
where α is the back stress, κ is the isotropic hardening and D is the scalar damage variable.
The algorithmic implementation of this functional follows from the energy considerations
encompassed by the fundamental principles of thermodynamics with internal variables. In
this approach, the first principle (energy balance)
23
ρe& = σ : d + r − div[q ] (63)
and the second principle

q  r
ρs& + div   − ≥ 0 (64)
θ  θ
are combined to give the Clausius-Duhem inequality
q  1
ρs& + div   − ( ρe& − σ : d + div[q ]) ≥ 0 (65)
θ  θ
The introduction of the specific free energy ψ (the thermodynamic potential or the
Helmholtz free energy per unit mass or the reversible energy that can be transformed into
work) defined as
ψ = e − θs (66)
and noting that

q  1 1
div   = div[q ] − 2 q ⋅ grad[θ ] (67)
θ  θ θ
yields the following expression of the Clausius-Duhem inequality
( 1
)
σ : d − ρ ψ& + sθ& − q ⋅ grad[θ ] ≥ 0
θ
(68)
If the Helmholtz free energy is defined as the following function

ψ = ψ (ε,θ , {V }) (69)
of the observable ( ε,θ ) and the internal ( {V }) variables which define the local material state
where
t
ε (t ) = ∫ d (τ ) d τ (70)
0
is the true strain (work conjugate to Cauchy stress), the Clausis-Duhem inequality can be
expressed as follows
 ∂ψ   ∂ψ & 1
σ − ρ  : d − ρ s + θ − AiV&i − q ⋅ grad[θ ] ≥ 0 (71)
 ∂ε   ∂θ  θ
The description of material inelastic behaviour follows from the multiplicative
decomposition of the deformation gradient into its elastic and inelastic parts
F = F el F inel (72)
24
and the consequent additive decomposition of the velocity gradient

l = F& F −1 = F& el F inel F −1 + F el F& inel F −1 ≡ l el + l inel (73)
from which it follows that

 d = d el + d inel
l= d+ w ⇒  (74)
w = w + w
el inel
where
 d inel = d pl + d θ + d Vi
 inel θ
(75)
w = w + w + w
pl Vi
in which the superscript pl denotes plastic or visco-plastic, the superscript θ denotes thermal
while the superscript Vi denotes the portion of the inelastic rate of deformation and inelastic
spin associated with dissipation processes described by the internal variables Vi. The
associated Helmholtz free energy density can then be expressed as follows
ψ = ψ (ε − ε inel ,θ , {V }) = ψ (ε el ,θ , {V }) (76)
thus resulting in the following form of the Clausius-Duhem inequality

 ∂ψ  el  ∂ψ & 1
 σ − ρ el  : d + σ : d − ρ  s + θ − AiV&i − q ⋅ grad[θ ] ≥ 0
pl
(77)
 ∂ε   ∂θ  θ
which represents the basis for development of constitutive models for real solid materials.
The thermoelastic constitutive laws (78) and (79) follow from the observation that the
Clausis-Duhem inequality holds for any particular d el thus
∂ψ
σ=ρ (78)
∂ε el
which in case of a thermal deformation defined by d pl = 0 , grad[θ ] = 0 and V&i = 0 yields
∂ψ
s=− (79)
∂θ
The constitutive models for non-linear material response to loading follow from the
application of thermoelastic constitutive laws (78) and (79) which reduce the Clausis-Duhem
inequality to the following form
Φ M + ΦT = {σ : d pl − AiV&i }−  q ⋅ grad[θ ] ≥ 0
1 
(80)
θ 
where Φ M is the mechanical (intrinsic) dissipation and ΦT is the thermal dissipation.
The introduction of a dissipation potential (a convex scalar valued function of flux
25
variables) or its complementary potential (a convex scalar valued function of thermodynamic

forces)
ϕ = ϕ  d pl , {V& }, 
 q
ϕ * = ϕ * (σ , {A}, grad[θ ]) (81)
 θ
which can be decoupled (if suitable)
ϕ = ϕ M (d pl , {V& }) + ϕT  
q
ϕ * = ϕ M* (σ , {A}) + ϕT* (grad[θ ]) (82)
θ 
enables the definition of corresponding equivalent evolution laws
∂ϕ M ∂ϕ M*
σ= d pl = (83)
∂d pl ∂σ
∂ϕ M ∂ϕ *
Ai = − − V&i = M (84)
∂V&i ∂Ai
∂ϕT q ∂ϕT*
grad[θ ] = − − =
q θ ∂ (grad[θ ]) (85)
∂ 
θ 
where Ai is the thermodynamic force associated to the internal variable Vi.
For a linear relation between the heat flux and temperature gradient the thermal dissipation
potential can be expressed as follows
1
ϕT* = C ⋅ grad[θ ] ⋅ grad[θ ] (86)
2
which yields
q ∂ϕT*
− = = C grad[θ ] (87)
θ ∂ (grad[θ ])
where C is the tensor of material diffusion properties (constants) which in the case of isotropic
materials reduces to a scalar k that varies inversely with respect to temperature thus
q = −k grad[θ ] (88)
As a result, the first principle of thermodynamics (the conservation of energy) expressed as

follows
ρ
d
dt
[( ) ]
ψ ε el ,θ , {V } + θs = σ : d + r − div[q ] (89)
yields
26
σ : d el + AiV&i + ρθ s& = σ : d + r − div[q ] (90)
By noting that
d  ∂ψ (ε el ,θ , {V }) ∂ 2ψ ∂ 2ψ & ∂ 2ψ &
s& = −   = − : d el
− θ− Vi (91)
dt  ∂θ  ∂ε el ∂θ ∂θ 2 ∂Vi ∂θ
from which it follows that
1 ∂σ ∂s & 1 ∂Ai &
s& = − : d el − θ− Vi (92)
ρ ∂θ ∂θ ρ ∂θ
and by introducing the specific heat
∂s
C =θ (93)
∂θ
and by noting that
div[q ] = − k div[grad[θ ]] = − k ∆θ (94)
where ∆ is the Laplacian operator, it follows that the first principle of thermodynamics can be
expressed as follows
 ∂σ el ∂Ai & 
AiV&i −θ  :d + Vi  + ρCθ& = σ : d pl + r − k ∆θ (95)
 ∂θ ∂θ 
Rapidly applied loading of materials results in adiabatic heating which can be described by
k ∆θ = 0
 (96)
r =0 
If the thermomechanical coupling is ignored
∂σ el 
: d = 0
∂θ
∂Ai &  (97)
Vi = 0 
∂θ 
and the non-recoverable energy stored in the material is considered constant
AiV&i = const. (98)
the temperature evolution law can be expressed as follows

σ : d pl
θ& = (99)
ρC
In the case of polycrystalline isotropic titanium alloy under consideration in this paper the
27
complementary dissipation potential was assumed to be equal to the following Von Mises
yield criterion with visco-plastic hardening and scalar damage
ϕ * = F (σ , α, κ , d pl , w pl ,θ , D, D& ) = 0 (100)
from which the evolution laws for plastic rate of deformation and spin as well as the evolution
of kinematic and isotropic hardening, adiabatic temperature and scalar damage are derived.
In order to allow for the simulation of deformation involving large rotations the adopted
updated Lagrangian formulation based upon explicit time integration of the underlying partial
differential equations of equilibrium employs the Jaumann objective stress rate formulation of
the Cauchy stress
T
   
σ tc+ ∆ t = σ t + ∆ t  w ∆ t  ⋅ σ t + σ t ⋅  w ∆ t  ∆ t (101)
 t+ 2   t+ 2 
and the back stress
T
   
α c
t +∆ t = αt + ∆ t  w ∆ t  ⋅ αt + αt ⋅  w ∆ t  ∆ t (102)
 t+ 2   t+ 2 
This enables the calculation of the trial deviatoric stress
  D& ∆ t  1  1   
σ t′+TR∆ t = σ tc+ ∆ t 1 − min1, t [ ]
 − tr σ tc+ ∆ t I + 2 µ (1 − Dt ) d ∆ t − tr d ∆ t  I  ∆ t
 t+  (103)
  1 − Dt  3  2 3  2 
t+
as well as the calculation of a hardening-minus-recovery form of the trial back stress

   
αtTR+ ∆ t = αtc+ ∆ t 1 −  rs (θ t ) + rd (θ t )d pl∆ t α tc+ ∆ t ∆ t  (104)
  
t+
2 
and the scalar valued trial isotropic hardening

   
κ tTR+ ∆ t = κ t 1 −  Rs (θ t ) + Rd (θ t )d pl∆ t κ t ∆ t  (105)
  t+
2  
where
2  pl 
d pl∆ t =  d ∆ t : d pl∆ t  (106)
t+
2 3  t + 2 t + 2 
is the Von Mises effective rate of plastic deformation and
α tc+ ∆ t =
2
(
3 c
αt + ∆ t : αtc+ ∆ t ) (107)
28
is the effective back stress, while I is the identity tensor, µ is the elastic shear modulus and ∆t
is the time increment between the reference and current configurations in the employed
explicit time integration algorithm.
The elastic trial effective stress is checked for yielding using the following yield criterion
ϕ * = F (σ , α, κ , d pl , w pl ,θ , D, D& ) =
2
(
3 TR
)
+∆ t − σ Y = 0
ξ t + ∆ t : ξ tTR (108)
in which
2 TR
ξ tTR ′TR
+∆ t = σt +∆ t − αt + ∆ t (109)
3
is the difference between the deviatoric trial stress and the trial back stress while
  d ∆t  
 −1   TR 
t+
σ Y = (1 − Dt )V (θ t )sinh  2
 + Y (θ t ) + κ t + ∆ t  (110)
 f (θ t ) 
   
is the yield stress and
2 
d =  d ∆t : d ∆t  (111)
t+
∆t 
3 2 + + 
2 
t t
2
is the effective rate of deformation.

The above formulation of yield stress uses the following scalar valued functions of
temperature
 C 
V (θ t ) = C1 exp − 2  (112)
 θt 
C 
Y (θ t ) = C3 exp 4  (113)
 θt 
 C 
f (θ t ) = C5 exp − 6  (114)
 θt 
to describe the initial yield stress and the following functions of temperature
 C 
rd (θ t ) = C7 exp − 8  (115)
 θt 
29
C 
h(θ t ) = C9 exp 10  (116)
 θt 
 C 
rs (θ t ) = C11 exp − 12  (117)
 θt 
to describe the kinematic hardening and these functions of temperature
 C 
Rd (θ t ) = C13 exp − 14  (118)
 θt 
C 
H (θ t ) = C15 exp 16  (119)
 θt 
 C 
Rs (θ t ) = C17 exp − 18  (120)
 θt 
to describe isotropic hardening processes. In these functions Ci are the parameters which
need to be determined from experiments. This formulation enables the incorporation of
thermal activation mechanisms in dislocation motion into the constitutive model.
If the trial stress violates the yield criterion the stress correction (radial return) is performed
by calculating the plastic multiplier
3 TR
2
( )
+∆ t − σ Y
ξ t + ∆ t : ξ tTR
∆γ = (121)
2
2 µ (1 − Dt ) + [h(θ t ) + H (θ t )]
3
and the flow vector
ξ
n=
3 (122)
(ξ : ξ )
2
after which it follows that
σ t + ∆ t = σ t′+TR∆ t − 2µ (1 − Dt ) ∆ γ n + pI (123)
αt + ∆ t = αtTR+ ∆ t + h(θ t ) ∆ γ n (124)
30
2
κ t + ∆ t = κ tTR+ ∆ t + H (θ t ) ∆γ (125)
3
where
]1 − min1, 1D− ∆Dt  + 3K (1 − D ) 13 tr d 

&
1
3
[
p = tr σ tc+ ∆ t t
t
t+
∆t  ∆t (126)
  t   2 
is the pressure in which K is the bulk modulus.
The change in temperature adopted in this constitutive model for simulation of short
lasting events follows from the assumption of the adiabatic (non-conducting) conditions in
which the inelastic work gets dissipated as heat. As a result the following evolution equation
for temperature is adopted
θ& =
C19
ρC20
(
σ t + ∆ t : d tpl+ ∆ t ) (127)
thus enabling temperature update

θ t + ∆ t = θ t + θ&t + ∆ t ∆ t (128)
where C19 is the fraction of inelastic work dissipated through heat (usually 0.9) and C20 is the
material’s specific heat coefficient while
d tpl+ ∆ t = ∆ γ n (129)
is the plastic rate of deformation.

The evolution of damage follows from the assumption of uncoupled spherical void growth
in a rate dependent material as follows
 2(2m − 1) p  1 
D& = sinh   − 1 + D  d tpl+ ∆ t (130)
 2(2m + 1)σ ef′   (1 − D )
m

where
3
σ ef′ = (σ t′+ ∆ t : σ t′+ ∆ t ) (131)
2
is the Von Mises effective stress and
2
∆γ
3 (132)
d tpl+ ∆ t = d tpl +
∆t
is the effective plastic rate of deformation while
31
m = C21 Lc + C22 (133)
is the void growth function and Lc is the characteristic length (the size of integration domain
in the direction of maximum principal rate of plastic strain).
The void growth function controls the rate at which damage grows such that the correct
amount of energy is dissipated per unit volume within the localisation zone thus correctly
representing the amount of fracture energy released per unit area of crack surface in the given
material. The incorporation of this function made it possible to further regularise the
behaviour of the presented constitutive model in the softening regime. The simulation of
experiments at high rates of strain prior to this modification of the constitutive model
exhibited enormous mesh dependency as illustrated in Figure 17 wherein the results of
simulations with several different mesh densities are presented. The results show that the
amount of energy dissipated per until volume in the coarse mesh is greater than the amount of
energy dissipated per unit volume of the fine mesh. This confirmed that the energy dissipated
per unit volume in the original formulation of the constitutive model would approach zero as
the element size decreased to zero value despite the implicit regularisation performed by
incorporating the visco-plastic response into the constitutive model.
The identification of non-measurable material parameters C21 and C22 was performed by
employing the inverse modelling technique presented in the next section.
Exp
10
Mesh 1
9 Mesh 2
8 Mesh 3
Mesh 4
Mesh 1
7 Mesh 5
6 Mesh 6
Force (kN)
5
4
3
2
1
0
0.0E+00 1.0E-03 2.0E-03 3.0E-03 4.0E-03 5.0E-03
Extension (m)
Mesh 2
Mesh 4
Figure 17: Mesh dependency prior to damage evolution regularisation
32
3.2 Inverse modelling for identification of non-measurable constitutive parameters

The application of constitutive models to real materials requires the identification of
material parameters to allow for quantitative representation of observed behaviour. These
parameters may be found in a number of ways. For some it is possible to design suitable
experiments to enable their measurement directly24 or indirectly2. For others it is not, or else
the experimental method inherently introduces too much uncertainty to allow sufficient
confidence in the results. For example, temperature effects associated with experiments at
high rates of strain are very difficult to quantify thus making accurate parameters
identification highly challenging28. Alternatively, the parameters can be identified by solving
the corresponding inverse problem17,29,30.
The inverse modelling approach explicitly tries to match the outcomes of laboratory
experiments and numerical simulations. The virtue of a given set of parameters is described
by the use of an objective function, which measures the difference between experimental
results and the corresponding numerical results. The problem of finding the optimal
parameters may then be cast as an inverse problem, where the objective function is to be
minimised. The desired parameter values are found at the objective function’s global
minimum. Evaluation of a single point of the objective function (a single solution of the
direct problem) requires the comparison of experimental and numerical results. The
generation of experimental results in the laboratory is a time-consuming and expensive
exercise, but only ever need be done once for each material under consideration. Numerical
results are obtained by performing numerical simulations, which may also be very time
consuming and need to be repeated many times throughout the parameters identification
process. A novel approach implemented to solve this problem, employs a heuristic, gradient-
based search method and is presented here in some detail.
Let the set of material parameters which need identifying be denoted by a point p in
parameter hyperspace P
p = ( p1 , p2 , K , pn ) ⊂ P (134)
Given a Dirichlet’s boundary condition (i.e. the velocity boundary condition v(t) measured
during the experiment) then the Neuman’s boundary condition (total resisting uniaxial force f
provided by the specimen under tensile loading) for the given domain is obtained numerically
and its variation with the parameters is expressed as follows
f ( p ) = f ( p, v (t ), t ) (135)
An objective function that measures the difference between the experimentally obtained
and numerically simulated reaction force histories in a least squares sense is formulated as
follows
1 2
J (p ) = f (p ) − fˆ (136)
2 Y
The goal of the inverse modelling exercise is the optimal set of parameters p * , achieved
when the objective function is at a global minimum and can be expressed by
33
p* = argmin[J ( p )] (137)
or by
∃p* ⊂ P ∈ ∀p ⊂ P : J ( p) ≥ J ( p* ) (138)
Initially, a value of the objective function is obtained at the starting point (a given
parameter set, estimated based on experience or guessed). From the initial position in
parameter space the computation towards the optimal one is carried out by the search
algorithm consisting of two phases at each increment. Firstly, an exploratory phase, where
the surrounding topology of the objective function is determined and secondly, an update
phase, where the current parameter set is updated and the size of the exploration subdomain is
modified if necessary.
The numerical and experimental results are linearly interpolated and the sum of squares of
differences is calculated. This is the objective function used to evaluate the parameters.
Inverse modelling simulations are set to terminate when a given tolerance is reached. This
gives reasonable execution times while reaching an accuracy far outweighed by other sources
of error in the experiment. Normally, optimisations are repeated from a variety of starting
positions in order to ameliorate the effects of local minima. The procedure is illustrated
schematically in two-parameter space in Figure 18.
(x*x+y*y)+AS*sin(oS*x)*AS*cos(oS*y)
200
180
160
140
250 120
200 100
150 80
100 60
40
50
20
0
0
-50
10
7 89
4 56
1 23
-10 -9 -8
-7 -6 -5 -2-10
-4 -3 -2
-1 0 1 -4-3
2 3 4 -7-6-5
5 6 7 -9-8
8 9 10 -10
Figure 18: Schematic illustration of the update phase in the optimization algorithm
The material parameters C21 and C22 which define the evolution of damage upon the mesh
size were identified by means of the outlined inverse modelling technique using the
34
information from uniaxial tensile experiments at high rates of strain. The results presented in
Figure 19 show that the developed regularisation algorithm is capable of ensuring that the
same amount of energy is dissipated per unit volume regardless of the size of RVE in the
simulation.
10
Exp
Mesh 1
9
Mesh 2
Mesh 3
8 Mesh 4
Mesh 5
Mesh 6
7
6
Force (kN)
0
0.0E+00 2.0E-04 4.0E-04 6.0E-04 8.0E-04 1.0E-03 1.2E-03 1.4E-03 1.6E-03 1.8E-03
Extension (m)
Figure 19: Results of developed regularization of softening by coupling visco-plasticity and damage mechanics
4 CONCLUSIONS
Elements of an integrated experimental-numerical approach to predictive modelling of
deformation and failure in ductile metallic alloys for aerospace applications have been
presented. The adopted methodology integrates several mature experimental and numerical
technologies and continues to be improved by incorporating state-of-the-art techniques.
This paper presents one of the routes for mechanical characterisation of ductile metallic
alloys involving observation and measurements during uniaxial tensile testing at high rates of
strain as well as the description of a suitable constitutive model and the identification of its
parameters. The developed algorithm for regularisation of material instability further
enhances the accuracy of predictive modelling of strain localisation, cracking and fracture in
ductile metallic alloys.
In the adopted integrated approach to predictive modelling of material response to impact
loading these activities are normally executed concurrently with a number of different types
of experiments followed by numerical simulations and inverse modelling of conducted
experiments. The aim of such elaborate research programmes is to enable the accurate
simulation of material states created in laboratory conditions in order to confidently
extrapolate the gained knowledge onto predictive modelling of the response of real
automotive and aerospace structures to impact loading.
35
ACKNOWLEDGMENTS
The support from the European Commission (grants G3RD-CT-2000-00186 and G4RD-
CT-2000-00395) involving major European civil aircraft and automobile manufacturers, the
UK Engineering and Physical Sciences Research Council EPSRC (grant GR/R92691/01) and
Rolls-Royce plc are gratefully acknowledged.
REFERENCES
[1] R.J. Asaro, Micromechanics of Crystals and Polcrystals, Advances in Applied
Mechanics, Vol.23, p.1-15, 1983
[2] D.J. Bammann, M.L. Chiesa, M.F. Horstemeyer and L.I. Weingarten, Failure in ductile
materials using finite element methods, in: Structural Crashworthiness and Failure, (N.
Jones and T. Wierzbicki, Eds. (Elsevier Applied Science, London, 1993) pp. 1-54.
[3] K.-J. Bathe, Finite Element Procedures, Prentice-Hall Inc, 1996, ISBN 0-13-301458-4
[4] Z.P. Bazant, T. Belytschko, Wave propagation in strain softening bar: exact solution,
J.Engrg.Mech., ASCE, Vol.111, p.381-389, 1985
[5] Z.P. Bazant, T. Belytschko, T.P. Chang, Continuum theory for strain softening,
J.Engrg.Mech., ASCE, Vol.110, p.1666-1692, 1984
[6] T. Belytschko, K.W. Liu, B.Moran, Nonlinear Finite Element Analysis for Continua
and Structures, John Wiley and Sons, 2000, ISBN 0471987743
[7] N. Bicanic, A. Munjiza, D.R.J. Owen, N. Petrinic, From Continua to Discontinua - A
Combined Finite/Discrete Element Modelling in Civil Engineering, Developments in
Computational Techniques for Structural Engineering, Topping,B.H.V. (ed.), p.1-13,
CIVIL-COMP PRESS, Edinburgh, 1995
[8] N. Bicanic (Editor), Proceedings of ICADD-4 Fourth International Conference on
Analysis of Discontinuous Deformation, Glasgow, 2001, ISBN 0-85261-735-6
[9] J. Bonet, R.D. Wood, Nonlinear Continuum Mechanics for Finite Element Analysis,
Cambridge University Press, 1997, ISBN: 052157272X
[10] P.W. Bridgman, Studies in Large Plastic Flow and Fracture, McGrow-Hill, 1952
[11] M.A. Crisfield, Non-linear finite element analysis of solids and structures, John Wiley
and Sons Ltd, 1996, ISBN 047197059X
[12] CRAHVI – Crashworthiness of Aircraft for High Velocity Impact, The European
Commission Community Research, Competitive and Sustainable Growth Programme,
New Perspectives in Aeronautics, www.crahvi.net, http://dbs.cordis.lu/fep-
cgi/srchidadb?CALLER=PROJ_GROWTH&ACTION=D&QM_EP_RCN_A=54244
[13] L.A. Cross, S.J. Bless, A.M. Ranjendran, E.A. Strader and D.S. Dawicke, New
Technique to Investigate Necking in a Tensile Hopkinson Bar, Exp.Mech., Vol 24
(No.3), p.184-186, 1984
[14] E.D.H. Davies and S.C. Hunter, The dynamic Compression Testing of Solids by the
Method of the Split Hopkinson Pressure Bar (SHPB), J.Mech.Phys.Solids, Vol.11,
p.155-179, 1963
[15] P.A. Cundall, R.D. Hart, Numerical Modelling of Discontinua, DEM 1st U.S.
Conference, CSM Press, Golden, Colorado, 1989
36
[16] G. Dieter, Mechanical Behaviour Under Tensile and Compressive Loads, Mechanical
Testing and Evaluation, ASM Handbook, Vol.8, p.99-108, ASM International, The
Materials Information Society, Ohio, 2000
[17] U.-J. Görke, S. Kretzschmar R. Kreiβig, Analysis of inhomogeneous displacement fields
for the identification of parameters for elasto-plastic deformation laws, In Sol H,
Oomens CWJ, editors. Material Identification Using Mixed Numerical Experimental
Methods. The Netherlands: Kulwer Academic Publishing, 1997
[18] G.T. Grey III, Classic Split-Hopkinson Pressure Bar Testing, Mechanical Testing and
Evaluation; ASM Handbook, Vol. 8, p.462-476; ASM International, The Materials
Information Society; Ohio; 2000
[19] J. Harding, E.O. Wood and J.D. Campbell, Tensile Testing of Materials at Impact Rates
of Strain, J.Mech.Eng.Sci., Vol.2, p.88-96, 1960
[20] IMPACT - Improved Failure Prediction for Advanced Crashworthiness of
Transportation Vehicles, The European Commission Community Research,
Competitive and Sustainable Growth Programme, http://dbs.cordis.lu/fep-
cgi/srchidadb?CALLER=PROJ_GROWTH&ACTION=D&QM_EP_RCN_A=52091
[21] W. Johnson, Impact Strength of Materials, Edward Arnold Publishers, 1972, ISBN 0-
7131-3266-3
[22] S. Kobayashi, S.-I. Oh and T. Altan, Metal Forming and the Finite Element Method,
Oxford University Press, 1989, ISBN 0-19-504402-9
[23] E. Kuhl, E. Ramm, Simulation of strain localization with gradient enhanced damage
models, Computational Materials Science, Vol.16, p. 176-185, 1999
[24] H. Kuhn and D. Medlin (eds.), Mechanical Testing and Evaluation, ASM Handbook,
Vol. 8, ASM International, The Materials Information Society, Ohio, 2000, ISBN 0-
87170-389-0
[25] J. Lemaitre, J.L. Chaboche, Mechanics of Solid Materials; Cambridge University Press;
1994; ISBN: 0521477581
[26] J. Lemaitre, A Course on Damage Mechanics, Springer Verlag, 1992, ISBN
0387536094
[27] J. Lemaitre (Editor), Handbook of Materials Behavior Models, Volume Set, Academic
Press; 2001, ISBN 0124433413
[28] D. Macdougall, J. Harding, A constitutive relation and failure criterion for Ti6Al4V
alloy at impact rates of strain, Journal of the Mechanics and Physics of Solids, Vol.47,
p.1157-1185, 1999
[29] D. Macdougall, N. Petrinic, A local optimisation method for obtaining material model
parameters, Engineering Computations, Vol.17(2&3), p.175-189, 2000
[30] R. Mahnken, E. Stein, Parameter identification for viscoplastic models based on
analytical derivatives of a least-squares functional and stability investigations,
International Journal of Plasticity, Vol.12(4), p.451-479, 1996
[31] M.A. Meyers, Dynamic Behaviour of Materials, John Wiley & Sons Inc, 1994, ISBN
047158262X
[32] M.A. Mayers, R.W. Armstrong, H.O.K. Kirchner, (eds.), Mechanics and Materials
Fundamental and Linkages, John Wiley & Sons, Inc., ISBN 0-471-24317-5
37
[33] A. Needleman, Material rate dependence and mesh sensitivity in localization problems,
Computer Methods in Applied Mechanics and Engineering, Vol.67, p.69-85, 1988
[34] D.R.J. Owen, E. Hinton, Finite Elements in Plasticity: Theory and Practice; Pineridge
Press; 1980; ISBN: 0906674052
[35] N. Petrinic, B.C.F. Elliott, B. Hackl, H. Engl, Multi-Action Inverse Modelling of
Material Response to Impact Loading, SFB Conference on Computational Methods for
Inverse Problems, Strobl, Austria, 2002
[36] N. Petrinic, DEST - Discrete Element Simulation Tools,
www.eng.ox.ac.uk/~ftgnp/dest.html
[37] A.J. Rosakis, G. Ravichandran, Dynamic failure mechanics, International Journal of
Solids and Structures, Vol.37, p.331-348, 2000
[38] J.W. Rudnicki, J.R Rice, Conditions for localization in pressure sensitive dilatant
materials, J.Mech.Phys.Solids, Vol.23, p.371-394, 1975
[39] E.A. de Souza Neto, D.R.J. Owen and D. Peric, Computational Modelling of Problems
in Single Crystal Plasticity, European Conference on Computational Mechanics
ECCM’99, München, Germany, 1999
[40] S.M. Walley, P.D. Church, M. Furth and J.E. Field, A High-Speed Photographic Study
of the Rapid Deformation of Metal Annuli: A Comparison of Theory with Experiment,
J.Phys. (France) IV Colloq. C3 (EURODYMAT 97), Vol.7, p.317-322
[41] O.C. Zienkiewicz and R.L. Taylor, The finite element method, McGraw Hill, Vol. I.,
1989, Vol. II, (1991).
[42] J.A. Zukas, T. Nicholas, H.F. Swift, L.B. Greszczuk and D.R. Curran, Impact
Dynamics; John Wiley and Sons, ISBN 0-471-08677-0
38

Petrinic (Impact)

Uploaded by

Copyright:

Available Formats

Petrinic (Impact)

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

Petrinic (Impact)

Uploaded by

Copyright:

Available Formats

VII International Conference on Computational Plasticity

AN APPROACH TO PREDICTIVE MODELLING OF FAILURE IN

Abstract: This paper introduces an approach to predictive modelling of failure in metals

Experimental observations Modify experiments

Theoretical abstraction Modify theory

Numerical algorithms Modify algorithms

Data for numerical analysis

Figure 1: Integrated experimental-numerical approach to predictive modelling of deformation and failure

A0 a) quasi-static testing b) medium rate testing

bearing input bar

velocity seal chamber B

The interpretation of data generated by a tensile Split-Hopkinson-Bar relies upon 1D

2.1 Analysis of longitudinal elastic waves in thin rods

Hence the following D’Alambert’s solution of the 1D wave equation (2)

In Figure 6 the function

2.2 Method of characteristics

results in the following expression for stress

Substituting into (33) and (34) it follows that

which determines the magnitude of backward travelling stress wave (β wave)

and the velocity of particles

2.3 Split-Hopkinson-Bar data analysis

Figure 11: Split-Hopkinson-Bar experimental setup

and its corresponding resistance force as the following functions of time

 β xMAX , t = σ xG 2 , t + (t3 − t2 ) − σ xG1 , t + (t3 − t 2 ) − (t2 − t1 )

 β xMIN , t = −α xMIN , t − (t8 − t5 )

Figure 12: Enhanced Split-Hopkinson-Bar analysis

A new software (ELAN2) for analysis of data generated by Split-Hopkinson-Bar systems

Nominal strain rate (1/s)

T11 Onset of specimen fracture L10

Figure 15: Microscopic observation of cracking in a ductile metallic alloy

3.1 The Discrete Element Method

Figure 16: Discrete element discretisation of solid bodies

ρ& + ρ div[u& ] = 0 (55)

- constitutive relation (rate form)

- contact interface impenetrability condition

ρe& = σ : d + r − div[q ] (63)

and the second principle

and noting that

If the Helmholtz free energy is defined as the following function

and the consequent additive decomposition of the velocity gradient

from which it follows that

thus resulting in the following form of the Clausius-Duhem inequality

variables) or its complementary potential (a convex scalar valued function of thermodynamic

As a result, the first principle of thermodynamics (the conservation of energy) expressed as

σ : d el + AiV&i + ρθ s& = σ : d + r − div[q ] (90)

the temperature evolution law can be expressed as follows

as well as the calculation of a hardening-minus-recovery form of the trial back stress

and the scalar valued trial isotropic hardening

is the Von Mises effective rate of plastic deformation and

is the yield stress and

is the effective rate of deformation.

αt + ∆ t = αtTR+ ∆ t + h(θ t ) ∆ γ n (124)

]1 − min1, 1D− ∆Dt  + 3K (1 − D ) 13 tr d 

thus enabling temperature update

is the plastic rate of deformation.

m = C21 Lc + C22 (133)

Figure 17: Mesh dependency prior to damage evolution regularisation

3.2 Inverse modelling for identification of non-measurable constitutive parameters

You might also like