Nonlinear Solid Mechanics

A.H. van den Boogaard

E.S. Perdahcıoğlu
V.T. Meinders

April 2020
course: 201400042
Contents

Preface vii

1 Introduction 1
1.1 Origins of nonlinear elasticity and inelasticity . . . . . . . . . . . . . . . . . . 1
1.2 The stress–strain curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.1 Metal plasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.2 Time dependent behavior . . . . . . . . . . . . . . . . . . . . . . . . . 3

2 Vectors and tensors 5

2.1 Vectors and vector operations . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2 Tensors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.3 Gradients and divergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.4 Material time derivatives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.5 The divergence theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.6 Transformation rules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.7 Invariants, eigenvalues and eigenvectors of second order tensors . . . . . . . . 15
2.8 Tensor derivatives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3 Kinematics and Strain 19

3.1 Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.1.1 Polar decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.2 Strain definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.3 Relation between Rate of Deformation and and the material time derivative of
the Deformation Tensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.4 Invariancy and Objectivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.4.1 Some Objective rates . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.5 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.5.1 Rotation in a plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.5.2 Simple stretching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.5.3 Simple shear . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.5.4 Rigid rotation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

4 Stress and Equilibrium 35

4.1 The true stress . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
4.2 Equilibrium equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
4.2.1 Force balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
4.2.2 Moment balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

iii
iv CONTENTS

4.3 Alternative stress representations . . . . . . . . . . . . . . . . . . . . . . . . . 39

4.3.1 Stress energy rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4.3.2 The Kirchhoff stress tensor . . . . . . . . . . . . . . . . . . . . . . . . 39
4.3.3 The first Piola–Kirchhoff stress tensor . . . . . . . . . . . . . . . . . . 39
4.3.4 The second Piola–Kirchhoff stress tensor . . . . . . . . . . . . . . . . . 40

5 Hyperelastic material models 41

5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5.2 General formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
5.2.1 Isotropic material . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
5.2.2 Stress calculations based on invariants . . . . . . . . . . . . . . . . . . 43
5.3 Nearly incompressible material . . . . . . . . . . . . . . . . . . . . . . . . . . 44
5.3.1 Stress calculations based on distortional invariants . . . . . . . . . . . 46
5.4 Specific material models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
5.4.1 Volumetric part of the strain energy density . . . . . . . . . . . . . . . 47
5.4.2 Neo-Hookean material model . . . . . . . . . . . . . . . . . . . . . . . 47
5.4.3 Mooney–Rivlin material model . . . . . . . . . . . . . . . . . . . . . . 47
5.4.4 Yeoh material model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.4.5 Generalized Rivlin model . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.4.6 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

6 Plasticity theory 51
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
6.2 One-dimensional plasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
6.2.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
6.3 Plastic deformation at multi-axial stress states . . . . . . . . . . . . . . . . . 53
6.3.1 The yield surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
6.3.2 Normality rule and the convexity of the yield surface . . . . . . . . . . 54
6.4 Elastic - ideal plastic constitutive equations . . . . . . . . . . . . . . . . . . . 59
6.5 Multi-axial yield functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
6.5.1 Tresca yield criterion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
6.5.2 Von Mises yield criterion . . . . . . . . . . . . . . . . . . . . . . . . . 62
6.5.3 Comparison of the Tresca and Von Mises yield surface . . . . . . . . . 64
6.5.4 The Hill yield criterion . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
6.5.5 Non-associated flow behavior . . . . . . . . . . . . . . . . . . . . . . . 65
6.5.6 Constitutive equations for elastic-ideal plastic material, using the Von
Mises yield criterion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
6.6 Constitutive equations for a hardening material . . . . . . . . . . . . . . . . . 68
6.7 Hardening . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
6.7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
6.7.2 Isotropic hardening . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
6.7.3 Kinematic hardening . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
6.8 Elastic-viscoplastic material . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
6.8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
6.8.2 Creep in metals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
6.8.3 Visco-plastic material behavior . . . . . . . . . . . . . . . . . . . . . . 80
CONTENTS v

7 Finite Element Method 83

7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
7.2 Virtual work/weak formulation of equilibrium . . . . . . . . . . . . . . . . . . 83
7.3 Finite element discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
7.4 Iterative Newton–Raphson procedure . . . . . . . . . . . . . . . . . . . . . . . 85
7.5 Determination of the stiffness matrix . . . . . . . . . . . . . . . . . . . . . . . 88
7.5.1 Geometrically linear case . . . . . . . . . . . . . . . . . . . . . . . . . 88
7.5.2 Geometrically nonlinear case . . . . . . . . . . . . . . . . . . . . . . . 89
7.6 Stress update . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
7.6.1 Cutting-plane algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . 94
7.7 A note on Finite Element implementation . . . . . . . . . . . . . . . . . . . . 96

8 Continuum Damage Mechanics 99

8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
8.2 Uncoupled damage models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
8.3 Coupled damage models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
8.3.1 Basic concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
8.3.2 Brittle fracture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
8.3.3 Ductile damage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

9 Tensile instability in sheet metal 105

9.1 Plastic deformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
9.2 Diffuse necking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
9.3 Localized necking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
9.3.1 Localization in the drawing region . . . . . . . . . . . . . . . . . . . . 108
9.3.2 Localisation in the stretching region . . . . . . . . . . . . . . . . . . . 111

A Principal stresses and invariants 115

Preface

In the course Nonlinear Solid Mechanics, we consider two phenomena that lead to nonlin-
ear relations between forces and displacements: large displacements and nonlinear material
behavior. Both play an important role in material forming processes, in which products are
shaped by plastic deformation. The simulation of forming processes is an important applica-
tion area of the models that are presented in this course. It is the primary research area of the
chair Nonlinear Solid Mechanics. A second application where both large deformations and
nonlinear material behavior are combined is the analysis of rubber or in general of elastomeric
products. For this purpose, hyperelastic models are used, as explained in Chapter 5.
A deeper understanding of the acting process forces is necessary in order to determine the
strength of the tooling and the required press forces. Insight in the material flow is important
in order to assess whether the product will acquire the prescribed shape.
Shape deviations are caused by fracture or tearing during the manufacturing process, but
also by deformations after the forming stage as a result of residual stresses – so-called elastic
spring back.
The methods that will be used for the analysis of forming processes are based on continuum
mechanics. Microscopic material behavior and micro structures are therefore not considered
in this course.
For the analysis of deformation processes, the mathematical description of the material behav-
ior during plastic deformation, by means of constitutive equations, is very important. These
equations present a relation between stresses and deformations. The most general constitutive
equations for anisotropic elastic-plastic material are quite complex and problems with real
deformation processes can often not be solved analytically. Even if we restrict ourselves to
isotropic elastic-plastic behavior analytical solutions are hard to get. Application of numerical
solution methods is then the only feasible option.
With the resulting equations describing the process, it is possible to derive approximating
analytical solutions. With this analytical approximation, a relatively accurate estimation of
the process forces can be made. Analysis of the material flow is not very accurate. A more
precise prediction can be obtained with numerical methods. Therefore, some attention will
be given to the application of the finite element method in the analysis of forming processes.
The relevant quantities in our analysis, like stress and deformation, are mathematically de-
scribed as tensors. In the first chapters, a short introduction to some topics of tensor algebra
is given and the index notation is introduced. Thereafter, we consider the stress tensor and
rate-of-deformation tensor. Next, the constitutive equations for hyperelastic and rigid-plastic
material, the equilibrium equations and the virtual work equation are elaborated. In Chap-
ter 7, we will treat the application of the finite element method. It is followed by a chapter on
the description of damage in a material and a chapter on the stability of plastic deformation
in sheet metal forming.

vii
viii Preface

The course, and also this reader, builds upon previous courses Plasticity and Creep and Solid
Mechanics 2, by Prof. Han Huétink and Dr. Timo Meinders.
Chapter 1

Introduction

1.1 Origins of nonlinear elasticity and inelasticity

Elastic deformation is by definition reversible meaning that once the stress on the material is
removed the original undeformed shape is recovered. Even when isolated from the geometrical
effects that cause nonlinearities (i.e. large deformation stress and strain measures) elasticity
of the material can be nonlinear. The nonlinear behavior however does not contradict with
the above statement requiring reversibility of deformation.
The remaining nonlinearity in elasticity is in general caused by phenomena that take place in
the microstructure of the material. Rubber for instance is a material made of polymer chains
that can be arbitrarily oriented in the microstructure. During deformation however because of
the alignment of chains towards the loading direction the apparent material stiffness changes.
However this is completely reversible and the originial configuration is recovered when stress
is released.
Inelasticity on the other hand is defined as all types of deformation that is not recoverable
upon releasing the load. In metals this is also referred to as plasticity of the material. Plas-
ticity in metals occurs because of dislocation motion in their microstructure. Dislocations
appear naturally in all metals and drastically reduce the energy required for slipping of crys-
tallographic planes over each other.
When there is load on the material the atomic bonds in the crystals are disturbed from their
equilibrium position resulting in a stress field. When the stress resolved on certain crystallo-
graphic planes and in certain directions reach a critical value the dislocations that reside on
that plane will start moving resulting in a net irreversible deformation. The cumulative sum
of all slip activities are referred to as plastic deformation.
As the dislocations move they interact with defects in the crystal structure, such as inclusions,
grain boundaries or other dislocations and through these interactions more dislocations are
generated. The generation of more dislocations in turn make it actually harder for the move-
ment of dislocations which means that higher resolved stresses become necessary for further
plastic deformation. This process is called strain (work) hardening.
As plastic deformation progresses further in the material due to the increased average stress
as well as locally concentrated high stresses around imperfections in the microstructure at
some point it becomes energetically more favorable to initiate voids rather then dislocation
movement. These voids naturally occur and during further loading they grow in size until
they start to touch each other and coalesce. The nucleation, growth and coalescence of

1
2 Introduction

voids in ductile materials is referred to as the damage mechanism. After coalescence usually
microcracks form inside the material which lead to fracture.

1.2 The stress–strain curve

The uniaxial tensile test is a fundamental tool to obtain information of the mechanical be-
havior of materials. In this test the material is elongated with a prescribed deformation and
the corresponding load that is required is recorded. The deformation that is applied is such
that it generates a uniaxial stress state in the material. The result of the test without any
manipulation is force vs elongation. By dividing the force to original cross-sectional area and
the elongation by the original length the engineering stress (s)– engineering strain (e) curve is
obtained. As an example of an engineering stress–strain curve for a metal we consider Figure
1.1.

s s
σuts
σfrac
σy
σp

(a) e (b) e

Figure 1.1: Typical engineering stress–strain curve for (a) metals (b) rubber-like materials.

1.2.1 Metal plasticity

The material behaves almost linear elastic until σy , the yield stress, after which plastic de-
formation starts. A small deviation from linearity before yielding is usually observed and
referred to as anelasticity. The point at which deviation from linearity starts is referred to as
the proportionality limit, σp . Although it is much easier to deform the material after yielding
compared to the elastic regime, the stress still increases after σy due to strain hardening.
On the other hand, due to the lateral contraction of the specimen as the plastic deformation
proceeds the cross sectional area decreases. These two competing factors result in a maximum
on the stress-strain curve where the hardening rate equals the decrease of the force due to
reduction of area. This point on the curve is referred to as the Ultimate Tensile Strength and
denoted as σuts . At this point instability starts in the specimen which leads to necking. The
plastic deformation localizes in the neck where the plastic strain increases rapidly. Due to
the large plastic strains damage occurs in the material which finally leads to the fracture of
the specimen at σfrac.
Upon unloading the material the stress drops almost linearly with strain until it vanishes.
The small amount of nonlinearity again is described by anelasticity. When the loading is
continued however on the compressive side it is observed that for most metals the absolute
1.2 The stress–strain curve 3

value of the (re-)yield stress is lower than the value for the original loading direction. This is
referred to as the Bauschinger effect and it is depicted in Figure 1.2.

Figure 1.2: The Bauschinger effect.

1.2.2 Time dependent behavior

The mechanical behavior of most materials is influenced by the speed of loading. In rubber-
like materials this is usually described as viscoelasticity and in metals as viscoplasticity.
Generally time dependency is important due to the fact that different stress-strain curves are
obtained for the same material when the speed of loading is changed. On the other hand at
slightly elevated temperatures the effects of the time show itself as solely time dependent de-
formation such as the creep phenomenon where additional plasticity is obtained in a material
without changing the stress. Additional phenomena such as stress relaxation and recovery
also come into picture at elevated temperatures.
Creep is usually attributed to microstructural phenomena such as dislocation activity and
grain boundary sliding at constant stress. At a constant load therefore the obtained strain
can be plotted as a function of time as in Figure 1.3a.
A similar behavior is observed when the material is constrained against deformation at ele-
vated temperatures. In this case it is observed that the stress in the material gets reduced
as a function of time. This is usually referred to as stress relaxation and recovery and it is
depicted in Figure 1.3b. Relaxation is commonly seen in polymers and is attributed to vis-
coelasticity. In metals however the stress usually decreases at elevated temperatures due to
the effect of dislocation annihilation. Dislocations can move more easily at high temperatures
and in order to reduce the energy stored in the material they annihilate each other thereby
reducing the current yield stress of the material.
4 Introduction

e s

(a) (b)
t t
Figure 1.3: Observed (a) creep at constant stress and (b) relaxation and recovery at constant
strain.
Chapter 2

Vectors and tensors

In this chapter tensor notation is introduced that generalizes the concept of vectors. Although
this notation may seem cryptic at first sight, it is very strong in describing complex relations
in a concise manner. We will need this notation in the following chapters to avoid lengthy
equations for general 3-dimensional stress-states and deformations.

2.1 Vectors and vector operations

✯
a

Figure 2.1: Representation of a vector without reference to a coordinate system.

First of all, consider the vector a as drawn in Figure 2.1. We will assume that this vector
represents a physical quantity, like e.g. a displacement, a velocity or a force. The physical
meaning is independent of our (arbitrary) choice of a coordinate system and it is represented
by the arrow.
Two different notation schemes are used in this reader: index-free tensor notation and index
notation. In index-free tensor notation vectors are denoted by bold letters1 , like a. In index
notation, the coordinates in a specific vector basis are used to denote a vector. This is close
to how we expressed a vector in linear algebra, using a column of coordinates:
 
a1 

 

{a} = a2 (2.1)

 
a 
3

1
In handwriting normal and bold symbols are hard to disinguish and we will put a tilde underneath the
symbol, like a.
˜

5
6 Vectors and tensors

Figure 2.2: Vectors in a coordinate system.

We will restrict ourselves to Cartesian coordinates in R3 (Euclidean space), with unit base
vectors e1 , e2 , e3 as represented in Figure 2.2. For an arbitrary vector a with components
(a1 , a2 , a3 ):
a = a1 e 1 + a2 e 2 + a3 e 3 (2.2)
or
3
X
a= ai e i (2.3)
i=1
We see that in (2.3) the summation takes place over the index that appears twice after the
summation sign for the coordinate directions 1, 2 and 3. In the sequel we will omit the
summation sign and its bounds, and implicitly take the sum over the indices that appear
twice:
a = ai e i dummy-index (2.4)
This is called the summation convention or Einstein convention. For the double index, it is
irrelevant whether it is i, j, k or any other symbol and it is therefore also called a dummy-
index. If a double index occurs it means that the value of the index is subsequently substituted
by 1, 2 and 3 and that the sum of the results of these three substitutions is taken.
In index notation, the vector is commonly denoted just by the components ai , where it is
implicitly assumed that i ∈ {1, 2, 3} and that these are the components with respect to a
Cartesian basis {e1 , e2 , e3 }. So formally, a is a vector and ai is a set of components of a
vector and they are not equal. The relation between them is given by the three equivalent
equations (2.2), (2.3) or (2.4). In practice, we often talk about the vector a in the index-free
notation and the same vector ai in index notation.
When e is a unit vector in the direction of a then

a = kake (2.5)

where kak is the length of vector a.

Vector operations
Scalar multiplication of vector a by α results in a vector b in the same direction but with a
length of αkak:
b = αa or bi = αai (scalar multiplication) (2.6)
 2.1 Vectors and vector operations 7

✒
b
✿
a
φ

Figure 2.3: Two vectors without reference to a coordinate system.

The sum of two vectors a and b gives a vector c that can be constructed by letting vector b
start at the end of vector a, c is then connecting the start of vector a to the end of b. Vector
addition is commutative.

(sum of two vectors) c = a + b = b + a or ci = ai + bi (2.7)

The inner product or dot product or scalar product of two vectors is a scalar quantity that is
defined as
a · b = kakkbk cos(φ) (inner product of two vectors) (2.8)
where φ is the smallest angle between a and b as presented in Figure 2.3. From this definition,
it is immediately clear that the inner product of two vectors is commutative:

a·b=b·a (2.9)

and that the inner product of two orthogonal vectors is zero. The definition also implies that
the inner product of a vector with itself is the square of the length of that vector:

a · a = kak2 (2.10)

In linear algebra, a distinction was made between row and column vectors. In tensor notation,
a vector has a particular geometrical interpretation, independent of the coordinate system.
The inner product is uniquely defined for two vectors and therefore there is no need to relate
to a representation in a column or row of coordinates or to use the transpose for the first
vector.
The base vectors ei are mutually orthogonal unit vectors and hence, by using definition (2.8)
and (2.10): (
1 if i = j
ei · ej = δij = (2.11)
0 if i 6= j
where we implicitly defined the Kronecker delta symbol δij .
Using the Kronecker delta symbol and the summation convention, we can write the vector
component ai as:
ai = δij aj (2.12)
where j is a dummy index taking the values 1, 2 and 3, but only if i = j, the value of aj
contributes to the result.
The inner product of two vectors a and b can then be written as:

(inner product of two vectors) a · b = ai ei · bj ej = ai bj ei · ej = ai bj δij = ai bi (2.13)

8 Vectors and tensors

The components of a vector in a particular Cartesian basis can easily be obtained using the
inner product.
ai = a · ei (components of vector a) (2.14)
This follows directly from (2.2) and orthogonality of the base vectors.
The cross product or vector product of two vectors results in another vector:

a × b = c (cross product of two vectors) (2.15)

where the length of c is equal to the area of the parallelogram spanned by a and b:

kck = kakkbk sin(φ) (2.16)

and the direction of c is perpendicular to the plane spanned by a and b, in the direction given
by the right-hand rule. Because of this right-hand rule, the cross product in not commutative,
but:
a × b = −b × a (2.17)
The triple product of 3 vectors a, b and c, takes the inner product of a with the cross product
of b and c. Since the cross product represents the area of the parallellogram spanned by b
and c and the inner product subsequently defines the size of a in the direction of the normal
to the parallellogram, the triple product represents the signed volume V of the parallellepiped
spanned by a, b and c.
V = |a · (b × c)| (triple product of 3 vectors) (2.18)
The sign of the triple product is positive or negative depending on the orientation of the
vectors with respect to each other.

2.2 Tensors
A vector is also called a first order tensor and a scalar is sometimes called a zero-th order
tensor. We can also define higher order tensors.
First we define the dyadic product or tensor product of two vectors a and b by the dyad2 ab
implicitly by the property

ab · c = a(b · c) and c · ab = (c · a)b (2.19)

A dyad can be expressed in index notation as

ab = ai ei bj ej = ai bj (ei ej ) (dyad) (2.20)

Here, we used the property that a scalar multiplication (with bj ) can be replaced freely within
the product. However this does not hold for the dyadic product ei ej ; the order of vectors
within a dyadic product must be maintained. Although a dyadic product is often just written
as ai bj in index notation (and this is equal to bj ai since the components are scalar values), it
must be reminded that this represents coordinates for a dyad of the vectors ei and ej in that
particular order.
2
We follow the convention that there is no operator symbol between the two vectors to denote the tensor
product. In another widely used notation set this dyad is denoted by a ⊗ b.
2.2 Tensors 9

The property given in (2.19) is easily demonstrated in index notation:

     
a1 b1 a1 b2 a1 b3   c1 
 
a 1 (b · c) 

     
(ab) · c = a2 b1
 a2 b2 a2 b3  · c2 = a2 (b · c) = a(b · c)
 (2.21)

   
a3 b1 a b a b c   a (b · c) 
3 2 3 3 3 3

Notice the similarity with matrix–vector multiplication from linear algebra.

A tensor product is a linear operator:

a(αb + βc) = αab + βac (2.22)

A second order tensor A is an operator that maps an arbitrary vector b to another vector c.
It can be expressed as a linear combination of dyads and is represented in this reader by bold
letters3 , e.g. A. A second order tensor can be expressed in index notation with reference to
a Cartesian coordinate system as:
A = Aij ei ej (second order tensor A) (2.23)
and the components can be obtained from:

Aij = ei · A · ej (components of second order (2.24)

tensor A)
The components Aij can be represented in a 3 × 3 matrix.
An inner product of a second order tensor A and a vector b can be interpreted as the
summation of all inner products of the last vectors in the dyads that constitute A with b.
The result is another vector. In index notation it is represented as:

c = A · b = Aij ei ej · bk ek = Aij bk ei δjk = Aij bj ei (2.25)

It is often shortly written as ci = Aij bj .

Similarly, an inner product by a vector and a second order tensor or by two second order
tensors can be defined. The inner product always works on the vectors (eventually as part of
a dyad) directly to the left and to the right of the dot. If second (or higher) order tensors are
involved, the inner product is not commutative, e.g.: A · b 6= b · A.
The transpose of a second order tensor swaps the order of the vectors in the constituting
dyads. By using the transpose, the pre-multiplication with a vector can be written as a
classical post-multiplication:
for transpose different order
d = b · A = bi ei · Ajk ej ek = Ajk ek ej · bi ei = AT · b (2.26)

because the inner product ei ·ej results in a scalar (δij ) and can therefore be placed at another
position in the product.
The second order identity tensor 1 maps any vector onto itself:

c=1·c=c·1 (2.27)

By the definition of a second order tensor as linear combination of dyads, it also holds for
second order tensors:
A=1·A=A·1 (2.28)
3
In handwriting we will put a bar underneath the symbol, like A.
10 Vectors and tensors

The second order identity tensor is defined by:

c = ci ei = c · ei ei = c · 1 ∀c ⇒ 1 = ei ei (2.29)

The components of the second order identity tensor are given by the Kronecker delta symbol:

1 = δij ei ej = ei ei (2.30)

For 2 dyadic products the double contraction or double dot product can be defined by4 :

(double contraction) ab : cd = (a · c)(b · d) (2.31)

and likewise for 2 second order tensors: dyadic product between vector a and b

A : B = Aij ei ej : Bkl ek el = Aij Bkl (ei · ek )(ej · el ) = Aij Bij (2.32)

The colon represents a summation over both indices (double contraction) and the result is
a scalar. Because the dot product of two vectors is independent of the coordinate system,
it also holds that the double contraction of two second order tensors is independent of the
coordinate system.
Next, the trace of a second order tensor can be defined as:

trA = Aii (2.33)

which equals the sum of the diagonal terms of the components of A when written as a matrix.
A double contraction of a second order tensor A with the second order identity tensor also
yields the trace of A:

1 : A = δij ei ej : Akl ek el = Akl ei ei : ek el = Akl δik δil = Aii = trA (2.34)

Analogous to the previous discussion we can also define third and fourth order tensors as
linear combinations of tryads of 3, repectively tetrads of 4 vectors. Third and fourth order
tensors will also be denoted by bold symbols5 . The order of the tensor must be derived from
the context. A third order tensor is a linear operator that maps a vector onto a second order
tensor or a second order tensor onto a vector.
A=D·a
(2.35)
Aij = Dijk ak

a=D:A
(2.36)
ai = Dijk Ajk
For a fourth order tensor similar operations hold, e.g.:

A=E:B (2.37)

Aij = Eijkl Bkl (2.38)

4
Be aware that in some ‘schools’ the double contraction is defined by ab : cd = (a · d)(b · c).
5
In handwriting we will put a bar and tilde underneath third order tensors like a and and a double bar
underneath fourth order tensors, like E ˜
2.3 Gradients and divergence 11

The fourth order identity tensor is defined by the condition that for any second order tensor
the following relation must hold:
A = 41 : A (2.39)
The fourth order identity tensor can be written as:
4
1 = δik δjl ei ej ek el (2.40)

Note that there is a large difference between the fourth order identity tensor and the dyadic
product of two second order identity tensors. This can easily be shown when contracting
these tensors with an arbitrary second order tensor A:
4
1 : A = δik δjl ei ej ek el : Apq ep eq = Apq ei ej ei ej : ep eq
= Apq δip δjq ei ej = Apq ep eq = A (2.41)
11 : A = δij δkl ei ej ek el : Apq ep eq = Apq ei ei ek ek : ep eq
= Apq δkp δkq ei ei = Akk ei ei = Akk 1 (2.42)

Note that the latter equation yields an identity tensor in which the diagonal terms are mul-
tiplied by value of the trace of the second order tensor A.

2.3 Gradients and divergence

The partial derivatives of a scalar function f (x) with respect to the coordinates form the
components of the gradient vector ∇f :
∂f
∇f = ei = f,i ei (2.43)
∂xi
In index notation, the partial derivative with respect to direction xi is represented by a
comma, followed by i. The symbol ∇ (nabla) represents an operator in three dimensions and
is therefore presented as a vector, using a bold symbol:
∂
∇= ei (2.44)
∂xi
Partial derivatives of vectors and tensors yield tensors of one order higher. The partial
derivatives of the components of vector a are:
∂ai
= ai,j
∂xj

they form the components of the post-gradient (second order tensor):

∂ ∂ai
a∇ = ai ei ej = ei ej = ai,j ei ej (post-gradient) (2.45)
∂xj ∂xj

Notice that the differentiation ∂/∂xj operates on the vector to the left of the operator, which
cannot easily be seen in the first part of this equation. In tensor notation we also know the
pre-gradient:
∂ ∂aj
∇a = e i aj e j = ei ej = aj,i ei ej (pre-gradient) (2.46)
∂xi ∂xi
12 Vectors and tensors

The pre-gradient is the transposed (conjugated) of the post-gradient. The gradient of a second
order tensor yields a third order tensor:

∂Aij
A∇ = ei ej ek = Aij,k ei ej ek (2.47)
∂xk

The divergence of a vector a yields a scalar (dot product of the gradient operator with a
vector):
∂ ∂
div a = ∇ · a = e i · aj e j = aj ei · ej = ai,i (2.48)
∂xi ∂xi

The divergence of a second order tensor yields a vector:

A · ∇ = Aij,k ei ej · ek = Aij,j ei (2.49)

2.4 Material time derivatives

The material time derivative gives the change of a property that is connected to a material
particle. In general, a property G is a function of time and place.

G = G(t, xi ) (2.50)

The position of the material particle is a function of time: xi = xi (t), therefore G = G(t, xi (t)).
The material derivative—denoted by a superimposed dot—is given by:

dG ∂G ∂G dxi
= Ġ = +
dt ∂t ∂xi dt

or, with the velocity vi = dxi / dt:

∂G ∂G
Ġ = + vi (2.51)
∂t ∂xi

For tensors it holds that

∂Gij ∂Gij
Ġij = + vk
∂t ∂xk
or
∂G ←
−
Ġ = + G∇ · v (2.52)
∂t
The arrow above the gradient symbol points to the quantity of which the derivative must be
taken (it is omitted if this cannot be mistaken). In index notation that is directly visible by
the summation over index k.
Note:
→
− ←
−
v · ∇G = G∇ · v (2.53)
2.5 The divergence theorem 13

The proof of this equation is given for the case that G is a vector (for simplicity):
 
f f f
o  1,1 2,1 3,1 
→
− n
v · ∇f = vi fj,i = v1 v2 v3 ·  f1,2 f2,2 f3,2 

f1,3 f2,3 f3,3
 
v1 f1,1 + v2 f1,2 + v3 f1,3
 
= v 1 f 2,1 + v 2 f 2,2 + v 3 f 2,3

 (2.54)
v1 f3,1 + v2 f3,2 + v3 f3,3
   
f1,1 f1,2 f1,3  v1  
    ←
−
= f2,1 f2,2 f2,3  · v2 = fi,j vj = f ∇ · v
 

 
f f f v  
3,1 3,2 3,3 3

2.5 The divergence theorem

The divergence theorem (or Gauss’s theorem or Ostrogradsky’s theorem) gives the relation
between the integral of the divergence of a vector field over a volume V and the integral of
the flux of that vector field over the boundary S enclosing that volume:
Z Z
(v · n) dS = v · ∇ dV (2.55)
S V

or in components: Z Z
vi ni dS = vi,i dV (2.56)
S V

It also applies to higher order tensors:

Z Z
(A · n) dS = A · ∇ dV (2.57)
S V

or in components: Z Z
Aij nj dS = Aij,j dV (2.58)
S V

2.6 Transformation rules

Imagine two sets of Cartesian base vectors e1 , e2 , e3 and e∗1 , e∗2 , e∗3 (Figure 2.4). In short, we
can call ei the original or reference coordinate system and e∗i the new or current coordinate
system. The relation between the two is fully determined by the angles between the original
and the new base vectors or equivalently between the cosines of these angles. Because the
base vectors are normalized, the cosine of the angle between two base vectors Qij = cos(ei , e∗j )
is equal to the inner product
Qij = ei · e∗j (2.59)
14 Vectors and tensors

Figure 2.4: Coordinate frames.

Since a = ai ei = (a · ei )ei we can write by substituting a = e∗j :

e∗j = (ei · e∗j )ei = Qij ei (2.60)

An orthogonal tensor Q can be defined with components Qij as

Q = Qij ei ej (2.61)

It then holds that

Q · ek = Qij ei ej · ek = Qik ei = e∗k (2.62)
So Q rotates the original base vector ek to the new base vector e∗k .
By definition of Qij , Q is orthogonal and the inverse is equal to the transpose:

QTij = Q−1 ∗
ij = Qji = ej · ei (2.63)

For vector a:
a = ai ei = a∗k e∗k = a∗k Qik ei
hence
ai = Qik a∗k (2.64)
Multiplication of (2.64) with Q−1
ji gives:

∗ ∗ ∗
Q−1
ji ai = Qij ai = Qji Qik ak = δjk ak = aj
−1

take again transpose, transpose vanishes

Hence:
T
a∗j = Q−1
ji ai = Qji ai = ai Qij (2.65)
Therefore, the components of a vector transform in the inverse way as the base vectors.
This can be understood by the fact that the components of a rotated vector in an unrotated
reference frame transform in the same way as the components of an unrotated vector in a
reference frame that is rotated in the opposite direction.
The transformation of components of a second order tensor can be derived from the operation
on a vector. Assume

b=A·a ⇐⇒ bi = Aij aj and b∗i = A∗ij a∗j (2.66)

2.7 Invariants, eigenvalues and eigenvectors of second order tensors 15

b∗i = Qki bk = Qki Akl al = QTik Akl Qlj a∗j (2.67)

because Equations (2.66) and (2.67) apply for arbitrary a∗j :

A∗ij = QTik Akl Qlj = Akl Qki Qlj (2.68)

So, for tensor components, the same transformation rule holds as for vector components, but
now applied twice, once for every index.
A reference frame transformation for third order tensors can be written as:
∗
Bijk = QTip QTjq QTkr Bpqr = Bpqr Qpi Qqj Qrk (2.69)

and for a fourth order tensor:

∗
Cijkl = QTip QTjq QTkr QTls Cpqrs = Cpqrs Qpi Qqj Qrk Qsl (2.70)

If for a product of second order tensors holds:

Pij = Aik Bkj (2.71)

then also:
Pij∗ = A∗ik Bkj
∗
(2.72)
Proof:
∗
A∗ik Bkj = QTil QTkm Alm QTkn QTjp Bnp (2.73)
QTkm QTkn = Qmk Q−1
kn = δmn (2.74)
therefore:
A∗ik Bkj
∗
= QTil QTjp Aln Bnp = QTil QTjp Plp (2.75)
We also have the relation:
Pij∗ = QTil QTjp Plp (2.76)
hence:
Pij∗ = A∗ik Bkj
∗
(2.77)
In general, relations between tensor components remain valid with respect to an arbitrary
base. The transformation rules for tensor components will hardly be used in this course,
because we will use only one coordinate system.

2.7 Invariants, eigenvalues and eigenvectors of second order

tensors
Components of second order tensors in general depend on the coordinate system. Some scalar
functions of second order tensors do not depend on the coordinate system and are therefore
called invariants. For a second order tensor in three dimensions, three independent invariants
can be defined. Notice that any linear combination of invariants is again an invariant.
The first invariant of a tensor A is the trace of A as defined in (2.33):

I1 (A) = tr(A) = A : 1 (2.78)

16 Vectors and tensors

It was already observed below Equation (2.32) that the double dot product of two second
order tensors is does not depend on the coordinate system and therefore also the trace of a
second order tensor is invariant.
The second invariant of a tensor A is defined as:
1
tr2 (A) − tr(A2 )
 
I2 (A) = 2 (2.79)

where A2 = A · A. Since the dot product of two second order tensors is fully described by
the vectors constituting the dyads and the trace of any second order tensor is invariant, the
result is also invariant. In some books, the second invariant is defined as tr(A2 ) but then the
equations below where I2 is used are not valid.
The third invariant of a tensor A, also called the determinant of A, is defined as:

(A · c1 ) · (A · c2 ) × (A · c3 )
I3 (A) = det(A) = (2.80)
c1 · c2 × c3

where {c1 , c2 , c3 } represents a vector basis. If we consider the definition and interpretation of
the triple product in (2.18) it is clear that the determinant represents the ratio of volume of
the parallellepiped spanned by A · c1 , A · c2 and A · c3 and the volume of the parallellepiped
spanned by c1 , c2 and c3 and this ratio is independent of the coordinate system.
The eigenvalues λ and the eigenvectors n of a tensor A are defined as those scalars and
vectors that satisfy:
A · n = λn with n 6= 0 (2.81)
It can also be written as (A − λ1) · n = 0. A nontrivial solution n 6= 0 is obtained when the
part between brackets is singular, or:

det(A − λ1) = 0 (2.82)

Elaboration of this equation leads to the characteristic equation:

λ3 − I1 (A)λ2 + I2 (A)λ − I3 (A) = 0 (2.83)

This equation has three roots λ1 , λ2 and λ3 . Because the eigenvalues do not depend on the
coordinate system, they are also invariant.
According to the Cayley–Hamilton theorem a tensor should satisfy its own characteristic
equation, so that (2.83) becomes

A3 − I 1 A2 + I 2 A − I 3 1 = 0 (2.84)

Multiplying this equation with A−1 gives an another useful expression:

A2 − I1 A + I2 1 − I3 A−1 = 0 (2.85)

Taking the trace of the resulting equation, the following relation between the second and third
invariant can be derived:
I2
tr(A−1 ) = (2.86)
I3
which shows that tr(A−1 ) can also be used as an invariant.
2.8 Tensor derivatives 17

2.8 Tensor derivatives

The time derivative of an inverse tensor A can be obtained by considering that A−1 · A = 1
is constant and therefore the time derivative is zero:
d dA−1 dA
(A−1 · A) = ·A+A· =0 (2.87)
dt dt dt
and therefore
dA−1
= −A−1 · Ȧ · A−1 (2.88)
dt
Three other useful relations give the derivatives of the invariants of a tensor with respect to
the tensor itself. They are given here without proof.

∂I1 (A)
=1 (2.89)
∂A
∂I2 (A)
= I1 1 − AT (2.90)
∂A
∂I3 (A)
= det(A)A−T (2.91)
∂A
Chapter 3

Kinematics and Strain

3.1 Kinematics
The kinematics of a deformable body concerns the motion of a material and coordinate system
from an initial state (reference state) to the final situation, see Figure 3.1. In the reference
configuration a material particle is located at position X, the current (spatial) position of the
material particle is defined as x. Several definitions of motion can be used. In the reference
description the only independent variables are the initial position X of a material particle in
an arbitrary reference configuration and time t. In this definition the observer can be thought
of as moving with the material particle. A special case is the Lagrangian description in which
the reference configuration is chosen to be the initial configuration at time t = 0. Another
description of motion is the Eulerian description. Here the only independent variables are the
current location x of a material particle and time t. In this definition the observer is located
at a fixed point in space.

Figure 3.1: Motion of a material point.

The Lagrangian location vector X determines the initial (reference) position of a particle
within a body at t = 0. In Figure 3.1 the motion history of a particle X is indicated by the
trajectory. The motion of each particle can be described by the Eulerian location vector x
at time t, which is a function of the Lagrangian location vector X (the location in the body)
and the time t:
x = x(X, t) (Langrangian description) (3.1)

19
 20 Kinematics and Strain

By means of these coordinates a motion is described by which the particle with the material
coordinates X is transferred to location x at t > 0. The inverse of Equation (3.1) reads:

X = X(x, t) (3.2)

which can be understood by considering: an observer situated at spatial position x watching

the particle with material coordinates X pass by. It should be stressed that at time t = 0
the material and spatial coordinates are identical: x(X, 0) = X(x, 0). Application of the
concept of motion in evolving processes necessitates the introduction of time derivatives of
field and state variables. The material time derivative of a certain field or state variable Θ in
the Lagrangian description is defined for a material particle as:

(Langrangian material time derivative) Θ̇(X, t) = ∂Θ(X, t) = dΘ = ∂Θ

(3.3)
∂t dt ∂t X

In this Lagrangian or material description the independent variable is the reference position
X. In the Eulerian or spatial description the only independent variable is the fixed position
in space. The material time derivative for a material particle at spatial position x is derived
to be:

(Eulerian material Θ̇(X, t) = ∂Θ(X, t) = d θ(x, t) = ∂θ(x, t) + ∂θ(x, t) · ∂x = ∂θ(x, t) + θ ←

−
∇·v (3.4)
time derivative) ∂t dt ∂t ∂x ∂t ∂t
←
−
where θ is a certain field or state variable in the Eulerian description, ∇ is the post-gradient
and v is the spatial velocity.
∂x(X, t) dx ∂x
v = ẋ = = = (3.5)
∂t dt ∂t X

The material time derivative is split into two parts. The first part of the right hand side of
Equation (3.4) is the spatial time derivative which is viewed by an observer in a fixed position.
The second term is the convective part; the value of θ at position x changes because another
particle with a different value of θ passes by. The split into a spatial part and a convective
part is necessary if θ is only known in Eulerian coordinates.
In order to describe the deformation of a body, a deformation gradient F is defined. The
deformation gradient F relates the spatial configuration to the reference configuration for a
material particle:
∂x ←
−
F(X, t) = = x∇ 0 (3.6)
∂X t
←−
where ∇0 is the post-gradient with respect to the initial coordinates. The deformation gra-
dient F can be considered as an operator which maps the initial configuration to the current
configuration:
dx = F · dX (3.7)
in which dx and dX represent the same line element in the deformed and undeformed state,
respectively, see Figure 3.2. Generally, F is not symmetric. The displacement u of a particle
is defined as the difference between the initial and the current position (or vice versa) and
reads in Lagrangian coordinates:

u(X, t) = x(X, t) − X(x, 0) = x(X, t) − x(X, 0) (displacement) (3.8)

 3.1 Kinematics 21

Figure 3.2: Change of infinitesimal line.

In Eulerian coordinates such a difference reads:

u(x, t) = x(X, 0) − X(x, t) = X(x, 0) − X(x, t) (3.9)

The latter can be understood if it is remembered that at one certain position in space, x,
the difference is calculated between the material coordinate of the particle which was initially
at a certain position and which is present at the current time t. Therefore, Equation (3.9)
does not represent the displacement of one specific particle. The gradient of the displacement
reads, in Lagrangian and Eulerian coordinates, respectively:
←
− ∂u ∂x
(Langrangian displacement gradient) u∇0 = = −1 = F−1 (3.10)
∂X ∂X
←− ∂u ∂X
u∇ = =1− = 1 − F−1 (3.11)
∂x ∂x
where 1 is the second order unity tensor.
The determinant of the deformation gradient J represents a change in volume:
dV
J = det F = (3.12)
dV0

3.1.1 Polar decomposition

The deformation tensor can be decomposed in a rotation and a stretch tensor:

F=R·U =V·R (3.13)

in which R is the rotation tensor (which is orthogonal: R−1 = RT ) and U and V the right
and the left stretch tensor, respectively. U as well as V are symmetric and positive definite
tensors. Both decompositions in the latter equation can be considered as a split of F into a
pure deformation part, U or V, and a rigid body rotation R. For the first decomposition in
Equation (3.13) the map consists of a pure deformation, after which a rotation is applied to
the line segment in Figure 3.2. For the second decomposition this sequence is reversed. Both
decompositions are visualized in Figure 3.3.
The right and left stretch tensors are both symmetric and hence have mutually perpendicular
eigenvectors. In the direction of these eigenvectors, the material stretches by a factor, equiva-
lent to the corresponding eigenvalue, respectively λ1 , λ2 and λ3 . In these directions, no shear
deformation is observed. Rotation of the material is represented in R and not in the stretch
tensor. In the undeformed state, the stretch tensors U = V = 1 and the principal stretches
λ1 = λ2 = λ3 = 1.
22 Kinematics and Strain

Figure 3.3: Polar decomposition of the deformation gradient.

3.2 Strain definitions

In order to be able to describe the behavior of materials it is necessary to have some measure
for the amount of deformation, i.e. the strain of the material. A well known strain measure is
based on the classical linear strain theory. However, this strain measure becomes inaccurate
in case of rotations. This will be shown in the last example of this chapter.
Therefore another strain measure is necessary.Some function of the deformation gradient is
utilized for this purpose. Since the deformation tensor is not symmetric two different strain
measures can be defined by multiplying the deformation tensors. The right and the left
stretch tensors can also be used to define these measures, respectively the right and the left
Cauchy–Green tensor :
(Right Cauchy-Green tensor) FT · F = U · RT · R · U = U · U = C (3.14)
(Left Cauchy-Green tensor) F · FT = V · R · RT · V = V · V = B (3.15)
The right and left Cauchy–Green tensors are called metric tensors because they describe the
deformed metric of a body. The right and the left Cauchy–Green tensors are indeed a measure
for the strain since they map the length of a line element in the initial configuration onto the
length of the same element in the current state Figure 3.2. For example, by using Equation
(3.14) and the definition of F it follows that the right Cauchy–Green tensor can be used to
express the square of the current arc length in the Lagrangian coordinates:

ds2 = dx · dx = (F · dX) · (F · dX) =∗ (FT · F) : dX dX = C : dX dX (3.16)

scalar

Proof of *:

(F · dX) · (F · dX) = Fij dXj Fik dXk = Fij Fik dXj dXk = FjiT Fik dXj dXk =
(FT · F) : dX dX (3.17)

The other way around, the square of the initial length can be expressed in terms of the
Eulerian coordinates:

ds20 = dX · dX = (F−1 · dx) · (F−1 · dx) = (F−T · F−1 ) : dx dx (3.18)

T −1
=∗ (F · F ) : dx dx = B−1
: dx dx
3.2 Strain definitions 23

Proof of *:
Assume the following two arbitrary relations: u = A · v and w = B · u. Substitution of the
second equation in the first equation and rearranging gives:

v = (B · A)−1 · w (3.19)

Both equations can also be written as: v = A−1 · u and u = B−1 · w. Again, substitution of
the second equation in the first equation and rearranging gives:

v = A−1 · B−1 · w (3.20)

Since Equations (3.19) and (3.20) have to hold for arbitrary w, it is shown that:

A−1 · B−1 = (B · A)−1 q.e.d. (3.21)

Note that the left and the right Cauchy–Green tensors equal the unity tensor in case of no
deformation. In order to overcome this inconvenience, strains are generally defined in terms
of relative changes of the lengths of arbitrary line elements:

ds2 − ds20 = dx · dx − dX · dX (3.22)

and not in terms of relative lengths themselves. For pure rigid body motions the strain
according to this measure then vanishes. Again, this measure can be evaluated in Lagrangian
and Eulerian coordinates. Using the expression for ds2 , Equation (3.16) and recognizing that
ds20 = 1 : dX dX, the deformation measure, Equation (3.22), can be written in Lagrangian
coordinates as:
ds2 − ds20 = (C − 1) : dX dX = 2E : dX dX (3.23)
in which E is the Lagrangian strain tensor, or the Green–Lagrange strain tensor:

E = 21 (C − 1) (3.24)

This tensor vanishes in case of rigid body motions. In terms of the displacements, the Green–
Lagrange strain tensor can be written as (using Equation (3.10)):
h−→ ←
− i
E = 12 (C − 1) = 21 (FT · F − 1) = 12 (∇ 0 u + 1) · (u∇0 + 1) − 1 (3.25)
h ←− →
− →
− ←− i
= 12 u∇0 + ∇0 u + (∇0 u) · (u∇0 )

in which the geometrically non-linear term can be recognized. Analogously, by recognizing

that ds2 = 1 : dx dx , the deformation measure, Equation (3.22) can be written in Eulerian
coordinates as:
ds2 − ds20 = (1 − B−1 ) : dx dx = 2e : dx dx (3.26)
in which e is the Eulerian strain tensor, or the Euler–Almansi strain tensor:

e = 21 (1 − B−1 ) (3.27)
24 Kinematics and Strain

Like the Green–Lagrange strain tensor, the Euler–Almansi strain tensor vanishes in case of no
deformation and it is geometrically nonlinear. In terms of displacements, the Euler–Almansi
strain tensor can be written as:
h →
− ←
− i
e = 21 (1 − B−1 ) = 12 (1 − F−T · F−1 ) = 12 1 − (1 + ∇u) · (1 + ∇u) (3.28)
h ←− − → →
− ←− i
= 21 u∇ + ∇u − (∇u) · (u∇)

Note that for the Euler–Almansi strain the gradients are defined with respect to the cur-
rent state, contrary to the Green–Lagrange strain. For convenience the strain definitions in
Lagrangian and Eulerian coordinates are summarized in Table 3.1. For small displacements

Table 3.1: Metric and strain tensor definitions.

Lagrange Euler
metric tensor undeformed 1 B−1 = (F · FT )−1
deformed C = FT · F 1
strain tensor E = 12 (C − 1) e = 21 (1 − B−1 )

the geometric nonlinear terms can be neglected. It can also be proven that the gradient
with respect to the initial coordinates almost equals the gradient with respect to the current
coordinates ( ∇0 ≈ ∇ ) in case of small deformations:
 
∂u ∂u ∂x ∂u ∂u ∂u ∂u ∂u ∂u
= · = · 1+ = + · ≈ (3.29)
∂X ∂x ∂X ∂x ∂X ∂x ∂x ∂X ∂x

This means that for small deformations both the Green–Lagrange and the Euler–Almansi
strain tensors are identical to the classical linear strain tensor ε :
 ←− → 
−  ←− − → 
ε = 12 u∇0 + ∇0 u ≈ 21 u∇ + ∇u ≈ E ≈ e (3.30)

For numerical incremental calculations of plastic deformation processes, time derivatives of

process quantities are needed. Therefore, descriptions of the time derivatives of the strains
have to be determined. For these time derivatives it is convenient to define the spatial velocity
gradient L:
∂v ←−
L= = v∇ (3.31)
∂x
Besides, L can be explicitly expressed as a function of the deformation gradient:

∂v ∂v ∂X
L= = · = Ḟ · F−1 (3.32)
∂x ∂X ∂x
Like all second order tensors, the velocity gradient can be decomposed into a symmetric and
a skew-symmetric part:
L = D+W (3.33)
3.2 Strain definitions 25

where D is the symmetric part, called the rate of deformation tensor :

 ←− − → 
D = 21 v ∇ + ∇v = 21 (L + LT ) (3.34)

and W is the skew-symmetric part, called the spin or vorticity tensor :

 ← − − → 
W = 12 v ∇ − ∇v = 21 (L − LT ) (3.35)

The Green–Lagrange strain tensor and the Euler–Almansi strain tensor are frequently used in
engineering mechanics. A relation between the rate of deformation tensor D and the Euler–
Almansi strain tensor can be found by taking the time-derivatives of Equations (3.22) and
(3.26):
d
( dx · dx − dX · dX) = 2 dv· dx = 2L : dx dx = 2(D+W) : dx dx = 2D : dx dx (3.36)
dt
(Try to proof yourself that W : dx dx = 0 by writing this tensor product in components)
d d
(2e : dx dx) = 2 ( dxi eij dxj )
dt dt
= 2 ( dvi eij dxj + dxi ėij dxj + dxi eij dvj )
= 2 (Lik dxk eij dxj + dxi ėij dxj + dxi eij Ljk dxk )
= 2 (ėij dxi dxj + LTki eij dxk dxj + eij Ljk dxi dxk )
= 2 (ė + LT · e + e · L) : dx dx (3.37)

Combining both equations yields:

D = ė + LT · e + e · L (3.38)

This expression is also known as the Cotter–Rivlin derivative of e or the lower convective
time derivative of e. The rate of deformation tensor D relates to the Green–Lagrange strain
tensor via the time derivatives of Equations (3.22) (and thus Equation (3.36)) and (3.23):

Ė = FT · D · F (3.39)

The derivation of this relation is seen easily when writing Equation (3.36) in original coordi-
nates.
Another, very commonly used, strain measure is the logarithmic strain (also known as the true
or natural strain). In case of proportional loading (i.e. principal directions of deformation do
not change), the rate of deformation tensor relates to the logarithmic strain as follows. The
starting point is Equation (3.36):
d d d 2
( dx · dx − dX · dX) = ( dx · dx) = ds = 2D : dx dx ⇐⇒
dt dt dt
2 ds dṡ = 2D : dx dx (3.40)

The right hand side of Equation (3.40) can be written as a function of ds in case of pro-
portional loading only. We consider the deformation of a tensile test. Assume that the first
principal axis is parallel to the tensile direction. Then the right hand side of (3.40) reads:

D : dx dx = d11 dx211 = d11 ds2 (3.41)

26 Kinematics and Strain

Combination of Equations (3.40) and (3.41) then gives:

dṡ
d11 = (3.42)
ds
The logarithmic strain is defined as:

Zds1 Zt1 ds Zt1

1 dt dt = dṡ
ǫl = d( ds) = dt (3.43)
ds ds ds
ds0 t0 t0

The logarithmic strain, also called true or natural strain, can be considered as the integral of
the length changes,
l    
dl∗ l ∆l
Z
ǫl = = ln = ln 1 + (3.44)
l0 l∗ l0 l0
Finally, Equation (3.42) and (3.43) define the relation between the logarithmic strain and the
rate of deformation tensor for proportional loading:

d11 = ǫ̇l (3.45)

3.3 Relation between Rate of Deformation and and the mate-

rial time derivative of the Deformation Tensor
The (rate of) the deformation tensor yields:
 
d ∂x ∂v ∂v ∂x
Ḟ = = = · =L·F (3.46)
dt ∂X ∂X ∂x ∂X
The velocity gradient can then be expressed as:

L = Ḟ · F−1 (3.47)

With Ḟ = Ṙ · U + R · U̇ and F−1 = U−1 · RT the expression for the velocity gradient results
into
L = Ṙ · U · U−1 · RT + R · U̇ · U−1 · RT
or
L = Ṙ · RT + R · U̇ · U−1 · RT (3.48)
but also:
L=W+D
In the case of rigid rotation only: L = Ṙ · RT = W

From this observation one may conclude that: D = R · U̇ · U−1 · RT

This conclusions is WRONG:
D 6= R · U̇ · U−1 · RT (3.49)
In the case of non proportional deformation; U̇ and U−1 have different principal directions
and hence U̇ · U−1 is NOT symmetric.
3.4 Invariancy and Objectivity 27

Even after continued deformation without additional rotation/spin, the rotation tensor may
change whereas the material is not rotating: W = 0 but Ṙ 6= 0 and L = D = Ṙ · RT + R ·
U̇ · U−1 · RT

Conclusion: There exists no simple relation between the multiplicative split of the deformation
tensor F = R · U and the additive split of the velocity gradient L = W + D.

3.4 Invariancy and Objectivity

Consider a superimposed rigid rotation Q (Q−1 = QT ):
dx∗ = Q · dx = Q · F · dX = F∗ · dX (3.50)
The right Cauchy–Green tensor remains unchanged:
C∗ = F∗T · F∗ = FT · QT · Q · F = FT · F = C (3.51)
A tensor that remains unchanged is called INVARIANT

The left Cauchy–Green tensor transforms according:

B∗ = F∗ · F∗T = Q · F · FT · QT = Q · B · QT (3.52)
A tensor that transforms similar to B is called OBJECTIVE

Invariant tensors are: C , U , E , Ċ , Ė , (does not change after rotation)

Objective tensors are: B , V , e , D (DOES change after rotation)
Neither invariant nor objective are: F , R , L , W , ė

As stated before, there exist no simple relation between the multiplicative polar decomposition
F=R·U
and the additive split
L=W+D

Note: The rotation tensor R may change even in the case of zero spin. In the case of non
proportional deformation; U̇ and U−1 have different principal directions and hence U̇ · U−1
is NOT symmetric.

3.4.1 Some Objective rates

The Cotter–Rivlin or the lower convective time derivative:
△
e = ė + e · L + LT · e = D (3.53)
The Jaumann derivative
◦
e = ė − W · e + e · W (3.54)
The Oldroyd or upper convective time derivative:
▽
e = ė − e · LT − L · e (3.55)
28 Kinematics and Strain

3.5 Examples
The theory in this chapter will be clarified with three 2-dimensional examples. All examples
are basic deformation modes. The first example is a simple stretching deformation, the
second example will be a simple shear test and the last example comprises a rigid rotation.
All examples make use of a constant strain rate α and a time t.
First, before the examples, the tensors R and U are determined for a 2-dimensional defor-
mation.

3.5.1 Rotation in a plane

The deformation gradient can be decomposed via polar decomposition. As already mentioned,
this decomposition is unique. For the 3-dimensional case, this decomposition is quite complex.
In two dimensions, the rotation tensor R can be easily expressed by the rotation angle β along
the z-axis: " #
cos β − sin β
R= (3.56)
sin β cos β
The right stretch tensor can be written as:
" #" #
T
cos β sin β F11 F12
U=R −1
·F=R ·F= (3.57)
− sin β cos β F21 F22

Knowing that U is symmetric, the angle β can be determined from the condition U12 = U21

F12 cos β + F22 sin β = −F11 sin β + F21 cos β (3.58)

Which yields after some manipulation:

F12 − F21
sin β = − p (3.59)
(F12 − F21 )2 + (F11 + F22 )2
F11 + F22
cos β = p (3.60)
(F12 − F21 )2 + (F11 + F22 )2

Proof:
If U is symmetric then U12 = U21 so:

− sin βF11 + cos βF21 = cos βF12 + sin βF22 ⇔

− sin β(F11 + F22 ) = cos β(F12 − F21 ) (3.61)

The quadratic form of the above equation reads:

sin2 β(F11 + F22 )2 = cos2 β(F12 − F21 )2 ⇔

(1 − cos2 β)(F11 + F22 )2 = cos2 β(F12 − F21 )2 ⇔
(F11 + F22 )2 = cos2 β (F11 + F22 )2 + (F12 − F21 )2


(3.62)

The latter equation gives an expression for the cosine of β

3.5 Examples 29

Figure 3.4: Deformation mode for the simple stretching example

3.5.2 Simple stretching

The displacement function can be written as:

x = X + αtX = (1 + αt)X (3.63)

y = Y

The deformation gradient for this deformation pattern is:

∂x ∂x
   
∂x 1 + αt 0
F= =  ∂X ∂Y  =   (3.64)
∂X ∂y ∂y 0 1
∂X ∂Y
The right and left Cauchy–Green stretch tensors now become:
" #
(1 + αt)2 0
C = FT · F = (3.65)
0 1
" #
(1 + αt)2 0
B = F · FT = (3.66)
0 1
The Green–Lagrange strain tensor relates to the right Cauchy–Green stretch tensor via:
" #
α 2 t2 + 2αt 0
E = 12 (C − 1) = 12 (3.67)
0 0

Note that the same expression can be found when using Equation (3.26).
The Euler–Almansi strain tensor relates to the left Cauchy–Green stretch tensor via:
" #
(α2 t2 + 2αt)(1 + αt)−2 0
e = 21 (1 − B−1 ) = 12 (3.68)
0 0

Note that the same expression can be found when using Equation (3.29).
The classical linear strain tensor reads:
" #
 ←− →
−  αt 0
ε = 21 u∇0 + ∇0 u = 12 ((F − 1) + (FT − 1)) = (3.69)
0 0

Figure 3.5 shows the xx-components of the three strain tensors as a function of the strain
(αt). The velocity field for this problem is:
30 Kinematics and Strain

0.8

0.6

xx-component strain tensor

0.4

0.2

0
-1 -0.5 0 0.5 1
-0.2

-0.4
E
-0.6
e
-0.8 epsilon

-1
engineering strain

Figure 3.5: Components of different strain measures for stretching

vx = αX (3.70)
vy = 0

The velocity gradient reads:

∂v ∂v
L = = · F−1 (3.71)
∂x ∂X
" #" # " #
α 0 (1 + αt)−1 0 α(1 + αt)−1 0
= =
0 0 0 1 0 0

which can be split into a symmetric and skew-symmetric part, respectively:

" #
α(1 + αt)−1 0
D = 21 (L + LT ) = (3.72)
0 0

" #
1 T
0 0
W= 2 (L −L )= (3.73)
0 0

The time derivatives of the Green–Lagrange strain tensor and the Euler–Almansi strain tensor
are:

Ė = FT · D · F (3.74)
" #" #" #
(1 + αt) 0 α(1 + αt)−1 0 (1 + αt) 0
=
0 1 0 0 0 1
" #
α(1 + αt) 0
=
0 0
3.5 Examples 31

ė = D − e · L − LT · e (3.75)
" # " #" #
α(1 + αt)−1 0 (α2 t2 + 2αt)(1 + αt)−2 0 α(1 + αt)−1 0
= −
0 0 0 0 0 0
" #
α(1 + αt)−3 0
=
0 0

Note that these expressions can also be derived by taking the time derivative of Equations
(3.67) and (3.68).
In this example, the rotation tensor R reads:
" #
1 0
R= (3.76)
0 1

The right stretch tensor now becomes:

" #
1 + αt 0
U = RT · F = (3.77)
0 1

We can easily see that indeed C = U · U.

3.5.3 Simple shear

Figure 3.6: Deformation mode for the simple shear example

The displacement function can be written as:

x = X + αtY (3.78)
y = Y

Notice that αt is related to the rotation φ of the vertical material lines in the origin by
αt = tan φ.
The deformation gradient for this deformation pattern is:
" #
∂x 1 αt
F= = (3.79)
∂X 0 1

The right and left Cauchy–Green stretch tensors now become:

" #
1 αt
C = FT · F = (3.80)
αt α2 t2 + 1
32 Kinematics and Strain
" #
T
α2 t2 + 1 αt
B=F·F = (3.81)
αt 1
The Green–Lagrange strain tensor relates to the right Cauchy–Green stretch tensor via:
" #
1 1 0 αt
E = 2 (C − 1) = 2 (3.82)
αt α2 t2

The Euler–Almansi strain tensor relates to the left Cauchy–Green stretch tensor via:
" # " # " #
1 1 1 0 1 1 −αt 1 0 αt
e = 2 (1 − B ) = 2
−1
−2 =2 (3.83)
0 1 −αt α2 t2 + 1 αt −α2 t2

In the definition of the Green–Lagrange strain and the Euler–Almansi strain tensile or com-
pressive terms evolve in the vertical direction, although only shear deformation is imposed.
The classical linear strain tensor reads:
" #
 ← − →
−  0 αt
ε = 21 u∇0 + ∇0 u = 21 ((F − 1) + (FT − 1)) = 12 (3.84)
αt 0

The velocity field for this problem is:

vx = αY (3.85)
vy = 0

The velocity gradient reads:

" #" # " #
∂v 0 α 1 −αt 0 α
L= · F−1 = = (3.86)
∂X 0 0 0 1 0 0

which can be split into a symmetric and skew-symmetric part, respectively:

" #
0 α
D = 21 (L + LT ) = 12 (3.87)
α 0
" #
1 T 1 0 α
W= 2 (L −L )= 2 (3.88)
−α 0
The time derivatives of the Green–Lagrange strain tensor and the Euler–Almansi strain tensor
are:
" #
T 1 0 α
Ė = F · D · F = 2 (3.89)
α 2α2 t

" #
1 0 α
ė = D − e · L − LT · e = 2 (3.90)
α −2α2 t
3.5 Examples 33

The deformation gradient can be decomposed via polar decomposition. The rotation tensor
and the right stretch tensor can be determined using Equations (3.56), (3.57), (3.59) and
(3.60). " #
1 2 αt
R= √ (3.91)
α2 t2 + 4 −αt 2
" #
T 1 2 αt
U=R ·F= √ (3.92)
α2 t2 + 4 αt α2 t2 + 2

3.5.4 Rigid rotation

Figure 3.7: Deformation mode for the rigid rotation example

The displacement function can be written as:

x = X cos αt − Y sin αt (3.93)

y = X sin αt + Y cos αt

The deformation gradient for this deformation pattern is:

" #
∂x cos αt − sin αt
F= = (3.94)
∂X sin αt cos αt

The right and left Cauchy–Green stretch tensors now become:

" #
1 0
C = FT · F = (3.95)
0 1
" #
1 0
B = F · FT = (3.96)
0 1
The Green–Lagrange strain tensor relates to the right Cauchy–Green stretch tensor via:

E = 12 (C − 1) = 0 (3.97)

The Euler–Almansi strain tensor relates to the left Cauchy–Green stretch tensor via:

e = 12 (1 − B−1 ) = 0 (3.98)

The classical linear strain tensor reads:

" #
←
− →
− cos αt − 1 0
ε = 12 (u∇0 + ∇0 u) = 1
2 ((F − 1) + (FT − 1)) = (3.99)
0 cos αt − 1
34 Kinematics and Strain

Note that in case of a rigid rotation, the Green–Lagrange strain tensor and the Euler–Almansi
strain tensor both equal the correct value 0, where the classical strain tensor shows non-zero
components. The error made becomes larger with increasing rotation, until αt = π, half a
revolution. The velocity field for this problem is:

vx = −αX sin αt − αY cos αt (3.100)

vy = αX cos αt − αY sin αt

The velocity gradient reads:

" #
∂v 0 −α
L = · F−1 = (3.101)
∂X α 0

which can be split into a symmetric and skew-symmetric part, respectively:

" #
0 0
D = 21 (L + L ) =
T
(3.102)
0 0
" #
1 T
0 −α
W= 2 (L −L )= (3.103)
α 0
The time derivatives of the Green–Lagrange strain tensor and the Euler–Almansi strain tensor
are:

Ė = FT · D · F = 0 (3.104)

ė = D − e · L − LT · e = 0 (3.105)

The deformation gradient can be decomposed via polar decomposition. The rotation tensor
and the right stretch tensor can be determined using Equations (3.56), (3.57), (3.59) and
(3.60). " #
cos αt − sin αt
R= (3.106)
sin αt cos αt
" #
1 0
U= (3.107)
0 1
which is trivial for a deformation gradient that represents rotation without deformation.
Chapter 4

Stress and Equilibrium

4.1 The true stress

unit normal vector n

Figure 4.1: Infinitesimal force

To define the true stress or Cauchy stress, we consider an infinitesimal surface element dS
such as given in Figure 4.1, with a unit normal vector n. An infinitesimal force dF acts on
dS:
dF = T dS = Ti ei dS (4.1)

hence, T has the dimension of a stress. The cross section areas of the surfaces perpendicular
to the xi -axes are called dSi (see Figure 4.2). It holds that:

dSi = ni dS (4.2)

proof: (see Figure 4.3):


dS = 21 ab 



dS3 = 12 ac


) =⇒ dS3 = 21 abn3 = n3 dS (4.3)
c = b cos(α) 
=⇒ c = bn3 


n3 = knk cos α = cos α



35
36 Stress and Equilibrium

Figure 4.2: Infinitesimal surfaces

analogously:
dS1 = n1 dS
dS2 = n2 dS
hence:
dSi = ni dS q.e.d.
Stress components σij act on the areas dSj . The first index refers to the surface element on

Figure 4.3: Normal on a surface element

which the force acts and the second index refers to the direction of the force (see Figure 4.4).
The equilibrium condition in the x1 -direction yields:
T1 dS = σ11 dS1 + σ21 dS2 + σ31 dS3 = (σ11 n1 + σ21 n2 + σ31 n3 ) dS
hence
T1 = σj1 nj
The same applies to the other directions:
T2 = σj2 nj
T3 = σj3 nj
hence:
Ti = σji nj (4.4)
4.2 Equilibrium equations 37

Figure 4.4: Stresses on a surface element

4.2 Equilibrium equations

The stress state must also fulfill the internal equilibrium conditions. These internal equilib-
rium conditions were treated in preceding mechanics courses. A short summary of force and
moment balance equations leading to the equilibrium equation and the symmetry of the stress
tensor are given below.

4.2.1 Force balance

The starting point of this section is a volume part V . At the surface of volume V , the forces
Ti per unit surface act.

Figure 4.5: An arbitrary portion of a body in equilibrium.

The total outer force is then:

Z Z
Ki = Ti dS = σji nj dS (4.5)
S S

with S the surface of the current volume part V . Using the theorem of Gauss, the next
equation is found:
38 Stress and Equilibrium

Z
Ki = σji,j dV (4.6)
V

The force Ki must be in equilibrium with the possible volume forces (e.g. gravity force, inertia
force, induction force). The volume forces are defined as Fi per unit of volume. Now we can
write:
Z Z
σji,j dV + Fi dV = 0 (4.7)
V V

The above equation applies for every volume part, so:

σji,j + Fi = 0 (4.8)

4.2.2 Moment balance

The starting point is the same body that is in equilibrium depicted in Figure 4.5.
The torque τ that is created by a force vector T at coordinates x around the origin is found
using the cross-product of these two vectors which gives three equations representing the
moment around each axis:

τ = x × T = ǫijk xj Tk (4.9)
is not equal, misses a vector
where ǫ is the permutation operator.
The total torque due to all the forces acting on the surface is represented by a surface integral
and those due to body forces are given in a volume integral. The sum of all must vanish in
equilibrium: Z Z
ǫijk xj Tk dS + ǫijk xj Fk dV = 0 (4.10)
S V

Using Equation (4.4) and the divergence theorem:

∂
Z Z Z
ǫijk xj Tk dS = ǫijk xj σmk nm dS = (ǫijk xj σmk ) dV (4.11)
∂xm
S S V

Using the chain rule to expand the differentiation and combining the two volume integrals in
one gives: Z
ǫijk δjm σmk + ǫijk xj (σmk,m + Fk ) dV = 0 (4.12)
V

The term in the parentheses vanishes when the forces in the body are in equilibrium due to
Equation (4.8). The remaining term in the integral must vanish for any arbitrary volume
which gives ǫijk σjk = 0. This results in the following three identities

σ12 = σ21 , σ23 = σ32 , σ13 = σ31 (4.13)

proving the symmetry of the Cauchy stress tensor, i.e. σij = σji
4.3 Alternative stress representations 39

4.3 Alternative stress representations

In large deformation problems, it is sometimes convenient to use other stress measures than
the Cauchy stress. This is especially true for numerical calculations, since the Cauchy stress is
defined in the ‘current’ state, which is not known beforehand in a calculation step. However,
stress measures can be transformed into each other, so the above argument does not make an
essential difference. The subject is related to the different definitions of strain in Chapter 3.

4.3.1 Stress energy rate

It can easily be derived that the internal mechanical energy rate (also called stress energy
rate or stress power) for a body is given by
Z Z
Ẇ = ẇ dV = σ : D dV (4.14)
V V

where the specific internal mechanical energy rate ẇ (per unit volume) is defined in the current
configuration as:
ẇ = σ : D (4.15)
Alternative stress measures are defined such that similar relations hold for different strain
measures. A pair of stress and strain measures for which the product results in the stress
energy rate are called work conjugate pairs.

4.3.2 The Kirchhoff stress tensor

Since the volume in the current volume changes upon deformation, it is more convenient to
express the stress power with respect to the undeformed volume dV0

ẇ dV = ẇ0 dV0 = σ : DJ dV0 (4.16)

For this reason, the Kirchhoff stress tensor τ is introduced, which is also defined on the
current configuration. This stress tensor is simply defined by

τ =Jσ (4.17)

The stress power with respect to the undeformed volume can be written as

ẇ0 = τ : D (4.18)

The Cauchy stress tensor and the Kirchhoff stress tensor are both frequently used in consti-
tutive relations.

4.3.3 The first Piola–Kirchhoff stress tensor

Another measure of stress can be obtained by realizing that the Kirchhoff stress is symmetric
and therefore, the multiplication with D can as well be performed with L and using (3.47):
 
ẇ0 = τ : L = τ : Ḟ · F−1 (4.19)
40 Stress and Equilibrium

By writing this equation in index notation, it can be proven that this can also be written as

ẇ0 = (τ · F−T ) : Ḟ (4.20)

The part between brackets is apparently work conjugated to the deformation gradient and is
called the first Piola–Kirchhoff stress P:

P = Jσ · F−T (4.21)

or

P = τ · F−T (4.22)
T
τ =P·F (4.23)

Clearly, P is non-symmetric. It is sometimes called the non-symmetric nominal stress tensor.

If a surface does not change orientation, the first Piola–Kirchhoff stress gets the meaning of
engineering stress, the force per unit undeformed area.
The stress power with respect to the undeformed configuration can now be written as

ẇ0 = P : Ḟ (4.24)

4.3.4 The second Piola–Kirchhoff stress tensor

The symmetric nominal stress tensor, called the second Piola–Kirchhoff stress tensor S, is
defined by
S = F−1 · P = F−1 · τ · F−T = JF−1 · σ · F−T (4.25)
Therefore
P=F·S (4.26)
and
σ = J −1 F · S · FT (4.27)
The second Piola–Kirchhoff stress tensor is a symmetric tensor which is fully defined on the
reference configuration.
The stress power with respect to the undeformed configuration can be written as

ẇ0 = P : Ḟ = (F · S) : Ḟ (4.28)

And with 2E = C − 1 = FT · F − 1:

Ċ = 2Ė = ḞT · F + FT · Ḟ (4.29)

finally, using the symmetry of the second Piola–Kirchhoff stress that easily follows from (4.25)
1
ẇ0 = S : Ė = S : Ċ (4.30)
2
This shows that the second Piola–Kirchhoff stress is work conjugated to the Green–Lagrange
strain.
Chapter 5

Hyperelastic material models

5.1 Introduction

Rubber and rubber-like materials remain elastic at deformations that easily exceed 100%
elongation. In Figure 5.1 e.g., the results are plotted of a tensile test on natural rubber in
which the a stretch ratio of 15 is reached. It is therefore essential that models for the me-
chanical behaviour of these elastomeric materials are valid for large deformations. The elastic
behaviour of elastomers is commonly described by so-called hyperelastic material models. In
these models a strain energy density function ψ is defined, from which the stresses can be
derived as a function of the deformation. The strain energy density represents the energy
that is stored in the material per unit undeformed volume. For elastic materials, the rate of
change of ψ is equal to the stress power ẇ0 in (4.18). The strain energy will be released if the
load is released.

16

14

12
engineering stress (Mpa)

10

0
0 2 4 6 8 10 12 14 16
engineering strain (-)

Figure 5.1: Stress–strain curve for a uniaxial tensile test on natural rubber.

41
42 Hyperelastic material models

5.2 General formulation

As a first approach, it can be considered that the (local) deformation is completely determined
by the deformation gradient F and therefore ψ must be a function of the deformation gradient
ψ(F). A rigid translation should not influence the stress.
Although a rigid rotation will influence the stress, it does not alter the deformation of the
material and should therefore not influence the strain energy. Since the deformation gradient
can be split in a rotational part and a deformation part by using the polar decomposition

F=R·U (5.1)

it can be concluded that ψ is only a function of the right stretch tensor U, which is symmetric.
In most models the square of U is used, referred to as the right Cauchy–Green tensor:

C = U · U = FT · F (right Cauchy–Green tensor) (5.2)

As shown above, the right Cauchy–Green tensor can easily be obtained without need for the
polar decomposition. Depending on the measure of deformation, one can write for the first
Piola–Kirchhoff stress P and the second Piola–Kirchhoff stress S with Equations (4.24) and
(4.30):

∂ψ
P= (5.3)
∂F
∂ψ ∂ψ
S=2 = (5.4)
∂C ∂E
The Cauchy stress can be calculated with Equation (4.27).

5.2.1 Isotropic material

A simplification is achieved by considering that the material is initially isotropic. This means
that the initial orientation of a specimen is not relevant for the stress–strain response. It does
not mean that a change in strain will give a stress change independent of the pre-strain. In
general the stress rate will depend on the pre-strain, but this is not excluded in the adopted
definition of isotropy.
If the material is isotropic, the strain energy density function can only be a function of the
invariants of the deformation measure. The invariants have the property that they do not
change value if the coordinate system is rotated. Therefore, many models use the three
invariants of the right Cauchy–Green stretch tensor as the deformation variables in a strain
energy density function.
ψ = ψ(I1 , I2 , I3 ) (5.5)
with the invariants:

I1 = tr(C) = C11 + C22 + C33 (5.6)

1 2
tr (C) − tr(C2 ) = C11 C22 + C11 C33 + C22 C33 − C122 2 2

I2 = − C23 − C13 (5.7)
2
2 2 2
I3 = det C = C11 C22 C33 + 2C12 C23 C13 − C11 C23 − C22 C13 − C33 C12 (5.8)
5.2 General formulation 43

It should be noted that the 3 invariants are sometimes defined differently. In fact any com-
binations of the above defined invariants is by itself an invariant of C. The invariants are
strongly related to the 3 eigenvalues.
Because of the straightforward physical interpretation also the square root of I3 is often used
instead of I3 . This equals the determinant of F and represents the ratio of volumes after and
before deformation:
p V
J = I3 = det F = (5.9)
V0
We then write
ψ = ψ(I1 , I2 , J) (5.10)

5.2.2 Stress calculations based on invariants

The 2nd Piola–Kirchhoff stress is determined from the strain energy density funtion by Equa-
tion (5.4):
∂ψ ∂I1 ∂ψ ∂I2 ∂ψ ∂J
S=2 +2 +2 (5.11)
∂I1 ∂C ∂I2 ∂C ∂J ∂C
The derivatives of I1 , I2 and J with respect to C are derived from Equations (2.89)–(2.91)
and (5.9) where the symmetry of C is used.

∂I1 ∂(1 : C)
= =1 (5.12)
∂C ∂C

∂ 21 tr2 (C) − tr(C · C)

 
∂I2
= = I1 1 − C (5.13)
∂C ∂C
∂I3 ∂ det C
= = det(C)C−1 = I3 C−1 (5.14)
∂C ∂C
∂J 1 ∂ det C 1
= (det C)−1/2 = JC−1 (5.15)
∂C 2 ∂C 2

so that
∂ψ ∂ψ ∂ψ
S=2 1+2 (I1 1 − C) + JC−1 (5.16)
∂I1 ∂I2 ∂J

With σ = J1 F · S · FT and F · FT = B we obtain for the Cauchy stress1 :

2 ∂ψ 2 ∂ψ
∂ψ
σ= B+ I1 B − B2 + 1 (5.17)
J ∂I1 J ∂I2 ∂J

The term B2 can be expressed in B and B−1 by using Equation (2.85):

2 ∂ψ 2 ∂ψ
∂ψ
σ= B+ I2 1 − J 2 B−1 + 1 (5.18)
J ∂I1 J ∂I2 ∂J
1
The invariants of B are equal to the invariants of C: I1 (C) = I1 (B), I2 (C) = I2 (B) and I3 (C) = I3 (B).
 44 Hyperelastic material models

Compressible Neo-Hookean material model (Do we need to know this for the exam?)

One possible formulation of the strain energy density function is

G λ
(compressible neo-Hookean material model) ψ= (I1 − 3) − G ln J + (ln J)2 (5.19)
2 2
with material constants G and λ. This model only depends on the first and third invariant of
C and is known as the compressible neo-Hookean material model. In the undeformed shape
C = 1, hence I1 = 3 and J − 1 and therefore ψ = 0. The second Piola–Kirchhoff stress can
now be obtained from (5.16):

S = G(1 − C−1 ) + λ(ln J)C−1 (5.20)

The Cauchy stress is expressed as:

G λ
σ= (B − 1) + (ln J)1 (5.21)
J J
As expected, for an undeformed neo-Hookean material, the stress is zero.
The derivation of the material stiffness relation is out of the scope of this course, but a
quick-and-dirty evaluation at the undeformed state shows the asymptotic behavior to the
well-known small strain linear elastic stiffness. Taking the time derivative of (5.20) yields

∂C−1 λJ˙ −1 ∂C−1

Ṡ = −G + C + λ ln J (5.22)
∂t J ∂t

By using (2.88), (5.15) and J˙ = ∂J : Ċ

∂C
1
Ṡ = GC−1 · Ċ · C−1 + λ (C−1 : Ċ)C−1 − λ ln JC−1 · Ċ · C−1 (5.23)
2

From (3.24) and (3.39) follows that Ċ = 2Ė = 2FT · D · F. Substituting this in the previous
equation gives

Ṡ = 2GC−1 · FT · D · F · C−1 + λ [C−1 : (FT · D · F)] C−1 + λ ln JC−1 · FT · D · F · C−1 (5.24)

For small deformations and rotations σ ≈ S, F ≈ 1, C ≈ 1 and J ≈ 1 and we obtain

σ̇ ≈ 2GD + λ1(1 : D) = 2GD + λ tr(D)1 (5.25)

The meaning of the material parameters G and λ is hereby demonstrated to be equal to the
E Eν
Lamé parameters for isotropic linear elastic material, with G = 2(1+ν) and λ = (1+ν)(1−2ν) .

5.3 Nearly incompressible material

Rubber is often considered incompressible, which is a reasonable assumption if the bulk mod-
ulus is compared to the shear modulus. Actually, the bulk modulus of rubber is lower than
that of steel, but the ratio with the shear modulus is much higher for rubber (see Table 5.1).
5.3 Nearly incompressible material 45

Table 5.1: Shear and bulk moduli

shear modulus (MPa) bulk modulus (MPa)

steel 75000–80000 160000
silicone rubber 0.3–20 1500–2000

For incompressible materials, the volume does not change and therefore the third invariant
of the Cauchy–Green stretch tensor I3 = J = 1. The strain energy density function is in that
case only a function of the first and second invariant:
ψ = ψ̂(I1 , I2 ) (5.26)
A problem, however, is that the stress can only be obtained up to an unknown hydrostatic
component. Or—in other words—if the deformation does not change upon a change in hy-
drostatic stress, the stress cannot be determined from the deformation only.
Although there are mathematical ways to calculate the hydrostatic stress component from
another set of equations and mathematically impose complete incompressibility, the common
solution is to allow some compressibility in the model. In fact this is physically correct.
Formally, this can be described by a general strain energy density function as in (5.10), but in
so-called nearly incompressible material models the part in the strain energy density function
that depends on the Jacobian is split-off in an additive way from the rest:
ψ = ψ̂(I1 , I2 ) + U (J) (5.27)
The second Piola–Kirchhoff stress can now be calculated from
∂ψ ∂ ψ̂ ∂U ∂J
S=2 =2 +2 volumetric strain (5.28)
∂C ∂C ∂J ∂C
A (minor) drawback of this formulation is that for a purely dilatational deformation I1 and I2
are still influenced, as can be seen in (5.6) and (5.7). Therefore, the hydrostatic stress is not
only a function of U , but also of ψ̂. It is therefore common to define a strain energy density
function that only depends on the distortional part of the right Cauchy–Green tensor:
−1/3
C̄ = I3 C = J −2/3 C (5.29)
Because an algebraic combination of invariants remains invariant, we can define distortional
invariants as
I¯1 = I1 (C̄) = J −2/3 I1 (5.30)
I¯2 = I2 (C̄) = J −4/3 I2 (5.31)
and define a—different—strain energy density function for nearly incompressible materials as
ψ = ψ̄(I¯1 , I¯2 ) + U (J) (5.32)
The model is now split in a purely distortional (isochoric) part and a purely dilatational
(volumetric) part and the second Piola–Kirchhoff stress is calculated from
∂ ψ̄ ∂U ∂J
S=2 +2 (5.33)
∂C ∂J ∂C
46 Hyperelastic material models

5.3.1 Stress calculations based on distortional invariants

The 2nd Piola–Kirchhoff stress is determined from the strain energy density funtion by Equa-
tion (5.33) where ψ̄ is a function of the distortional invariants I¯1 and I¯2 .

∂ ψ̄ ∂ I¯1 ∂ ψ̄ ∂ I¯2 ∂U ∂J
S=2 ¯ +2 ¯ +2 (5.34)
∂ I1 ∂C ∂ I2 ∂C ∂J ∂C
The derivatives with respect to C follow from:

∂ I¯1 ∂ I¯1 ∂I1 ∂ I¯1 ∂J

= + (5.35)
∂C ∂I1 ∂C ∂J ∂C

∂ I¯2 ∂ I¯2 ∂I2 ∂ I¯2 ∂J

= + (5.36)
∂C ∂I2 ∂C ∂J ∂C
From (5.30) and (5.31) follows

∂ I¯1
= J −2/3 (5.37)
∂I1
∂ I¯2
= J −4/3 (5.38)
∂I2
∂ I¯1 2
= − J −5/3 I1 (5.39)
∂J 3
∂ I¯2 4
= − J −7/3 I2 (5.40)
∂J 3

and with Equations (5.12), (5.13) and (5.15) derivatives of the distortional invariants become:

∂ I¯1
 
1
= J −2/3 1 − I1 C−1 (5.41)
∂C 3

∂ I¯2
 
2
= J −4/3 I1 1 − C − I2 C−1 (5.42)
∂C 3

and
   
∂ ψ̄ −2/3 1 ∂ ψ̄ −4/3 2 ∂U
S=2 ¯ J 1 − I1 C −1
+2 ¯ J I1 1 − C − I2 C −1
+ JC−1 (5.43)
∂ I1 3 ∂ I2 3 ∂J

With σ = J1 F · S · FT and F · FT = B we obtain for the Cauchy stress:

   
∂ ψ̄ −5/3 1 ∂ ψ̄ −7/3 2 2 ∂U
σ=2 ¯J B − I1 1 + 2 ¯ J I1 B − B − I2 1 + 1 (5.44)
∂ I1 3 ∂ I2 3 ∂J

With help of the Cayley–Hamilton equation (2.84) this equation can also be written as:
   
∂ ψ̄ −5/3 1 ∂ ψ̄ −7/3 1 2 −1 ∂U
σ=2 ¯J B − I1 1 + 2 ¯ J I2 1 − J B + 1 (5.45)
∂ I1 3 ∂ I2 3 ∂J
5.4 Specific material models 47

5.4 Specific material models

5.4.1 Volumetric part of the strain energy density
The volumetric part of the strain energy density is often chosen to be
1
U = K(J − 1)2 (5.46)
2
To determine the effect on the stress the following derivative is determined:

∂U
= K(J − 1) (5.47)
∂J
An alternative formulation for U that was also used in (5.19) is:

1
U= K(ln J)2 (5.48)
2
For this function the derivative is:
∂U K
= ln J (5.49)
∂J J

5.4.2 Neo-Hookean material model

The distortional part of the strain energy density function for the Neo-Hookean model is

ψ̄ = C1 (I¯1 − 3) (5.50)

The Neo-Hookean model can be obtained from statistical mechanics when vulcanized rubber
is regarded as a random network of long molecular chains, pinned at a few positions by
cross-links.
The following derivatives are required in the calculation of the stress:

∂ ψ̄
= C1 (5.51)
∂ I¯1
∂ ψ̄
=0 (5.52)
∂ I¯2

5.4.3 Mooney–Rivlin material model

The isochoric part of the strain energy density function for the Mooney–Rivlin model is

ψ̄ = C1 (I¯1 − 3) + C2 (I¯2 − 3) (5.53)

With C2 = 0 this model degenerates to the Neo-Hookean model.

The following derivatives are required in the calculation of the stress:

∂ ψ̄
= C1 (5.54)
∂ I¯1
∂ ψ̄
= C2 (5.55)
∂ I¯2
48 Hyperelastic material models

5.4.4 Yeoh material model

The isochoric part of the strain energy density function for the nearly incompressible Yeoh
material model is

2
3
ψ̄ = C1 I¯1 − 3 + C2 I¯1 − 3 + C3 I¯1 − 3


(5.56)

In order to calculate the constitutive relations the following derivatives are determined:

∂ ψ̄
2
= C1 + 2C2 I¯1 − 3 + 3C3 I¯1 − 3


¯ (5.57)
∂ I1
∂ ψ̄
=0 (5.58)
∂ I¯2

5.4.5 Generalized Rivlin model

All above models are polynomials in the distortional invariants. A generalized form is at-
tributed to Rivlin as
n X
X n
ψ̄ = Cij (I¯1 − 3)i (I¯2 − 3)j with C00 = 0 (5.59)
i=0 j=0

with derivatives
n X
n
∂ ψ̄ X
¯1 − 3 i−1 I¯2 − 3 j



= = iC ij I (5.60)
∂ I¯1 i=1 j=0
n X
n
∂ ψ̄
i
j−1
= jCij I¯1 − 3 I¯2 − 3
X
= (5.61)
∂ I¯2 i=0 j=1

For the determination of the material stiffness tensors, the reader is referred to e.g. [1, 3].

5.4.6 Examples
For an almost incompressible material, we can easily estimate the stress–strain response of a
uniaxial tensile test by considering the stretch ratios λ1 , λ2 and λ3 . For an incompressible
material J = √λ1 λ2 λ3 = 1 and for a uniaxial tensile test in 1-direction of an isotropic material:
λ2 = λ3 = 1/ λ1 . The deformation gradient, right and left Cauchy–Green tensor and their
inverses then become:  
λ1 0 0
 √ 
(Deformation gradient)
F=  0 1/ λ1 0  (5.62)
√

0 0 1/ λ1

and    
λ21 0 0 λ−2
1 0 0
   
C=B=
0 1/λ1 0  C−1 = B−1 = 
 0 λ1 0
 (5.63)
0 0 1/λ1 0 0 λ1

Right Cauchy-Green tensor

5.4 Specific material models 49

The invariants can be calculated from (5.6) and (5.7):

2
I1 = λ21 + λ22 + λ23 = λ21 + (5.64)
λ1
1
I2 = λ21 λ22 + λ21 λ23 + λ22 λ23 = 2λ1 + (5.65)
λ21
The Cauchy stress follows from (5.45) and can be written in components:
   
4 ∂ ψ̄ 2 1 4 ∂ ψ̄ 1 ∂U
σ1 = ¯ λ1 − + ¯ λ1 − 2 + (5.66)
3 ∂ I1 λ 3 ∂ I2 λ1 ∂J
   
2 ∂ ψ̄ 1 2 ∂ ψ̄ 1 ∂U
σ2 = − ¯ λ21 − − ¯ λ1 − 2 + (5.67)
3 ∂ I1 λ 3 ∂ I2 λ1 ∂J

The value of ∂U
∂J can be derived from the boundary condition that σ2 = 0 independent of the
shape of U . We make a small error because for a particular function U , some compressibility
is allowed and λ1 λ2 λ3 is not exactly equal to one.
Substituting the obtained value for ∂U in σ1 gives:
∂J
   
∂ ψ̄ 2 1 ∂ ψ̄ 1
σ 1 = 2 ¯ λ1 − + 2 ¯ λ1 − 2 (5.68)
∂ I1 λ ∂ I2 λ1
The engineering stress can be derived by multiplication with the ratio of current and initial
cross section area λ2 λ3 which is here equal to 1/λ1 . Note that this is also the value of the
first Piola–Kirchhoff stress P1 .
   
∂ ψ̄ 1 ∂ ψ̄ 1
P1 = 2 ¯ λ1 − 2 + 2 ¯ 1 − 3 (5.69)
∂ I1 λ ∂ I2 λ1
This equation can directly be compared with the engineering stress–strain curve, derived from
experiments. The derivative with respect to λ1 is
∂ 2 ψ̄ ∂ 2 ψ̄
       
dP1 ∂ ψ̄ 2 2 1 ∂ ψ̄ 3 1
=2 ¯ 1+ 3 +2 λ1 − +2 ¯ +2 λ1 − 2 (5.70)
dλ1 ∂ I1 λ ∂λ1 ∂ I¯1 λ ∂ I2 λ41 ∂λ1 ∂ I¯2 λ
For the initial material stiffness (comparable to the Young’s modulus) at λ1 = 1 this yields
 
dP1 ∂ ψ̄ ∂ ψ̄
(λ1 = 1) = 6 + (5.71)
dλ1 ∂ I¯1 ∂ I¯2
The relation between shear modulus and Young’s modulus is 2G(1 + ν) = E, which for an
incompressible material (ν = 0.5) gives 3G = E and therefore
 
∂ ψ̄ ∂ ψ̄
G=2 + (5.72)
∂ I¯1 ∂ I¯2
In the generalized Rivlin model, this means that G = 2(C10 + C01 ) in order to model the
initial stiffness correct.
We can now fit the material parameters for the Neo-Hookean and Mooney–Rivlin models to
the experimental stress–strain values as presented in Figure 5.1, where we try to match the
initial shape of the curve as good as possible. The results are presented in Figure 5.2. The
parameter for the Neo-Hookean model is C1 = 0.2 and the parameters for the Mooney–Rivlin
model are C1 = 0.03 and C2 = 0.017. Only strains up to 100% are plotted. It can be
observed that both models represent the same initial stiffness and that the Mooney–Rivlin
model represents the curve well up to strains of about 50%.
50 Hyperelastic material models

0.7
experiment
Neo-Hookean
0.6 Mooney-Rivlin
engineering stress (Mpa)

0.5

0.4

0.3

0.2

0.1

0
0 0.2 0.4 0.6 0.8 1
engineering strain (-)

Figure 5.2: Stress–strain curve for a uniaxial tensile test on natural rubber.
Chapter 6

Plasticity theory

6.1 Introduction
In metals the observed inelastic deformation above a critical stress is referred to as plasticity.
Most metals have an atomic structure which consists of crystals (grains) oriented arbitrarily
(isotropic) or in a preferred way (textured) which affects the observed plastic behavior in
macroscale. Plasticity occurs by dislocation motion within this crystallographic arrangement
therefore it is expected to be isotropic or anisotropic with regard to how the grains are
oriented. In order for dislocations to move there has to be stress resolved in the crystal
structure meaning that the total strain has to be made up of an elastic part and a plastic part.
It is therefore necessary to split the total strain tensor accordingly. This split however cannot
be defined without choosing a certain direction for the strain tensors. As the dislocations
become mobilized they lead to creation of more dislocations which in turn start to hinder the
mobility of themselves. Therefore larger stresses are required as the amount of accumulated
plastic strain increases resulting in the strain hardening phenomenon.
In order to model the phenomena described very briefly above this chapter starts with a
one-dimensional formulation of the elastic-plastic split and a yield criterion. Afterwards
a general three dimensional formulation is developed by making use of Drucker’s postulate.
Commonly used yield surfaces are introduced. In order to deal with strain hardening different
formulations are given and the solution procedure is introduced.

6.2 One-dimensional plasticity

The uniaxial tensile test is an example of one-dimensional plasticity due to the fact that
in that configuration there is only one component of the stress tensor. Although lateral
components of the strain tensor exist since the yield criterion is defined in the stress space
these components are irrelevant in the formulations.
The motivation for plasticity in this case is based on a frictional device sketched in figure
6.1. On this device when an external force σ starts to act first an elongation in the spring
according to the spring constant E will be observed. When the external force reaches the
frictional force σy sliding will occur. The elongation in the spring is reversible however the
deformation due to slip is not.
The resulting force displacement curve is sketched in figure 6.2. It is also shown that when
the force is decreased after a certain amount of sliding the deformation will be reversed due

51
52 Plasticity theory

E
σ
σy

Figure 6.1: One-dimensional frictional device demonstrating elastic-rigid plastic behavior.

to unloading of the spring with a proportion of E.

σ
σy

−σy

Figure 6.2: Stress-strain response corresponding to the frictional device.

6.2.1 Governing equations

We can formulate the observations on the frictional device into equations that will be used
to solve the system. The aim is now to solve for the resulting stress (or strain) for a given
strain (or stress) increment and also the relation between the stress and strain increment.
First of all it is clear that the total strain is composed of two parts, namely the elastic and
plastic parts:
ε = εe + εp (6.1)
This relation is also used for calculating the stress in the material. Since in a strain driven
formulation the total strain is known, the stress can be related to strain as follows:

σ = E εe = E (ε − εp ) (6.2)

This means that in order to find the stress, the plastic strain has to be found first.
It is clear considering the device that the stress level cannot exceed the frictional stress in
either direction. To incorporate the bidirectional response we can define a yield criterion
which describes the admissible set of stresses:

φ(σ) = |σ| − σy ≤ 0 (yield criterion) (6.3)

Using (6.3) we can formulate conditions for plastic strain:

6.3 Plastic deformation at multi-axial stress states 53

• ε̇p = 0 if |σ| < σy

• ε̇p > 0 if σ = σy > 0

• ε̇p < 0 if σ = σy < 0

where the rate of plastic strain is defined as ε̇p = ∂εp /∂t and describes the change in the
plastic strain instead of its total value. These conditions can be described as follows: If the
absolute value of the stress is less than the yield stress there is no change in the plastic strain;
If the stress is equal to the yield stress on the positive side, there will be a positive change in
the plastic strain; If the stress is equal to the yield stress on the negative side there will be a
negative change in the plastic strain.
Note that conditions 2 and 3 can actually be represented by a single equation:

ε̇p = λ̇ sign(σ) (6.4)

signum function
where λ̇ is introduced as the change of the absolute accumulated value of the plastic slip.
Note that according to this definition although ε̇p can be negative λ̇ is always positive.
Combining the above we obtain the following set of conditions that will be used to find the
value of the plastic strain rate:

λ̇ ≥ 0
φ(σ) ≤ 0 (Kuhn-Tucker conditions)

λ̇ φ(σ) = 0 (6.5)

These conditions in plasticity literature are usually referred to as Kuhn-Tucker conditions.

One additional condition that we can obtain by differentiating the last Kuhn-Tucker condition
when φ(σ) = 0 is:
λ̇ φ̇(σ) = 0 (consistency condition) (6.6)
which is called the consistency condition.
Using the consistency condition (6.6) we can solve for the unknown plastic strain rate which
is required to update the value of the stress. From (6.3) we obtain:
∂φ
φ̇(σ) = σ̇ = sign(σ) E (ε̇ − ε̇p )
∂σ
= sign(σ) E ε̇ − E λ̇ = 0
⇒ λ̇ = ε̇ sign(σ) (6.7)

Equation (6.7) reveals that once the yield stress is reached all strain increment becomes plastic
which is the solution that we are after.

6.3 Plastic deformation at multi-axial stress states

In the preceding section plasticity is formulated in a one-dimensional case in which certain
simplifications occur naturally such as that a straightforward yield criterion can be used and
the direction of plastic strain rate is equal to the sign of stress.
When multi-axial stress states are concerned though important considerations must be made
that satisfy physical as well as mathematical requirements. A three dimensional yield criterion
54 Plasticity theory

must be defined to account for arbitrary stress states and the direction of the plastic strain rate
must be postulated in order to find a solution. Other considerations such as the Kuhn-Tucker
and the consistency conditions can be generalized to multi-axial form in a straightforward
manner by replacing the scalar components with tensorial counterparts.

6.3.1 The yield surface

For the uniaxial loading case, the transition from elastic to plastic deformation is defined by
a yield criterion which is a function of a scalar component of stress. For multi-axial stress
states this transition is characterized by a yield surface where more components of stress play
a role. In order to illustrate the definition of a yield surface consider first a plane stress state.
The plane stress state can be characterized by two principal stresses. Plotting these two
principal stresses along the x- and y-axis and assuming an arbitrary yield surface definition
gives figure 6.3. The bounding line represents the yield surface for which it holds that when
the stress state is in the interior of the yield surface the deformation stays elastic and plastic
deformation takes place when the stress state is on the yield surface (by definition the stress
state can never be outside the yield surface). The shape and size of the yield surface can
change during plastic deformation due to different types of hardening. Note however that the
intersection point of the yield surface and one of the principal axis represents the uniaxial
yield stress.

When stress state is in the interior of the

yield surface -> deformation is elastic
When the stress state is ON the yield surface
-> Plastic deformation

Figure 6.3: An arbitrary yield surface

6.3.2 Normality rule and the convexity of the yield surface

The arbitrary yield surface that is defined in the previous section has to have some proper-
ties in order to represent important observations in metal plasticity. Additionally from the
formulations that will follow we can obtain a definition for the direction of the plastic strain
tensor.
The uniaxial tensile curve of a metal possesses some important characteristics that is not
yet considered in the above treatment regarding the stability of deformation. In figure 6.4
two tensile curves are plotted which represent instable mechanical behavior. The reasoning
is that in thees cases a negative strain rate accompanies a positive stress rate or vice-versa.
Therefore for a stable stress-cycle the work done by the external agency (loads) is positive:
6.3 Plastic deformation at multi-axial stress states 55

s a s b

(instable mechanical behaviour)

Softning

e e
Figure 6.4: Two instable stress-strain curves

I
σ : dε ≥ 0 (6.8)

<-- Stable cycle

Figure 6.5: The stress cycle

A continuation of this reasoning is given by Drucker in which a more restrictive definition of

stable and instable deformation is given:
• σ̇ : ε̇ > 0 ⇒ stable (hardening),

• σ̇ : ε̇ = 0 ⇒ neutral (rigid plastic),

• σ̇ : ε̇ < 0 ⇒ instable (softening).

Since the inequalities will not change when multiplied with a positive time increment it can
be further elaborated that for a stable deformation we need:

dσ : dε ≥ 0 (6.9)

but since dσ : dεe is always positive Drucker has postulated that:

dσ : dεp ≥ 0 (6.10)

In words what (6.10) says is that the plastic strain increment cannot oppose the stress incre-
ment. It should be noted that this postulate applies generally for metals but other materials
for instance soil do not have to possess this property.
The consequences of Drucker’s postulate are very significant in plasticity theory. Consider
for instance the arbitrary yield surface once again with a closed stress path starting from σ a
56 Plasticity theory

onto the yield surface σ y to an increment of ∆σ on the yield surface (consider hardening)
and then back to σ a as shown in figure 6.6 which would result in a one-dimensional results
as shown in figure 6.5.

Figure 6.6: Stress path on an arbitrary yield surface

The mechanical work done by an external agency in this cycle according to Figure 6.6 is
calculated as:
I I
(σ − σ a ) : dε = (σ − σ a ) : d(εe + εp ) ≥ 0 (6.11)

Carrying out the elastic-plastic gives:

I I
(σ − σ a ) : dε + (σ − σ a ) : dεp ≥ 0
e
(6.12)

The elastic part of the integral can be shown to be zero at a closed stress cycle. Since σ : εe
gives the elastic strain energy it vanishes in a cycle. Furthermore σ a is constant therefore
what remains in the integral is only dεe which also vanishes in a closed stress cycle. WHat
remains is the plastic part which only exists in the path traveling with the yield surface by
∆σ. This part can be linearized for small increments and using trapezoidal rule the integral
can be written as:
1
2 [(σ y − σ a ) + (σ y + ∆σ − σ a )] : dεp = (σ y − σ a ) : dεp + ∆σ : dεp (6.13)

Note that for a hardening material the last term is always positive (equation (6.10)). The
remaining term is the important result that we are after:

(σ y − σ a ) : dεp ≥ 0 (6.14)
The consequences of (6.14) are twofold. First, let’s consider the direct result the double
contraction operation. The difference of the stress at the yield surface and an arbitrary stress
inside is projected on the plastic strain increment resulting in a positive value. This can only
happen if the projection is done with a sharp angle (−π/2 ≤ α ≤ π/2). This is explained
further with the help of Figure 6.7, in which the vectors σ a and σ ∗a are arbitrary vectors.
On the figure the dashed areas are the admissible directions for the plastic strain increment
due to Drucker’s postulate given in equation (6.14). But at the current state of deformation
which happens on the yield surface the previous stress (σa ) is irrelevant, the postulate should
6.3 Plastic deformation at multi-axial stress states 57

Figure 6.7: Planes of the deformation increment

hold for any admissible stress. The only direction of plastic strain that allows this is the
normal to the yield surface as shown in figure 6.8. This result is known as the principle of
normality.

Figure 6.8: the deformation increment is perpendicular to the yield surface

In case dεp is not perpendicular to the yield surface, an arbitrary σa would exist for which
the postulate is violated:

(σ y − σ a ) : dεp < 0 (6.15)

which would mean that plastic deformation generates mechanical energy, see figure 6.9.
The normal to the yield surface can be found considering that the yield function is a potential
and therefore the outward normal vector is proportional to the gradient of the function. In a
∂φ
space spanned by σ1 , σ2 and σ3 therefore ∇φ = ∂σ is perpendicular to the yield surface.
Since the above applies for the plastic deformation increment, one can write:

∂φ (plastic deformation increment)

dεp = dλ (6.16)
∂σ
The direction of both vectors are equal, the magnitude of the plastic deformation increment
is scaled with the scalar dλ which is similar to its one-dimensional counterpart.
58 Plasticity theory

Figure 6.9: non-perpendicular deformation increment

∂φ
It can be shown that ∇φ = ∂σ with two examples. First, assume a circular function
∂φ
2 2 2
φ = x + y − R . The gradient of this function is ∂x = ( ∂φ ∂φ
∂x , ∂y ) = (2x, 2y). The latter
vector is perpendicular to the circular function φ.
Second, assume a vertical line at a distance c from the origin φ = x − c = 0. The gradient of
∂φ
this function is ∂x = ( ∂φ ∂φ
∂x , ∂y ) = (1, 0). The latter vector is perpendicular to the vertical line.

The principle of normality can be used to derive the demand for convexity of the yield surface.
For a visual explanation refer to figure 6.10. Equation (6.14) must be valid for any arbitrary
stress in the elastic domain. If the yield surface is non-convex though this might conflict
with the normality rule where the resulting stress increment would be instable for a plastic
deformation direction that is normal to the yield surface. As a result it is required that the
yield surface is a convex function.

Figure 6.10: Bending points in yield surface

Note that this only applies to regular points. For singular points, in case of an intersection of
two yield surfaces, see Figure 6.11, the convexity of the yield surface is guaranteed if φ1 = 0
6.4 Elastic - ideal plastic constitutive equations 59

and φ2 = 0 are both convex.

(Interesection of two yield surfaces)

Figure 6.11: Singular points

The plastic deformation increment in the singular point is defined as:

∂φ1 ∂φ2
dεp = dλ1 + dλ2 (6.17)
∂σ ∂σ

6.4 Elastic - ideal plastic constitutive equations

In the previous section the requirements from a multi-axial yield surface as well as the direction
of plastic flow is derived. A similar treatment as in the one-dimensional scenario can now be
applied in obtaining consitutive relations that describe the mechanical behavior in terms of
classical plasticity formulations.
The starting point is the decomposition of total strain rate into an elastic and a plastic part
and Hooke’s law:
ε̇ = ε̇e + ε̇p (total strain rate) (6.18)
σ̇ = E : ε̇e = E : (ε̇ − ε̇p ) (Hooke's law) (6.19)
The plastic strain increment can be written as:

∂φ
ε̇p = λ̇ (6.20)
∂σ
Substitution of equation (6.20) in (6.19) gives:
 
∂φ
σ̇ = E : ε̇ − λ̇ (6.21)
∂σ
Kuhn-Tucker conditions can be written in tensorial form as follows:

λ̇ ≥ 0
φ(σ) ≤ 0
λ̇ φ(σ) = 0 (6.22)
60 Plasticity theory

Consistency condition becomes:

λ̇ φ̇(σ) = 0 (6.23)
Using (6.23) and (6.19):
∂φ
φ̇ = : σ̇ = 0
∂σ
  (6.24)
∂φ ∂φ
⇒ : E : ε̇ − λ̇ =0
∂σ ∂σ

λ̇ can be now written as:

Result is a second order tensor
∂φ
: E : ε̇
λ̇ = ∂σ (6.25)
∂φ ∂φ
:E:
∂σ ∂σ
Substituting the expression for λ̇ in (6.21), gives the constitutive equation for elastic ideal
plastic material behavior in its most general form:

∂φ ∂φ
 
E: :E
σ̇ = E −
 ∂σ ∂σ 
: ε̇ = (E − Y) : ε̇ (6.26)
∂φ ∂φ 
:E:
∂σ ∂σ

6.5 Multi-axial yield functions

An arbitrary multi-axial yield surface, normality rule and convexity requirement are intro-
duced previously. In this section, some specific yield criteria will be treated, i.e. Tresca, Von
Mises and Hill. The first two are used to describe isotropic material behavior whereas the
Hill yield criterion can describe anisotropy. All yield criteria have to fulfill the demands of
the previous subsection:

• the yield surface reduces to the yield stress in case of a uniaxial stress state,

• the yield surface has to be convex.

For metals the following are also observed in experiments: hydrostatic pressure does not cause
plastic deformation and the plastic deformation is incompressible. The first condition implies
that the yield surface does not depend on pressure. The second condition can be expressed
as:

εpii = 0 =⇒ ε̇pii = 0
With (6.16) follows:

∂φ
ε̇p = λ̇ (6.27)
∂σ
In general λ̇ 6= 0 , which gives:
trace
∂φ
=0 (6.28)
∂σii
6.5 Multi-axial yield functions 61

6.5.1 Tresca yield criterion

In the yield criterion of Tresca, plastic deformation occurs when the maximum occurring shear
stress exceeds a certain value. This value is related to the yield stress via Mohr’s circles, which
is explained below. Starting point of Figure 6.12 is a plane stress state; the third principal
stress σ3 is set to 0. Extension to three-dimensional stress states is straightforward. Following
Mohr’s circles, the maximum shear stress equals the radius of the largest circle of Mohr. In
case of a uni-axial stress state the principal stress equals the yield stress. The accompanying
shear stress will be half the principal stress and thus half the yield stress. Extension to a

(Plane stress state)

Figure 6.12: Creation of Tresca yield surface

multi-axial stress state means that the diameter of the largest circle of Mohr equals the yield
stress. The diameters of both other circles have to be smaller or equal to the yield stress. In
other words, the Tresca yield criterion is nothing more than:
|τmax | = 12 σy (6.29)
In three dimensions, the Tresca yield criterion has the following form:
φ = (σ1 − σ2 )2 − σy2 (σ2 − σ3 )2 − σy2 (σ3 − σ1 )2 − σy2 = 0
   
(6.30)
With this equation it is easy to prove that the Tresca yield criterion does not depend on the
pressure. Since the pressure does not influence the yield criterion, the yield criterion is an
62 Plasticity theory

infinite prismatic, along the axis in {1,1,1}-direction.

   


 σ 1 
 1

 

σ2 = ξ 1 (6.31)

   
σ   1
 
3

The cross section of the plane perpendicular to the {1,1,1,}-direction is always identical. The
plane σ1 + σ2 + σ3 = 0 is also known as the π-plane. The cross section of the π-plane with
the Tresca yield surface is represented in figure 6.13.

Figure 6.13: Tresca yield surface in the π-plane

6.5.2 Von Mises yield criterion

In metals, a uniform pressure in all directions only results in elastic deformations and no
plastic deformations. For plastic deformations, it is necessary to have a deviating stress state,
not only a pressure (p). The deviation from a uniform pressure is called the deviatoric stress.
The principal deviatoric stresses are defined by:
s1 = σ1 − 13 (σ1 + σ2 + σ3 ) = σ1 − p
s2 = σ2 − 13 (σ1 + σ2 + σ3 ) = σ2 − p
s3 = σ3 − 13 (σ1 + σ2 + σ3 ) = σ3 − p
with
p = 13 (σ1 + σ2 + σ3 )
In previous courses, the yield criterion of Von Mises is already discussed in terms of principal
stresses; plastic deformation occurs if the sum of squares of the principal deviatoric stresses
reaches a certain critical value:
s21 + s22 + s23 = 2k2 (6.32)
The principal stresses and principal deviatoric stresses are components of tensors with respect
to a coordinate system of which the coordinate axes coincide with the principal axes:
   
σ1 0 0 S1 0 0
   
σij = 
 0 σ 2 0  , sij =  0 S2 0 
   (6.33)
0 0 σ3 0 0 S3
6.5 Multi-axial yield functions 63

The relation between deviatoric stresses and the total stresses with respect to an arbitrary
reference frame, is given by:

sij = σij − 13 σkk δij = σij − pδij (6.34)

in which p = 13 σkk = 13 (σ11 + σ22 + σ33 ). The elastic boundary can be represented by:

sij sij = 2k2 (6.35)

The term sij sij is invariant, it does not change if the tensor components are transformed to
another reference frame (because it is a scalar quantity). Therefore, Equation (6.35) holds for
every arbitrary Cartesian coordinate system. The equivalent stress according to Von Mises is
defined by:
q
σ = 32 sij sij (6.36)

The factor 32 is included in order to achieve that in a uni-axial tensile test the uni-axial stress is
equal to the equivalent stress (see section 6.5.6). This section also shows the relation between
k and σy . Hence, the Von Mises yield surface in terms of deviatoric stresses is defined as:

φ = 23 s : s − σy2 = 0 (6.37)

In case of using stresses instead of deviatoric stresses, the Von Mises yield surface reads,
according Huber and Hencky:

(σx − σy )2 + (σy − σz )2 + (σz − σx )2 + 6τxy

2 2
+ 6τyz 2
+ 6τzx − 2σy2 = 0 (6.38)

The Von Mises yield surface is also pressure independent, and as a result, will be a prismatic
in the three dimensional stress space. Writing the Huber and Hencky equation in principal
stresses, the equation describes a cylinder with its axis in the {1,1,1}-direction. The cross
section of the yield surface with the σ1 − σ2 - plane is elliptical, the cross section with the
π-plane is a circle, see Figure 6.14. In three dimensions, the Von Mises yield surface is
represented in Figure 6.15. Besides, it is simple to prove that the Von Mises yield surface
fulfills all demands, as described in the introduction.

Figure 6.14: Von Mises yield surface

64 Plasticity theory

s1
von Mises

p-plane

s2 s3

Tresca

Figure 6.15: Von Mises yield surface in 3 dimensions

6.5.3 Comparison of the Tresca and Von Mises yield surface

Since a yield surface has to be convex and has to degenerate to the yield stress in case of a
uni-axial stress state, the Tresca yield surface is a lower bound. Also, an upper bound with 6
sides can be defined, see Figure 6.16. Von Mises yield surface lies in between the bounds. A

s1 A: inner bound (Tresca)

B B: outer bound
C: von Mises
A

s2 s3

Figure 6.16: The Von Mises criterion, bounded by the lower bound and the upper bound

disadvantage of the Tresca yield criterion compared to the Von Mises yield criterion, is that
it contains singular points.

6.5.4 The Hill yield criterion

The Tresca yield criterion and the Von Mises yield criterion cannot be used in case of
anisotropic material behavior. The Hill yield criterion is an extension of the Von Mises
yield criterion, which is able to capture anisotropy. However, 6 extra material constants have
to be determined experimentally as input for the Hill yield criterion:

φ = F (σx − σy )2 + G(σy − σz )2 + H(σz − σx )2 + 2Lτxy

2 2
+ 2M τyz 2
+ 2N τzx − 32 (F + G + H)σy2 = 0
(6.39)
In general, sheet metals show anisotropy. If the in-plane material properties are similar, some
terms in the yield criterion vanish. This material behavior is called planar isotropic. In the
6.5 Multi-axial yield functions 65

principal stress space, the Hill yield criterion is then represented as follows:

φ = R(σ1 − σ2 )2 + (σ2 − σ3 )2 + (σ3 − σ1 )2 − (R + 1)σy2 = 0 (6.40)

where R is defined as:

ε2
R= (6.41)
ε3
In case of a tensile test, ε2 is the in-plane strain perpendicular to the tensile direction, and

σ2

σ1

R=1/2
R=1
R=2
R=3

Figure 6.17: Hill 48 yield functions for various R-values

ε3 is the strain in thickness direction. The Hill yield criterion degenerates to the Von Mises
yield criterion if R equals 1 (isotropic material behavior). Besides, it is simple to prove that
the Hill yield surface fulfills all demands, as described in the introduction. Examples of the
Hill yield surface for various R-values are given in Figure 6.17.

6.5.5 Non-associated flow behavior

Not all materials show pressure independent flow behavior (i.e. equation (6.28) does not apply
anymore). This applies to the material models used in soil mechanics as well as polymers
which show pressure dependent flow behavior. In these cases, the yield surface has the shape
of a parabola or cone in the principal stress space. For example, the behavior of sand can be
characterized by:

• yield under tensile stresses;

• an increasing yield stress with an increasing pressure.

Ideally, the material is assumed incompressible. In two dimensions, the material behavior can
be characterized by a material which shows Coulomb friction:
The stress-strain relation is:

τ = Gγ (|τ | < τmax ) (elastic behavior till the yield stress)

|τ | = −µσn (|τ | = τmax ) (pressure dependent flow behavior)
If the principle of normality should apply here, then the plastic deformation rate is normal
to the yield surface. Lets focus on the slant part of Figure 6.19. Normality means that the
66 Plasticity theory

Figure 6.18: Coulomb-model

plastic deformation rate has a component in positive direction, which also means that the
plastic deformation rate will be positive and hence, that the volume will increase. In other
words, applying a pressure should lead to an increase in volume, which is strongly arguable.
Therefore plastic deformation is only allowed to occur in τ -direction, see Figure 6.19.

Figure 6.19: No normality

In general, for this kind of material behavior, it is assumed that the plastic deformation rate
is determined by a potential which not equals the yield surface (this in contrast with the
former material behavior where both coincide). This is called non-associated flow.

6.5.6 Constitutive equations for elastic-ideal plastic material, using the

Von Mises yield criterion (Do we have to know this for the exam?)
In the following formulations certain assumptions will be made: linear isotropic elastic re-
sponse, Von Mises yield criterion, normality rule, no hardening.
Von Mises yield criterion is introduced previously as:

sij sij − 2k 2 ≤ 0

with sij the deviatoric stress:

sij = σij − 13 σkk δij (6.42)
6.5 Multi-axial yield functions 67

and with p = 13 σkk :

sij = σij − pδij (6.43)
However the parameter k, which is a material constant, is not known yet. One of the condi-
tions of the yield surfaces has been introduced as the reduction to a uniaxial yield stress for
uniaxial loading.
In case of a tensile test and ideal plasticity, σ11 = σy and σ22 = σ33 = σ12 = σ13 = σ23 = 0.
Hence:

s11 = σ11 − 31 σ11 = 23 σ11 = 32 σy

s22 = s33 = − 13 σ11 (= − 21 s11 )

sij = 0 if i 6= j
Furthermore:
sij sij = s211 + (− 12 )2 s211 + (− 12 )2 s211 = 32 ( 23 σ11 )2
in other words:
sij sij = 2k2 = 23 σy2 (6.44)
hence:
σy
k=√ (6.45)
3
The yield surface is obtained as in equation (6.20) using the above:

φ = 32 s : s − σy2 ≤ 0 (6.46)

Often an alternative version of the yield surface is used based on the equivalent stress concept:
VM
φ∗ = σeq − σy ≤ 0 (6.47)
q
with VM
σeq = 32 s : s.
With this result at hand now the plastic strain can be related using normality to the yield
surface using the previously obtained result:

∂φ
ε̇p = λ̇
∂σ
and consequently using (6.46):
ε̇p = 3 λ̇ s (6.48)
Using the consistency condition λ̇ can be determined as:

1 s : E : ε̇
λ̇ = (6.49)
3s:E:s
Based on the isotropic elasticity assumption this result can be further reduced to:
1 n
λ̇ = √ : ε̇ (6.50)
3 s:s
√
where n is the direction of the deviatoric stress: n = s/ s : s.
68 Plasticity theory

 
E : ss : E
σ̇ = E− : ε̇ = (E − 2G nn) : ε̇ (6.51)
s:E:s
In this case, the deformation history is given by the logarithm of the strain ε. By introducing
a scalar deformation value, the equivalent plastic strain, this equation can be generalized to
an arbitrary deformation state (like the equivalent stress)

Zt q
2 p p
ε= 3 ε̇ij ε̇ij dt (6.52)
t0

Note that the equivalent plastic strain cannot decrease. The factor 23 is introduced to make
the equivalent plastic strain equal to the logarithmic strain in case of a uni-axial tensile test
(ε̇11 = ε̇11 , ε̇22 = ε̇33 = − 21 ε̇11 ) . Proof:

Zt q Zt q Zt q Zt
2 2 2 2 1
ε= 3 ε̇ij ε̇ij dt = 3 (ε̇11 + ε̇222 + ε̇233 ) dt = 3 (1 + 4 + 14 )ε̇211 dt = ε̇11 dt = ε
t0 t0 t0 t0
(6.53)

6.6 Constitutive equations for a hardening material

A material with arbitrary hardening behavior is assumed. Time dependent effects like creep
and relaxation are not taken into account. The yield surface can then be written as:

φ(σij , βij ) = 0 (6.54)

component of a second order tensor
where tensor β denotes hardening. The constitutive equations for elastic non-ideal plastic
material behavior are derived similar to the strategy followed in section 6.4. Following Hooke’s
law and the expression for the plastic strain increment, equation (6.21) was derived:
 
∂φ
σ̇ij = Eijkl ε̇kl − λ̇
∂σkl

The consistency condition is used to find an expression for λ̇:

∂φ ∂φ
φ̇ = σ̇ij + β̇ij = 0 (6.55)
∂σij ∂βij

If no plastic deformation occurs, (see (6.16)):

λ̇ = 0

as well as:
β̇ij = 0
It is assumed that the relation between β̇ and ε̇p can be written as:

β̇ij = Pijkl ε̇pkl (6.56)

6.6 Constitutive equations for a hardening material 69

with Pijkl a 4th order tensor. Using (6.20) gives:

∂φ
β̇ij = Pijkl λ̇ (6.57)
∂σkl
and
∂φ ∂φ ∂φ
β̇ij = Pijkl λ̇ (6.58)
∂βij ∂βij ∂σkl
∂φ ∂φ
The term ∂βij Pijkl ∂σkl is a scalar and is called ’f’ from now on. So:

∂φ
β̇ij = f (σij , βij , T, . . .)λ̇ (6.59)
∂βij
Now equation 6.55 can be written as:
∂φ
φ̇ = σ̇ij + f λ̇ = 0 (6.60)
∂σij
In the subsequent subsections the function f will be derived for both isotropic and kinematic
hardening. Substitution of (6.21) in (6.60) yields the next expression:
 
∂φ ∂φ
φ̇ = Eijkl ε̇kl − λ̇ + f λ̇ = 0 (6.61)
∂σij ∂σkl
which gives the following expression for λ̇:
∂φ
Eijkl ε̇kl
∂σij
λ̇ = (6.62)
∂φ ∂φ
Epqrs −f
∂σpq ∂σrs
The general constitutive equations for elastic non-ideal plastic material behavior can be de-
rived by combining equations (6.62) and (6.21).
 
∂φ ∂φ
Eijmn Etukl 
 ∂σmn ∂σtu
σ̇ij = E
 ijkl −  ε̇kl (6.63)
∂φ ∂φ 
Epqrs −f
∂σpq ∂σrs
Or in tensor notation:
∂φ ∂φ
 
E:
:E
σ̇ = E −
 ∂σ ∂σ  : ε̇

(6.64)
∂φ ∂φ
:E: −f
∂σ ∂σ
A brief version of the expression is:
σ̇ = (E − Y ∗ ) : ε̇ = (E − (1 − h)Y) : ε̇ (6.65)
with:
∂φ ∂φ
E: :E f
Y= ∂σ ∂σ and h =
∂φ ∂φ ∂φ ∂φ
:E: f− :E:
∂σ ∂σ ∂σ ∂σ
with h a scalar which depends on the hardening. The yield behavior is expressed by the
tensor Y in (6.65).
70 Plasticity theory

6.7 Hardening
6.7.1 Introduction
Previously a constitutive equation was derived for elastic ideal plastic material behavior.
Since hardening does not occur, the yield surface will not change. In practice, most materials
show hardening. As a result the yield surface will change in size or position in the stress
space, see Figure 6.20. Two extreme types of hardening can be distinguished:

Figure 6.20: Hardening

1. Isotropic hardening The shape and position of the yield surface remain the same, the
size of the yield surface uniformly increases, see Figure 6.21.

Figure 6.21: Isotropic hardening

2. Kinematic hardening The shape and the size of the yield surface remain the same, the
position of the yield surface shifts in the stress space, see Figure 6.22.

The difference between both types of hardening can also be clearly spotted when looking at
the stress-strain curve, see Figure 6.23. In case of kinematic hardening the elastic part during
unloading equals the initial elastic part. In case of isotropic hardening, the elastic part has
increased during unloading. The hardening of a material is not purely kinematic or isotropic.
Some parts of the material behavior are better described with the isotropic hardening and
6.7 Hardening 71

Figure 6.22: Kinematic hardening

Figure 6.23: Stress-strain curves for both types of hardening

72 Plasticity theory

vice versa, see last subsection on hardening. First, a constitutive equation will be derived to
take into account hardening. Then the general form of this equation will be applied to both
types of hardening.

6.7.2 Isotropic hardening

The Von Mises yield criterion will be used as a starting point for isotropic hardening:

φ = 32 sij sij − σy2 (εp ) = 0 (6.66)

It is assumed that the yield stress is not a constant, but depends on the equivalent plastic
strain (which is a measure for the plastic strain). Note: in the general form (6.54), β is
defined as the hardening tensor. Equation (6.66) is a special case of (6.54) when it is assumed
that βij = εp δij . The relation between the yield stress and the equivalent plastic strain is
p
determined by a tensile test. During this tensile test, the stress σ11 and the strain ε11 (so also
p
the plastic strain ε11 ) in the tensile direction are measured. It is defined that the equivalent
plastic strain equals the plastic strain in the tensile direction εp = εp11 . To determine the
unknown function f in equation (6.60), the derivative of the yield surface is taken (6.66):

dσy p
φ̇ = 3sij ṡij − 2σy ε̇ = 0 (6.67)
dεp
The first term in (6.60) equals the first term in (6.67) since in case of Von Mises yield behavior:

∂φ
σ̇ij = 3sij σ̇ij = 3sij ṡij (6.68)
∂σij

∂φ ∂φ ∂skl ∂(σkl − 31 σpp δkl )

= = 3skl = 3skl (δik δjl − 13 δip δjp δkl ) = 3sij − 3skk δip δjp = 3sij
∂σij ∂skl ∂σij ∂σij
(6.69)

An expression for f is found, using equations (6.67), (6.68) and (6.60):

dσy p
f λ̇ = −2σy ε̇ (6.70)
dεp

This expression can be elaborated further, using the relation between ε̇p and λ̇, equation
(6.20):
ε̇pij = 3λ̇sij = µ̇sij (6.71)
First a relation between the equivalent plastic strain rate ε̇p and the plastic strain rate tensor
ε̇pij has to be determined. This relation is defined by the tensile test. The plastic volume
changes are zero, see subsection 6.5. So, in case of a tensile test, the following relations can
be determined:
ε̇p22 = ε̇p33 = − 12 ε̇p11 (6.72)
The following norm is used to define the equivalent plastic strain rate:
q
ε̇p = const ε̇pij = const ε̇pij ε̇pij
6.7 Hardening 73

Filling in (6.72) gives:

q q
ε̇p = constε̇11 1+ 1
4 + 1
4 = constε̇11
3
2

It was already defined that εp = εp11 , so ε̇p = p

ε̇11 , which means
that the constant can be
determined. The relation between the equivalent plastic strain rate and the plastic strain
rate tensor can now be defined as:
q q
ε̇p = 23 ε̇pij = 23 ε̇pij ε̇pij (6.73)

Filling in the Prandtl-Reusz equation (6.71) gives a relation between ε̇p and λ̇.
q
ε̇p = 3λ̇ 23 sij sij (6.74)
Again, the expression under the square root is found, using the tensile test. The deviatoric
stress components are written as:
s11 = σ11 − 31 σ11 = 23 σ11 = 32 σy (6.75)
s22 = s33 = − 13 σ11 = − 31 σy (6.76)
Filling in the above stress components in (6.74) gives:
ε̇p = 2λ̇σy (6.77)
The expression for f in case of isotropic hardening, using the Von Mises yield criterion, is
found by substituting equation (6.77) in (6.70):
dσy
f = −4σy2 (6.78)
dεp
Finally, note that the equivalent plastic strain is determined as follows (and not calculated
directly from the plastic strains):
Zt1 Zt1 q
p p 2 p p
ε = ε̇ dt = 3 ε̇ij ε̇ij dt (6.79)
t0 t0
This strategy is followed since the total strain can be zero during cyclic loading, while the
equivalent plastic strain can not decrease by definition.

Nadai model
For many steels, a good approximation of a large part of the true stress–strain curve can be
obtained with a power law, e.g. the Ludwik–Nadai or Nadai relation1
σy = Cεn (6.80)
C is the coefficient of hardening and n is the exponent of hardening.ε should formally be the
plastic strain in the tensile direction, but for large deformations the elastic part is usually
neglected. To determine the yield stress according to the Nadai hardening law, the following
equation is used:2

σy = σ0 + C (εeq + ε0 )n (6.81)
with σ0 the offset yield stress, ε0 the offset plastic strain, and C and n are power law variables.
1
Also called the ‘Hollomon’ relation.
2
Also called the ‘Swift’ relation.
74 Plasticity theory

Figure 6.24: Nadai hardening

Voce model
If the hardening rate is plotted as a function of the stress, often a linearly decreasing relation
can be observed. This part of the hardening curve is often denoted as Stage III hardening.
A relation that represents this behaviour is as follows:

dσf  σ 
θ= = θ0 1 − sat (6.82)
dε σ
where σ sat is the saturation stress. The resulting flow stress as a function of the strain is
written as   
ε
σy = σy0 + ∆σ 1 − exp − (6.83)
ε0
where ∆σ = σ sat − σy0 and ε0 = σ sat /θ0 . This model is known as the Voce model.
The Voce hardening law is evaluated as:
(
σy = σy0   if εeq ≤ ε1

εeq −ε1 (6.84)
σy0 + ∆σ 1 − exp − ε0 if εeq > ε1

with σy0 the initial yield stress and ∆σ the increase of the yield stress until saturation is
reached.ε0 is the plastic strain after the yield point plateau was reached; the plastic strain
parameter ε1 indicates the length of the yield plateau.

Effects of temperature and strain rate

Increasing the temperature, decreases the yield strength of many metals. It also affects their
ductility, but not in a simple way. With increasing T /Tm (Tm is the melting temperature)
the ductility increases or decreases depending on the crystal structure of the metal. E.g.
the ductility of Al–Mg alloys increases significantly, while for pure aluminum the effect is
marginal.
6.7 Hardening 75

Some metals with small grains, at temperatures above 0.5 Tm and at specific strain rates
exhibit so-called super-plastic behavior. Under these conditions they can reach engineering
strains up to 1000 %. E.g. an aluminum alloy 2004, when deformed at 430 ◦ C and a strain
rate of 2 × 10−3 s−1 , gives et = 800 % at a flow stress of 10 MPa. This same alloy when heat
treated to the T6 condition has a tensile strength of 450 MPa and is used for the manufacture
of a range of aircraft components [8, 127].
The superplastic deformation behavior is often described by the stress–strain-rate relation

σ = kε̇m (6.85)

and, in the superplastic range, 0.5 < m < 1. The parameter m is itself a function of the
strain rate (among others) and it is the task of the manufacturer of products to control the
process as to assure that the value of m is near a maximum. For high values of m localization
is postponed and high strains can be reached.
Combining, for lower temperatures, the strain hardening and strain-rate hardening equations,
we get
σ = Kεn ε̇m (6.86)

On a microscopic scale, plastic deformation in metals is caused by dislocation movement. The

physics of dislocation movement suggests a thermally activated process. The relation between
strain rate and temperature can then be derived from statistical mechanics. This approach
is the basis of the Zener–Hollomon parameter [10] defined as
 
U
Z = ε̇ exp (6.87)
kT

where U is an activation energy k is Boltzmann’s constant and T is the absolute temperature.

This parameter can be used to describe the combined effects of strain rate and temperature. In
Equation (6.86), the strain rate ε̇ must be replaced by Z to account for variable temperatures.

Physically based hardening models

In physically based hardening models, the evolution of the flow stress is predicted by consider-
ing the physical mechanisms of plastic deformation. Ideally, the influence of work hardening,
strain rate and temperature, and their interactions, are included.
Plastic deformation in f.c.c. and b.c.c. metals is mainly the result of dislocation movement on
the active slip systems. Hence, the resistance to dislocation movement determines the yield
stress. In a physically based model the yield stress is often decomposed into a strain and
strain rate independent stress σ0 , a dynamic stress σ ∗ that depends on the strain rate and
temperature and a term σw that incorporates the work hardening:

σf = σ0 + σ ∗ (ε̇, T ) + σw (ρ, T ) (6.88)

The dynamic stress σ ∗ is the driving force that is required to move dislocations through the
atomic lattice at a certain velocity. The velocity of a dislocation depends on the rate at
which atoms change position in the atomic lattice. This rate is described with a probabilistic
approach in statistical thermodynamics. An approximation of the resulting implicit relation
76 Plasticity theory

for σ ∗ for relatively low temperatures yields the description of Krabiell and Dahl [6] that is
often used in simulation programs:

0  if ε̇ < ε̇0 exp(−∆G0 /kT )

 p
∗ ∗ kT ε̇
σ (ε̇, T ) = σ0 1 + ∆G0 ln ε̇0 if ε̇0 exp(−∆G0 /kT ) < ε̇ < ε̇0 (6.89)


σ ∗ if ε̇ > ε̇
0 0

where ∆G0 is an activation energy.

The dynamic stress causes a translation of the stress–strain curves. For increasing temper-
atures, the dynamic stress decreases and vanishes at an absolute temperature of about half
the melting temperature.

The work hardening part of the model, σw , takes the evolution of the micro-structure into
account. A relatively simple one-parameter model only considers the dislocation density ρ
to be responsible for the hardening. The relation between the dislocation density and σw is
given by the Taylor equation:
√
σw = αGb ρ (6.90)
where α is a scaling parameter of order 1 and G is the shear modulus. This equation is sup-
ported by a dimensional analysis and numerous experiments and is used in a large number of
physically based models. The essential part in these models is the evolution of the dislocation
density ρ. Several models are available, and the reader is referred to literature.

Equivalent stress and strain

The hardening relations are usually derived from one-dimensional experiments. For the anal-
ysis of sheet metal forming, the relations should be expanded to at least two-dimensional
deformations.
For many metals it is found that for proportional loading the hardening is a function of the
(plastic) work done per unit volume. We want to relate the work done in a uniaxial stress
state to the work done in a multi-axial stress state. For this purpose we define an equivalent
plastic strain κ and find for the plastic work per unit volume:

σy dκ = σ · dε (6.91)

For a uniaxial loading process κ is taken equivalent to the uniaxial plastic strain and con-
sequently, the uniaxial stress σun is equal to the flow stress σf . The hardening found in a
uniaxial test is presented as the uniaxial stress σun as function of the uniaxial strain εun .
The equivalent stress σy in the numerical description of the yield functions is now considered
to vary in the same way depending on κ. For proportional loading paths, the ratio between
components in the stress tensor remains constant (α = σ2 /σ1 ), as well as the ratio between
components in the strain tensor (β = ε2 /ε1 ). From this consideration, relations can be de-
rived for the expression of any component of the stress tensor as function of any component
of the strain tensor. For the Nadai relation this yields simply:

σy = Cκn ⇒ σij = C ′ εnij (6.92)

where
C ′ = Cf (β, n) (6.93)
6.7 Hardening 77

6.7.3 Kinematic hardening

In case of kinematic hardening, the center of the yield surface f moves through the stress
space, see Figure 6.25. The size and shape of the yield surface do not change (Prager-Ziegler).
The Von Mises yield surface can be written as, using Prager-Ziegler (equation of a circle with

Figure 6.25: Kinematic hardening

moved center):
2
φ = 32 (sij − αij )(sij − αij ) − σy,0 =0 (6.94)
Now the yield stress is constant since the size of the yield surface does not change. To
determine the unknown in equation (6.60), equation (6.94) is differentiated first.

φ̇ = 3(sij − αij )ṡij − 3(sij − αij )α̇ij = 0 (6.95)

According Prager, the velocity of the center can be written as:

α̇ij = h∗ ε̇pij (6.96)

The displacement of the yield surface in the stress state has the same direction as the plastic
deformation. The first term in (6.60) equals the first term in (6.95), using the Von Mises
yield surface and the assumption of Prager (6.96):

∂φ
σ̇ij = 3(sij − αij )σ̇ij = 3(sij − αij )ṡij (6.97)
∂σij

A comparison of the second term in (6.60) and the second term in (6.95) gives an expression
for f :
f λ̇ = −3(sij − αij )α̇ij (6.98)
Using the assumption of Prager (6.96) and the principle of normality gives:

∂φ
f λ̇ = −3(sij − αij )h∗ λ̇
∂σij
78 Plasticity theory

The derivative of the yield surface to the stress is given in (6.97). This gives:

f = −9h∗ (sij − αij )(sij − αij ) (6.99)

Substitution of (6.94) gives:

2
f = −6h∗ σy,0 (6.100)
with σy,0 the initial yield stress. The factor h∗ can be determined from (6.95) and (6.96),
using the tensile test:
α̇11 = ṡ11 = 23 σ̇11 = h∗ ε̇p11 (6.101)
Now the following relations for h∗ and f can be determined:

2 dσy
h∗ = 3 (6.102)
dεp
2 dσy
f = −4σy,0 (6.103)
dεp
In case of a bilinear approximation of the stress-strain curve, as given in Figure 6.26, the
factor h∗ is written as:
h∗ = 23 Et (6.104)

Figure 6.26: bilinear stress-strain curve

The model described above is a linear hardening model. The evolution of the backstress,
Equation (6.96), can be written in incremental form as:
∂φ
∆α = ∆λh∗ (6.105)
∂σ
Another well known kinematic hardening model is the Armstrong–Frederick hardening
model. The evolution of the backstress is described by:
 
∗ ∂φ
∆α = ∆λ h − hl α (6.106)
∂σ
The second term in the brackets gives a smooth transition from the elastic to the plastic
region.
6.8 Elastic-viscoplastic material 79

6.8 Elastic-viscoplastic material (Do we have to know this for the exam?)

6.8.1 Introduction
Creep is the occurrence of plastic deformation during prolonged constant loading. On the
other hand, relaxation means that at a constant deformation, the stress decreases with in-
creasing time. Shape recovery after load removal is related to both phenomena. Finally, there
is a relation with strain-rate dependent material behavior. The material behavior of metals
is almost independent of the strain-rate, when this is not too large. Therefore metals do have
a so called static stress-strain curve. This behavior can be described with a Bingham model
with a spring in series. However, plastic strain occurs during elongated loading. This behavior
can be modeled with a Maxwell model (area of secondary creep); this model does not describe
the plastic behavior for stresses higher than the yield stress. In case of large strain rates, the
influence of the strain rate on the material behavior of metals can not be neglected, like in the
rolling process. An increase of the strain rate yields an increase of the stress. For each fixed
strain rate, the stress-strain curve is higher than the static stress-strain curve, but with the
same shape. Polymers are sensitive for creep and are strain rate dependent for all strain rates.
Polymers can be described with the same model in the both the creep area as in the plastic
deformation area. The model looks like the Maxwell model, but has a non-linear damper. In
general, the logarithm of the strain is proportional with the stress. Strain rate dependency
in metals and polymers and creep in polymers is called viscoplastic material behavior and is
treated in subsection 6.8.3. Subsection 6.8.2 treats creep in metals.

6.8.2 Creep in metals

Creep behavior is determined by loading a specimen at a constant temperature and a constant
load. The strain is recorded as a function of the time, see Figure 6.27. The temperature must

Figure 6.27: Creep behavior

be as constant as possible since a small change in temperature can have a large effect on the
results. The next areas can be distinguished in Figure 6.27:

Primary creep (area I):

The strain rate decreases according a logarithmic function. At primary creep, slip occurs
along the favorable slip planes in the grains. The dislocation density will increase, hardening
occurs and consequently the strain rate decreases.
80 Plasticity theory

Secondary creep (area II) :

If the activation energy becomes high enough, dislocations will move along their obstacles and
secondary creep occurs. This high level of energy can be gained by sufficient high stresses or
temperatures. the mechanism of deformation is a combination of slip within the grains and
the in between grains.

Tertiary creep (area III):

The strain rate increases until fracture due to creep occurs. Inter-crystalline holes are formed
due to slip in the grains. In the tertiary area necking occurs, which means that it is not
possible anymore to keep the load constant.

The creep curves depend on the temperature and the stress. Creep is a time dependent
phenomenon; plastic deformation occurs while the yield stress is not exceeded. This is called
viscoelastic material behavior. When the yield stress is exceeded, the material behavior can
become both time and history dependent: plastic deformation is caused by plasticity and
creep. This phenomenon is called viscoplastic material behavior. We restrict ourselves to
the area with the constant strain rate, secondary creep. Starting point will be the Maxwell
model, see Figure 6.28.

Figure 6.28: Maxwell model

6.8.3 Visco-plastic material behavior

The strain tensor can be divided in a pure elastic and a pure viscoplastic part:

εtot = εe + εvp (6.107)

Visco-plastic behavior is characterized by a change in strain at a constant stress. The strain

rate is a function of the creep strain, stress and temperature. We assume a constant temper-
ature, so:
ε̇vp = ε̇vp (σ, εvp ) (6.108)
6.8 Elastic-viscoplastic material 81

Besides, isotropic material behavior is assumed and it is assumed that the relation is also
valid for other stress states. The effective stress becomes then:
q 1
σeq = 32 (s : s) 2 (6.109)

and the equivalent viscoplastic strain:

1
q
ε̇vp = 2 vp
3 (ε̇ : ε̇vp ) 2 (6.110)

The direction of the plastic strain is defined by the deviatoric stress:

ε̇vp = λs (6.111)

Where the scalar λ can be written as:

ε̇vp
λ= 2 (6.112)
3 σeq

To obtain a stress strain relation, Hooke’s law is used in combination with equation 6.107:

σ = E : εe = E : (ε − εvp ) (6.113)

Since we only have an expression for the viscoplastic strain rate, the derivative of the above
equation is taken:
σ̇ = E : (ε̇ − ε̇vp ) (6.114)
substitution of equations 6.111 and 6.112 in 6.114 yields and expression for the viscoplastic
stress state:
ε̇vp
σ̇ = E : (ε̇tot − 2 s) (6.115)
3 σeq
Note that, contrary to the elastic-plastic model, the equivalent viscoplastic strain rate does
not have to be eliminated since the value is known through equation 6.110
Chapter 7

Finite Element Method

7.1 Introduction
Currently the most common numerical method that is used for solving mechanics problems
involving nonlinearity in material behavior is the Finite Element Method. In previous courses
this method was introduced (Introduction to the Finite Element Method) and further elabo-
rated (Numerical Methods in Mechanical Engineering). In the following, first, a brief summary
of the method is supplied with a focus on material nonlinearities. Afterwards the necessary
integration algorithms for the constitutive models that were introduced in Chapter 6 will be
given.

7.2 Virtual work/weak formulation of equilibrium

Both the virtual work method and method of weighted residuals are based on the equations
of equilibrium (4.8) and act as a means of integrating them. The result in both cases is
the integral form of the equations with weaker requirements. The derivation based on the
weighted residuals method is given below. However note that the weight functions are denoted
as virtual displacements (velocities).
Taking the dot product of the vector field resulting from the equilibrium equation with an
arbitrary vector field implies:
σij,j + Fi = 0 ⇔ (σij,j + Fi ) δvi = 0 ∀ δvi (7.1)
Integrating (7.1) over the whole volume gives:
Z
(σij,j + Fi )δvi dV = 0 ∀ δvi (7.2)
V

Integrating (7.2) by parts:

Z Z Z
(σij δvi ),j dV − σij δvi,j dV + Fi δvi dV = 0 ∀ δvi (7.3)
V V V

and using Gauss theorem we obtain:

Z Z Z
σij δvi nj dS − σij δvi,j dV + Fi δvi dV = 0 ∀δvi (7.4)
S V V

83
84 Finite Element Method

The first integral is nothing else than the contribution of surface forces Ti . In the second
integral, only the symmetric part (δDij ) of δvi,j contributes to work due to the required
symmetry of Cauchy stress:
δDij = 12 (δvi,j + δvj,i ) (7.5)

Equation 7.4 can be rewritten as:

Z Z Z
δW = σij δDij dV − Fi δvi dV − Ti δvi dS = 0 ∀δvi (7.6)
V V S

Equation (7.6) is usually referred to as the weak form of the equilibrium equation. However
there are alternative ways it can be interpreted. If the chosen arbitrary weight functions,
δvi , are considered to be virtual, kinematically admissible infinitesimal displacements then
the resulting equation is an expression of the virtual work theorem. In this case the first
term represents the internal virtual work and the last two represent the virtual work done
by external forces. And since the volume change would also be infinitesimal in this case this
equation can be solved using the linear finite element method.
In case of large deformations it is valuable to consider the weight functions, δvi , as kinemati-
cally admissible virtual velocities that are not necessarily infinitesimal. In this case equation
(7.6) represents the virtual power. It also implies that the derivatives are with respect to
current coordinates and that the volume is not constant. Accordingly, the tensor δDij can be
interpreted as the virtual rate of deformation.
In tensor notation equation (7.6) is given as:
Z Z Z
δD : σ dV = δv · F dV + δv · T dS ∀δv (7.7)
V V S

with δv a kinematically admissible virtual velocity field.

7.3 Finite element discretization

The finite element discretization of a domain V consists of the definition of subdomains with
a finite size Ve , the elements, and the respective interpolation functions. Evaluation of the
integrals in the weak form of mechanical equilibrium is performed over these elements. Later,
the contribution of each element is assembled into a large system representing the whole do-
main, yielding an approximate solution for the whole domain. An element is defined by means
of pre-defined points, the nodes. The interpolation functions provide continuous functions for
the displacements within an element. The interpolation functions relate the displacements at
any point in the interior of the element to the nodal point values. For isoparametric elements
the same interpolation is applied for both the geometry and the independent field variables:
n
X
x= N α xα (7.8)
α=1

n
X
ξ= N α ξα (7.9)
α=1
7.4 Iterative Newton–Raphson procedure 85

where N α is the time independent interpolation function for node α of the element with n
the number of nodes per element and ξ represents any independent field to be interpolated in
space. In the following the summation operator will be omitted and the summation convention
will be assumed.
In order to discretize the weak mechanical equilibrium, all continuous field variables and their
gradients are determined by the interpolation functions and the nodal values, leading to:

v = N α vα with N α = N α (x, y, z) (7.10)

or:        
v v 1  v 2  α
1 1 1  v1

 
 
  
  
 

      
v2 = N1 v21 + N2 v22 + . . . + Nα v2α (7.11)
       

 v    v1 
   v2 
   vα
 

3 3 3 3

The collection vα is a (multi-dimensional) vector of nodal velocities.

In a similar way, the virtual velocity can be written as:

δv = N α δvα (7.12)

The virtual rate of deformation tensor can be written as:

δDij = 12 (δvi,j + δvj,i ) = 12 N,jα δviα + N,iα δvjα = δvkα 21 δki N,jα + δkj N,iα = δvkα Bkij
α



(7.13)

δD = δvα · Bα (7.14)
where Bα is a third order tensor that relates the rate of deformation tensor to the nodal
velocity vector. Using Equation (7.12) and (7.14) in Equation (7.7) gives the discretized form
of the virtual work (power):
 
Z Z Z
δvα · Bα : σ dV = δvα ·  N α F dV + N α T dS 
V V S
α
⇒ δv · Fαint = δv α α
· Fext (7.15)

where Fαint is the internal force vector and Fαext the external force vector. Since Equation
(7.15) must hold for all kinematically admissible virtual velocity fields, the following equality
is obtained for the internal and external force vector:

Fαint = Fαext (7.16)

Note that this equation has to hold for every node.

7.4 Iterative Newton–Raphson procedure

The discrete equilibrium Equation (7.16) is generally non-linear and can therefore be solved
with the iterative Newton–Raphson procedure. A Newton–Raphson iteration consists of
two stages which are referred to as the predictor step and the corrector step, respectively.
In the predictor step an estimate of the solution is obtained based on linearization in the
current iteration. The corrector step provides the resulting reaction force for this iteration.
86 Finite Element Method

Repeatedly solving the linearized set of equations can lead to an approximation of the desired
solution, provided a converging process.
The starting point of the Newton–Raphson procedure is a displacement vector ∆u(0) . Next,
define:
∆u = ∆u(0) + ∆∆u (7.17)
Consequently, Equation (7.16) can be written as:

Fint (∆u) = Fint (∆u(0) + ∆∆u) = Fext (7.18)

Note that the superscript α, indicating the nodes, is omitted in this section to enhance
readability. Applying a Taylor series expansion to this equation yields:
∂Fint
Fint (∆u) = Fint (∆u(0) ) + | (0) · ∆∆u + O(∆∆u2 ) = Fext (7.19)
∂u ∆u
Neglecting the higher order terms O(∆∆u2 ) and rearranging gives:
∂Fint
K(0) · ∆∆u(1) = Fext − Fint (∆u(0) ) = R(0) with K(0) = | (0) (7.20)
∂u ∆u
where K(0) is the tangent stiffness matrix and R(0) is the residual force vector which van-
ishes when the exact solution is found. The unknown displacement increment ∆∆u(1) can be
solved using Equation (7.20). However, ∆u(0) + ∆∆u(1) is in general not equal to the exact
solution, since the higher order terms are omitted. It is likely that the new approximation
of the displacement vector, ∆u(1) = ∆u(0) + ∆∆u(1) , is a better approximation of the exact
solution than the old approximation ∆u(0) . The procedure is repeated again using ∆u(1) as
a new initial approximation. With the new approximation ∆u(1) , a new stiffness matrix K(1)
and internal force vector Fint (∆u(1) ) are determined. Subsequently, ∆∆u(2) is calculated
using Equation (7.20) (with K(1) and Fint (∆u(1) )) and the new displacement vector is set to
∆u(2) = ∆u(1) + ∆∆u(2) . This procedure is repeated until a user-defined accuracy is reached,
see algorithm below:

initialize
∆u(0) = 0
R(0) = Fext − Fint (∆u(0) )
k=0
WHILE (kR(k) k > tolerance) DO
determine K(k)
K(k) · ∆∆u(k+1) = R(k)
∆u(k+1) = ∆u(k) + ∆∆u(k+1)
determine Fint (∆u(k+1) )
R(k+1) = Fext − Fint (∆u(k+1) )
k := k + 1
END
In the above algorithm, it is assumed that Fext does not depend on ∆∆u. However, if Fext
7.4 Iterative Newton–Raphson procedure 87

depends on ∆∆u (e.g. in pressure loads), this can easily be incorporated in the algorithm.
In one dimension, the iterative Newton–Raphson process can be graphically represented as
shown in Figure 7.1. The exact solution is represented by û in this figure.
To avoid the creation of a new system of equations, the modified Newton–Raphson method
can be applied. In this method, the stiffness matrix in iterate k equals the stiffness matrix in
the first iteration (K(k) = K(0) ). While the convergence of this method will be slower, time
is saved by not computing the tangent stiffness matrices within each iteration, which in the
end can lead to a more economical method.
In general, the Newton–Raphson method converges, provided that the initial guess is close

Figure 7.1: Newton–Raphson iteration scheme.

to the exact solution. Divergence can occur for a too distant initial guess, see Figure 7.2.
If convergence cannot be obtained for a certain external force, this force can be subdivided
in several smaller loadsteps, and subsequently solving the system of equations for each of
the load substeps. The internal force vector Fint (∆u(k) ) and the stiffness matrix K(k) must

Figure 7.2: Diverging Newton–Raphson iteration process.

be determined in order to reach equilibrium using the iterative Newton–Raphson procedure.

88 Finite Element Method

To determine the internal force vector, according Equation (7.15), the stress at the end of
the time increment must be calculated. This is done by a stress update algorithm, which is
explained in Section 7.6. The derivation of the stiffness matrix is given in Section 7.5.

7.5 Determination of the stiffness matrix

The starting point for the derivation of the stiffness matrix is the time derivative of the
internal force vector:
dFαint dFαint duβ
Ḟαint = = = Kαβ · u̇β (7.21)
dt duβ dt

with α and β being nodal counters.

At this point we have established that due to material nonlinearities a nonlinear solution
scheme needs to be followed. However regarding geometrical nonlinearities, a choice can be
made at this stage. The difference between a geometrically linear calculation and a geomet-
rically non-linear one can be explained using the definition of the internal force vector, see
Equation (7.15):
Z
Fint = Bα : σ dV
α
(7.22)
V

Upon deformation, a start geometry (at time = t(n) and a final geometry (at time = t(n+1) ) can
be distinguished. Hence, the integral to compute the internal force vector can be evaluated
on the start geometry or on the final geometry. The geometrically linear equations are based
on the starting geometry whereas the geometrically nonlinear equations are based on the
final geometry. Since by definition equilibrium has to be fulfilled on the final geometry, the
geometrically nonlinear case represents a more realistic solution.
On the other hand for some applications it can be the case that nonlinear regime in the
material is reached before geometrical nonlinearities become significant, especially when an
incremental scheme is used with small time steps where rotations are also negligible. In the
following sections both cases are treated.

7.5.1 Geometrically linear case

The derivative of the internal force vector reads (see, Equation (7.15)):
Z
Ḟαint = Bα0 : σ̇ dV0 (7.23)
V0

where the subscript 0 refers to the start geometry. According to Equations (6.26) or (6.65)
the constitutive equation can be written as:

σ̇ = Ly : ε̇ (7.24)
7.5 Determination of the stiffness matrix 89

Substitution of Equation (7.24) in (7.23) and assuming that ε̇ = D gives:

Z
Ḟαint = Bα : Ly : D dV
V
Z
= Bα : Ly : Bβ · vβ dV
V
Z
= Bα : Ly : Bβ dV · vβ
V
= Kαβ · u̇β (7.25)
Hence, the stiffness matrix reads:
Z
αβ
K = Bα : Ly : Bβ dV (7.26)
V

7.5.2 Geometrically nonlinear case

In the geometrically nonlinear case, the internal force vector is defined in the current con-
figuration which changes in time. To determine the derivative of the internal force vector,
the integral is transformed into an integral with respect to the initial configuration where the
rate of the integral only consists of differentiating the terms within the integral with respect
to time. After determining the time derivatives in the initial configuration, the integral is
transformed to the current configuration.
First, without proof, it can be shown that when BαT is double contracted with a symmetric
second order tensor φ, the product degenerates to:
→
−
BαT : φ = 1∇N α : φ (7.27)
The transformation of the internal force vector from the current configuration to the initial
configuration (represented by subscript 0 ) reads:
→ α
− ∂N α ∂Xl →
−
Z Z Z  
α −T
Fint = 1∇N : σ dV = ei δij σjk dV = 1∇0 N α : σ · F J dV0 (7.28)
∂Xl ∂xk
V V V0

where J is the Jacobian (determinant of the deformation gradient F). The time derivative of
the internal force vector becomes:
→ α 
− →
−
Z  Z  
α −T −T
Ḟint = 1∇0 N : σ̇ · F J dV0 + 1∇0 N α : σ · Ḟ J dV0 (7.29)
V0 V0

where the time derivative of J is neglected since the volume changes are only elastic and
very small. Again without proof, it is stated that when the time derivative of the density is
incorporated in the constitutive equations, the time derivative of J cancels out with the extra
term concerning the density.
−T
An expression must be found for the time derivative of Ḟ :
0 = 1̇ = (F · F−1 )˙ = Ḟ · F−1 + F · Ḟ−1 ⇔
−1 −1 −1 −1 −1
Ḟ = −F · Ḟ · F = −F ·L·F·F = −F−1 · L ⇒
−T T −T
Ḟ = −L · F (7.30)
 90 Finite Element Method

Substitution of Equation (7.30) in (7.29) gives:

→
− →
−
Z   Z  
α α −T α T −T
(initial configuration) Ḟ = 1 ∇ 0 N : σ̇ · F J dV 0 − 1 ∇ 0 N : σ · L · F J dV0 (7.31)
int
V0 V0

The transformation of the time derivative of the internal force vector from the initial to the
current configuration is obtained through Equation (7.28):
→
− α →
−
Z Z  
α T
(current configuration) Ḟint = 1∇N : σ̇ dV − 1∇N α : σ · L dV (7.32)
V V

Next, the derivative of the Cauchy stress must be defined. In the geometrically nonlinear case
we have to account for large changes in geometry. First, a situation is considered in which
the strain remains elastic but the displacements and rotations are large. In that case Hooke’s
law still applies:
σ = E : εe (7.33)
Taking the material derivative gives:

σ̇ = E : ε̇e (7.34)

However, from chapter 3, it is known that the classical strain tensor is sensitive to rotations;
the strain tensor does not remain 0 upon rotation. In chapter 3 it is shown that the Green
Lagrange strain tensor is insensitive for rotations. However, if ε̇ is replaced by Ėgl in the
above equation, then the equation is not valid anymore in case of rotations. This can be
shown with the next example.
Assume a bar placed along the x-direction and load it with a small tensile stress σ 0 = σx0 .
Then the bar is rotated until it is aligned with the y-axis, without applying any further
deformation, i.e. Ėgl = 0. However, the stress components do change since these are related
to the x- and y-direction, see Figure 7.3. This means that, if we want to replace the strain

Figure 7.3: Stress change at rigid rotation

rate tensor with the Green Lagrange strain tensor, we have to define a stress measure which
is also independent on rotation. The Cauchy stress can be expressed in components referring
either to the global base vectors ei or local, co-rotating, base vectors gi :

σ = ei σij ej = gi τij gj (7.35)

7.6 Stress update 91

with gi = R · ei , and consequently:

T T
σ = R · ei τij ej · R = R · τ · R (7.36)

with τ = ei τij ej . The stress tensor τ is invariant under superimposed rigid rotation. Hence,
the derivative of Hooke’s law in invariant measures reads:

τ̇ = E : Ėgl (7.37)

with E the elasticity tensor and Ėgl the derivative of the Green Lagrange strain, as defined
by Equation (3.39):
T
Ėgl = F · D · F (7.38)
or, using polar decomposition and assuming small elastic deformations (U = I):
T T T
Ėgl = U · R · D · R · U ≈ R · D · R (7.39)

Next, the time derivative of the Cauchy stress in Equation (7.32) must be expressed as a func-
tion of the invariant stress rate to incorporate the constitutive equations in the formulation:
T T T
σ̇ = Ṙ · τ · R + R · τ̇ · R + R · τ · Ṙ (7.40)
T
The rotation tensor R is orthogonal (R = R−1 ), hence:
T T T
σ̇ = Ṙ · (R−1 · R) · τ · R + R · τ̇ · R + R · τ · (RT · R) · Ṙ
T T
= Ṙ · R−1 · σ + R · τ̇ · R + σ · R · Ṙ (7.41)

or:
T T T
R · τ̇ · R = σ̇ − Ṙ · R · σ − σ · R · Ṙ (7.42)
The above rate formulation is known as the Green-Naghdi rate. The Green-Naghdi rate equals
0 in case of rigid rotation (τ̇ = 0). In fact, the Green-Naghdi rate is a material derivative in
a co-rotating reference configuration.
Summarizing, the derivative of the Cauchy stress tensor can be expressed as a function of
the Cauchy stress tensor and the derivative of the invariant stress tensor τ , which, on its
turn is defined by the constitutive Equation (7.37). Now, it is possible to determine the time
derivate of the internal force vector, Equation (7.32), and similar to Equation (7.26), the
stiffness matrix for geometrically nonlinear behavior can be determined. (The derivation of
the stiffness matrix is beyond the scope of this course).

7.6 Stress update

As described in section 7.4 during a global Newton-Raphson iteration (Finite Element level)
two quantities must be determined. In order to check for convergence, i.e. Fint = Fext , the
internal force must be updated for a given displacement increment and the tangent stiffness
matrix must be computed in order to determine the next guess for the displacement increment.
The plasticity constitutive relations that were presented in the previous chapters were formu-
lated in rate form. In order to use these relations the time derivatives must be integrated using
a suitable numerical method for solving ODEs which usually involves discretizing these equa-
tions also in the time domain. This is also in line with the Newton-Raphson iterations that
92 Finite Element Method

were introduced for the nonlinear finite element method where the concept of (pseudo-)time
was already introduced.
Therefore during an iteration of the Newton-Raphson scheme a trial solution for the displace-
ment increment is found. Based on the displacement increments strain increments at each
integration point in the Finite Element model is found using the strain-displacement matrix
of each element. The reason that only the values at the integration points are required is
due to the fact that the internal forces and the tangent stiffness matrix are computed by
taking the integrals given in equations (7.22) and (7.26), respectively, for which numerical
integration is used.
Within a time increment therefore the stress increment that corresponds to the given strain
increment must be computed that also satisfies all the constitutive relations that are for-
mulated for the chosen material. This procedure is called ‘stress update’. After computing
the stresses at the integration points directly the internal force vector can be computed and
convergence can be checked. If it is not attained then a new trial solution for displacement
must be found using the tangent stiffness matrix. To be able to do that however the relation
between the stress and strain increment is also required as shown in equations (7.26) and
(7.24).
In order to formulate a stress update algorithm first the discrete form of the constitutive
relations are presented. For elastoplastic material behavior, the elastic strain tensor is related
to the Cauchy stress tensor via the elasticity tensor E. When the total strain is decomposed
into a plastic part and an elastic part, the elastic stress-strain relation in incremental form
yields:
∆σ = E : ∆εe = E : (∆ε − ∆εp ) (7.43)
The plastic strain rate is integrated by a generalized trapezoidal rule, yielding the following
expression for the incremental plastic strain:
t(n+1)
Z
∆εp = ε̇p dt ≈ (1 − α)ε̇p(n) ∆t + αε̇p(n+1) ∆t (7.44)
t(n)

where t(n) and t(n+1) bound a time increment. The direction of the incremental plastic strain
is determined by interpolating the plastic strain rate directions at the start and end of the
time increment with the interpolation value α. Several well known integration techniques are
related to a fixed value of the interpolation value α. Taking α = 0 will lead to the Euler
forward method which is based on values at the start of the step only. This algorithm is
an explicit algorithm since the strain or stress of the current step has no influence on the
integration algorithm. Taking α = 1 will lead to the Euler backward method which is based
only on values at the end of the step and thus is a fully implicit method. In the remainder of
this section, the focus will be on the Euler Backward method since this integration scheme is
unconditionally stable, contrary to the Euler forward method, see Figure 7.4.

An expression for the plastic strain increment is obtained by substituting the definition of the
∂φ
flow rule (i.e. ε̇p = λ̇ ) into Equation (7.44) and taking α = 1:
∂σ
∂φ ∂φ
∆εp = ε̇p(n+1) ∆t = λ̇(n+1) ∆t ≈ ∆λ(n+1) (7.45)
∂σ(n+1) ∂σ (n+1)
7.6 Stress update 93

Then, substitution of Equation (7.45) in (7.43) yields the expression for the current stress
state:

σ (n+1) = σ(n) + ∆σ
∂φ
= σ(n) + E : ∆ε − ∆λ(n+1) E : (7.46)
∂σ (n+1)

The new stress state σ (n+1) is determined using an elastic predictor / plastic corrector method.
The basic idea is to treat the total strain increment as elastic, calculate the corresponding
stress, and then project this calculated stress onto the yield surface. The elastic predictor
defines a trial stress state σ t :
σt = σ (n) + E : ∆ε (7.47)

σ (n) is the stress state at time t(n) , see Figure 7.4. If the trial stress lies within the yield
surface, then φ(σ t , εpeq(n) ) ≤ 0, and the trial stress will be the new stress state: σ (n+1) = σ t .
If the trial stress state lies outside the yield surface then φ(σ t , εpeq(n) ) > 0. Since the yield
function can never be larger than 0, a plastic corrector has to be used to determine a new
stress state σ (n+1) which lies on the yield surface φ(σ (n+1) , εpeq(n+1) ) = 0. In case of hardening,
the yield surface changed due to plastic deformation, see dashed line in Figure 7.4.
Note: This figure also shows why the Euler forward method is conditionally stable. The

Figure 7.4: Return mapping

Euler forward method is based on values at the start of the increment, so the direction in
which the trial stress is scaled to the yield surface equals ∂σ∂φ . Therefore it is possible that
(n)
the return mapping algorithm cannot find a solution for the new stress state σ̃ (n+1) , using
the Euler forward stress update.
The equations presented above can be solved using the conventional backward-Euler approach
by defining a set of non-linear equations to be satisfied at time (n+1). However this requires a
non-linear solution scheme such as the Newton-Raphson method for which second derivatives
of the yield function with respect to stress and strain are required. However other algorithms
also exist that circumvent this thereby reducing the complexity of implementation. In the
following one of these algorithms is presented.
94 Finite Element Method

7.6.1 Cutting-plane algorithm

This method can be applied in general to different types of yield functions and hardening rules.
The new stress state σ (n+1) can be determined by writing Equation (7.46) as a function of
the trial stress state:
∂φ
σ(n+1) = σt − ∆λ(n+1) E : (7.48)
∂σ (n+1)
Besides, the new stress state must be on the yield surface:

φ(σ (n+1) , εpeq(n+1) ) = 0 (7.49)

It is clear that the plastic multiplier ∆λ must be known to determine the new stress state.
The value of ∆λ can be determined by a Taylor expansion of the yield function around the
trial stress state σ t :

∂φ ∂φ
φ(σ (n+1) , εpeq(n+1) ) = φ(σ t , εpeq(n) ) + : ∆σ + p ∆εp + O(∆2 ) = 0 (7.50)
∂σ t ∂εeq(n) eq(n)

Neglecting the higher order terms (O(∆2 )) and taking into account that ∆σ = σ (n+1) − σt
gives an expression for the plastic multiplier ∆λ:

∂φ ∂φ ∂φ
φ(σ t , εpeq(n) ) − ∆λ :E: + p ∆εpeq(n) = 0 (7.51)
∂σ t ∂σ t ∂εeq(n)

The equivalent plastic strain can be written as a function of the plastic multiplier. Starting
point is that the incremental plastic deformation energy can be defined by the product of
the incremental strain tensor and the stress tensor, as well as by the incremental equivalent
plastic strain multiplied by the equivalent stress:

∂φ ∆λ ∂φ
σy ∆εpeq = σ : ∆εp = ∆λσ : ⇔ ∆εpeq = σ: (7.52)
∂σ σy ∂σ

Substitution of Equation (7.52) into (7.51) gives an expression for the plastic multiplier. The
plastic multiplier can be calculated and subsequently, the new stress state σ (n+1) is obtained
by Equation (7.48).

Due to the inaccuracy introduced by the linearization given in (7.50) however it is expected
that the yield function will be violated after one iteration. However, replacing the trial stress
with the newly found stress and repeating the process will give an iterative approach in
searching for the solution without requiring second derivatives. Below an example is given
for the cutting-plane algorithm applied to Von Mises plasticity.

Von Mises plasticity

The Von Mises yield criterion will be applied in root form given by:
q √
φ = (σxx − σyy )2 + (σyy − σzz )2 + (σzz − σxx )2 + 6τxy
2 + 6τ 2 + 6τ 2 − 2σ
xz yz y
√
= σf − 2σy = 0 (7.53)
7.6 Stress update 95

The derivative of φ with respect to σ reads:

 


 2σ xx − σ yy − σ zz 


 
2σ − σ − σ
 
yy xx zz 

 

 

 
∂φ 1 2σzz − σyy − σxx 

= (7.54)
∂σ σf 
 6τxy 


 

6τ
 
xz

 


 

 
 6τyz 

Substitution of Equation (7.54) into Equation (7.52) yields a relation between the plastic
multiplier ∆λ and the equivalent plastic strain increment ∆εpeq :
√
∆εpeq = 2∆λ (7.55)

Substitution of Equation (7.55) into Equation (7.51) and taking into account the chain rule
∂φ ∂σy
for ∂φ
∂ε = ∂σy ∂ε gives:

∂φ ∂φ ∂φ ∂σy √
φ(σ (k−1) , εp(k−1)
eq ) − ∆λ(k) : E : + 2∆λ(k) = 0 (7.56)
∂σ (k−1) ∂σ (k−1) ∂σy(k−1) ∂εp(k−1)
eq

Note that the trial stress is replaced by a stress state with superscript (k−1) . This is because
the above expression will be used in an iterative procedure where (k) denotes the current
iteration. The stress σ(k−1) for k = 1 equals the trial stress σ t . Rewriting the above equation
gives an expression for the plastic multiplier in iteration k:
p(k−1)
φ(σ (k−1) , εeq )
∆λ(k) = (7.57)
∂φ ∂φ ∂σy
:E: + 2 p(k−1)
∂σ (k−1) ∂σ (k−1)
∂εeq

The total plastic multiplier is then:

q
X
∆λ = ∆λ(k) (7.58)
k=1

A new stress state in the kth iteration is determined analogous to equation (7.48):

∂φ
σ (k) = σt − ∆λE : (7.59)
∂σ (k−1)
Note: the direction of the stress projection (direction of the plastic strain increment) depends
on the stress state of the previous iteration (σ (k−1) ). In the limit, the direction of the
plastic strain increment depending on the new stress state (σ (k) ) and the one depending
on the stress state of the previous iteration will be almost equal. This means that the new
stress state depends on the direction of plastic strain increment of the new stress state,
which is characteristic for an Euler Backward method. Also note that the new stress state
is is determined by projecting the stress state from the trial stress, using the total plastic
multiplier. If the new stress state σ (k) fulfills the yield criterion, then the new stress state
96 Finite Element Method

σ (n+1) = σ (k) and the iterative process can be stopped. Summarizing, the iterative stress
update algorithm looks as follows:
initialize
σy = σy(n−1)
dσy
determine
dεpeq
k=0
σ (k) = σ t
determine φ(σ (k) , σy )
 
p(k)
WHILE (kφ σ (k) , σy (εeq ) k > tolerance) DO
k := k + 1
determine ∆λ(k)
∆λ(k)
P
determine ∆λ =
p(k) √
determine ∆εeq = 2∆λ
p(k) p(k)
determine εeq = εpeq(n) + ∆εeq
p(k)
determine σy (εeq )
dσy
determine p(k)
dεeq
∂φ
determine σ (k) = σ t − ∆λE :
∂σ (k−1)
END
Note that the subscript (n) is omitted to enhance readability.

7.7 A note on Finite Element implementation

This section summarizes the incremental procedure and the iterative loops of the finite element
method in case of non-linear material behavior (plasticity). A vector-matrix representation
will be used in this section. Hence, a second order tensor will be represented by a vector (e.g.
T
{σ} = {σxx σyy σzz τxy τyz τxz } ) while a fourth order tensor will be represented by a
matrix. Superscripts denote the current iteration belonging to the iterative Newton–Raphson
procedure, the subscript denotes the time increment (the load history of the problem at hand
is taken into account). Based on the equilibrium Equation (7.16), the following linearized
equation has to be solved (see Equation (7.20):

(k) (k+1) (k)

[K(n) ]{∆∆u(n) } = {R(n) } (7.60)

with the residual force:

(k) (k)
{R(n) } = {F(n) }ext − {F(n) }int (7.61)

and the stiffness matrix (geometric linearized case):

Z
(k) (k) T y(k) (k)
[K(n) ] = [B(n) ] [L(n) ][B(n) ] dV (7.62)
7.7 A note on Finite Element implementation 97

Note that the initial residual force vector is the external force vector itself, because the
internal force vector is {0} in the first iteration (k = 1). Additionally, the deviation between
the internal force vector and external force vector of time step (n) can be accounted for in
time step (n+1), by introducing this difference as an additional external force. This process is
y(k)
called load correction. The material matrix [L(n) ] represents elastic or elastic–plastic material
behavior. The B-tensor is here formulated as a matrix, contrary to the formulation given in
Equation (7.14) where it is defined as a 3rd -order tensor. In the intermezzo below, it is shown
how the [B]-matrix can be constructed.

As an illustration, the [B]-matrix will be derived for a bi-linear rectangular element The
element has a length of 2a in x-direction and a length of 2b in y-direction. The origin of
the coordinate system is located in the middle of the element. In matrix-vector notation,
Equation (7.10) can be written as:

{v} = [N ]{v̄}
 


 vx1 


 
vy1 

 
( ) " #
 

vx N1 0 N2 0 N3 0 N4 0  
= vx2 (7.63)
vy 0 N1 0 N2 0 N3 0 N4 
 .. 

. 

 

 

 

vy4 
 

with {v̄} the nodal velocities and:

(a − x)(b − y) (a + x)(b − y)
N1 = N2 =
4ab 4ab
(a + x)(b + y) (a − x)(b + y)
N3 = N4 = (7.64)
4ab 4ab
Note that the vector with nodal velocities and the interpolation functions are defined differ-
ently compared to Section 7.3.
The rate of deformation vector reads:
   ∂ 
0
 Dx   ∂x

  ( )
 v
∂  x
Dy =  0 (7.65)

∂y 
   ∂ ∂  vy

D   1 1
xy 2 ∂y 2 ∂x

Substitution of Equation (7.63) in (7.65) gives the expression for the [B]-matrix ({D} =
[B]{v̄}):
 
−(b − y) 0 (b − y) 0 (b + y) 0 −(b + y) 0
1  
[B] =  0 −(a − x) 0 −(a + x) 0 (a + x) 0 (a − x) 
4ab  
−(a − x) −(b − y) −(a + x) (b − y) (a + x) (b + y) (a − x) −(b + y)
(7.66)
(k+1)
Solving Equation (7.60) gives an expression for {∆∆u(n) }.
98 Finite Element Method

Next, the displacements and the coordinates are updated:

(k+1) (k) (k+1)
{∆u(n) } = {∆u(n) } + {∆∆u(n) } (7.67)

(k+1) (k+1)
{x(n) } = {x(n−1) } + {∆u(n) } (7.68)
(k)
The incremental displacement {∆u(n) } is {0} in the first iteration (k = 1).
Then, the incremental strain is obtained via:
(k+1) (k+1) (k+1) (k+1) (k+1) (k+1)
{∆ε(n) } = {D(n) }∆t = [B(n) ]{v(n) }∆t = [B(n) ]{∆u(n) } (7.69)

and subsequently, the accompanying stress is calculated by a stress update algorithm (see
Section 7.6):  
(k+1) (k+1)
{σ(n) } = f {∆ε(n) }, {h} (7.70)

with {h} a vector with history variables. Hence, the internal reaction force is obtained
through: Z
(k+1) (k+1) T (k+1)
{F(n) }int = [B(n) ] {σ(n) } dV (7.71)
V
(k+1) (k+1)
If the norm of the residual force vector {R(n) } = {F(n) }ext −{F(n) }int is larger than a user
specified tolerance, the Newton–Raphson iterative procedure will be continued: k := k + 1.
The Newton–Raphson loop is continued until the norm of the residual force vector is smaller
than the user specified tolerance. Subsequently, the data has to be prepared for the next
time increment, i.e. the results at the end of increment (n) will now be used as the starting
conditions for increment (n + 1).
Chapter 8

Continuum Damage Mechanics

8.1 Introduction
In previous chapters, we considered elastic and (visco)plastic deformations. According to the
presented elastic or plastic material models, infinitely large deformations can be imposed.
A real material will fail earlier by fracture. Before final fracture the material locally looses
cohesion and microscopic voids can be detected. Upon further loading these voids grow in
size and neighbouring voids coalesce and form micro-cracks. These voids and micro-cracks
are also known as damage.
In Figure 8.1 a micrograph is presented of a cross section through the necked region of a tensile
specimen from a dual phase steel. The black dots are voids. You can see that there is not

Figure 8.1: Voids (black) in a DP600 steel, close to a tensile fracture on the right-hand side
of the picture.

99
 100 Continuum Damage Mechanics

one void where the final fracture initiates, but the voids are distributed over a large region.
The models that are discussed in this chapter consider these distributed voids as part of
the continuum and are therefore classified as Continuum damage mechanics (CDM) models.
CDM models distinguish commonly between 3 stages: initiation, growth and coalescence of
voids. After coalescence, micro-cracks exist that grow into macroscopic cracks that lead to
fracture. Macroscopic cracks are typical discrete discontinuities in a material that cannot be
modelled with continuum models.

8.2 Uncoupled damage models

Uncoupled damage models are also called fracture models or fracture criteria. The evolution
of damage is monitored by updating a damage variable as function of stress and deformation,
but the value of this variable has no direct effect on the local stress–strain relation. A critical
damage value is included in these models that corresponds to the full desintegration of the
material. The damage variable is used as a post-processing value to show the damaged state
of the structure, or it is used to delete elements from a finite element simulation if the critical
value is reached. In the latter case, the evolution of damage has influence on the calculation
results, because the load that was previously transferred through the element must now be
redistributed to other parts of the structure.
It is typically observed that the amount of void growth depends on the triaxiality of the stress
state. The triaxiality η is defined as
p
η= (8.1)
σeq
where p = σii /3 is the hydrostatic stress component and σeq is the (Von Mises) equivalent
stress. Rice and Tracey fitted fracture strains at different triaxialities and formulated the
fracture strain εpf as
εpf = 2.48 exp(−1.5η) (8.2)
This formulation can only be used for deformation at constant triaxiality. For a general
deformation path a damage variable D is defined as an integral over the path:
Z ε̄
D(ε̄) = f (σ, . . .) dεp (8.3)
0

such that a critical value Dc is reached when the equivalent plastic strain reaches the fracture
strain:
D(εpf ) = Dc (8.4)
Some famous fracture criteria use the following functions in (8.3):

(based on highest principal stress) Cockroft and Latham (1968) f (σ) = σ1 (8.5)
Rice and Tracey (1969) f (σ) = exp(Cη) (8.6)
Oyane (1980) f (σ) = 1 + Cη (8.7)
p B
Goijaerts (1999) f (σ) = h1 + Aηi (ε ) (8.8)

In the Cockroft and Latham criterion, the highest principal stress is used and integrated over
the equivalent strain. The Rice and Tracey model directly use the triaxiality and gives it
an exponential weight in the integration. The Oyane model assumes a linear influence of
8.3 Coupled damage models 101

the triaxiality. Goijaerts adapted this model by using the Macaulay brackets (hxi equals x
if x > 0 and equals 0 if x ≤ 0). This assumes that no damage evolution takes place under
negative triaxiality. Furthermore, a better fit could be obtained by weighing with a power
function of the equivalent plastic strain.
Because the uncoupled damage models are not directly taking the damage into account while
plastic deformation takes place, these models are especially suited for ductile materials, where
significant plastic deformation takes place before damage evolves.

8.3 Coupled damage models

8.3.1 Basic concepts
The basic concept of continuum damage mechanics is to describe the relation between nominal
stresses σ and strains ε as function of the effective stresses σ eff and strains εeff in between
the voids and a damage variable D. The basic idea for describing damage in this way is
due to Kachanov in 1958 [5]. We will restrict ourselves in this course to isotropic damage
formulations, in which the damage variable D is a scalar.
The simplest form of damage mechanics models considers the linear elastic deformation of
material, but assuming that the externally observed (nominal) stiffness, depends linearly on
the effective stiffness of the (undamaged) material and the damage variable D. In terms of
the Young’s modulus:
E = (1 − D)Eeff (8.9)

Using this model, the material will have the full, undamaged, stiffness if the damage variable
has the value zero and it will have no stiffness left if the damage variable has the value one.
It then has no remaining load bearing capacity.
It can be expected that voids in a material will affect both the stress and strain, compared
to an unvoided material. For modelling purposes, simplifying assumptions are made. Two
different strategies are known as the stress equivalence principle and the strain equivalence
principle. According to the stress equivalence principle, the nominal stress and the effective
stress are equal and according to the strain equivalence principle, the nominal strain and the
effective strain are equal. The strain equivalence principle is mostly used and also adopted
here, hence:
ε = εeff (This is a second order tensor) (8.10)

In the undamaged material, we assume the Hookean relation σ eff = Eeff : εeff , where Eeff is
proportional with the undamaged Young’s modulus Eeff . For a damaged material we then
get
σ = E : ε = (1 − D)Eeff : εeff = (1 − D)σ eff (8.11)

An alternative interpretation is that the nominal stress and effective stress are related by the
effective load bearing area. In a uniaxial tensile test, the effective stress is the load divided
by the effective load bearing area, while the nominal stress is the load divided by the nominal
area. The damage variable can be considered as the void volume fraction.
 
(Void volume fraction)
Vd
Veff = V − Vd = 1− V = (1 − D)V (8.12)
V
102 Continuum Damage Mechanics

σf

σ Eeff

ε0 εmax
ε

Figure 8.2: Brittle fracture modelled by a linear damage evolution.

where Vd is the void volume and V is the nominal volume. For spherical voids the void volume
fraction equals the void area fraction in a cross section:

(Void area fraction) Aeff = (1 − D)A (8.13)

such that
Aeff
F = σA = (1 − D)σeff = σeff Aeff (8.14)
1−D

8.3.2 Brittle fracture

A brittle fracture model can be derived from purely elastic behaviour and damage evolution.
The damage variable D has an initial value of 0 and can only grow, up to a value of 1. In a
uniaxial tensile test, it could be postulated that the material behaves linear elastic up to a
maximum value of the stress σf , after which the damage grows linearly with the strain, up to
a maximum strain of εmax where the damage variable reaches the value of 1. The strain at
which damage starts to evolve is ε0 = σf /Eeff .
For a monotonic test, the damage variable can be defined as:

0
 if ε < ε0
D= 1 if ε > εmax (8.15)
 ε−ε0

εmax −ε0 otherwise

Damage, however, can only grow and if the material is unloaded, the damage remains constant.
Therefore a parameter κ is introduced that represents the maximum obtained non-elastic
strain.
κ(t) = max(ε(τ )) 0 ≤ τ ≤ t (8.16)
In Equation (8.15), ε is then substituted by κ
The nominal Young’s modulus is derived from Equation (8.9). If the material is unloaded
from a damaged state, the stress–strain relation follows the line of reduced stiffness E, down
to the origin. This behaviour is depicted in Figure 8.2.
8.3 Coupled damage models 103

As long as the stress has not reached the fracture stress σf , the material is undamaged. When
the fracture stress is reached, micro-cracks initiate and evolve, leading to a reduced load
bearing capacity. Only when the strain reaches a value of εmax the load bearing capacity is
completely lost.
In the decreasing part of the stress–strain curve, the stress is derived from
 
κ − ε0
σ = (1 − D)Eeff κ = 1 − Eeff κ (8.17)
εmax − ε0

It should be noted that a negative slope in the stress–strain relation will result in mesh size
dependent results in finite element simulations and it requires a regularization or the use of
mesh size dependent ‘material’ parameters to get physically relevant results.
The approach for a uniaxial brittle damage model can easily be extended to the general case
by making the damage evolution dependent on an equivalent strain measure. It is however
questionable if damage due to micro-cracks should be modelled with isotropic damage or with
a more complex anisotropic damage model. In the latter case the damage variable D should
be written as a second, or even a fourth order tensor. This is beyond the scope of this course.

8.3.3 Ductile damage

In materials that show ductile fracture, plastic deformation dominates before final rupture.
The same approach can be taken as in the previous section, but now the basic material model
is formed by an elastoplastic model, usually with work hardening. The stress is then directly
calculated from
σ = (1 − D)σ eff (8.18)
where σ eff is calculated using the elastoplastic model. For many steels, it is observed that voids
are created at, or just before the maximum force in a tensile test. Damage models therefore
usually take this moment for initiation of damage. Including damage in the material model
avoids that a too large post-necking strain is predicted. With damage, the stress quickly
drops after the load maximum, like also presented in the brittle fracture model above.
A popular approach for modelling of ductile damage was proposed by Gurson. Originally, he
investigated void growth in an ideally plastic material. The void growth is then incorporated
in a yield function that presents the same behaviour as a matrix material modelled with
the Von Mises model with embedded voids. The presence of voids makes the model com-
pressible and dependent on the hydrostatic stress p, unlike the classical plasticity models for
metals. Modifications by Tvergaard and Needleman, resulted in the so-called GTN (Gurson,
Tvergaard, Needleman) model. The yield function can be written as:
2
σeq
 
q2 p
φ= + 2q1 f ∗ cosh − (1 + (q1 f ∗ )2 ) = 0 (8.19)
σf2 2σf

with f ∗ the effective void volume fraction

(
∗ f if f ≤ fc
f = 1/q1 −fc (8.20)
fc + ff −fc if f > fc

σf is the current flow stress, fc and ff are the void volume fractions at start of coalescence
of voids and at final rupture respectively. The model parameters q1 and q2 were added for a
104 Continuum Damage Mechanics

better fit. Setting them equal to one and using f = f ∗ degenerates the GTN model to the
original Gurson model.
The void volume evolution follows from the change of volume of the material, because the
plastic deformation of the matrix material is considered incompressible. This results in
 
˙ 2 3p
f = λ̇(f − f )σf sinh (8.21)
2σf

This term describes the growth of existing voids. Often an additional term for void initiation
is used.
Contrary to the original damage model by Kachanov, in the case of ductile damage, the
Young’s modulus is often taken constant and is not affected by damage.
Chapter 9

Tensile instability in sheet metal

In this chapter the combined effect of material hardening and geometrical changes on the
stability of deformation is considered for planar stretching of sheet metal. The stability
analysis resembles the analysis for a tensile bar, but it is now extended to stretching of sheet
under bi-axial stress states. The analysis is relevant for determination of the formability limits
in sheet metal.
The limiting values for the deformation of sheet metal depend on many parameters. One
important type of failure is the fracture of sheet during the forming process. This fracture is
preceded by localized necking and depends on the current state of deformation. In the plane
of the sheet, the deformation in two directions influence the failure behavior, e.g. εx = 0.2
and εy = 0.2 is a completely different situation than εx = 0.2 and εy = −0.2. The first is a
state of equi-biaxial stretch and the second a state of pure shear.
In a uniaxial tensile test localized necking is preceded by diffuse necking. After a small
introduction to plasticity, this type of necking is treated first. Thereafter localized necking
for different proportional strain paths is described.

9.1 Plastic deformation

For this chapter only very little theory on plastic deformation is needed. In this section we
define the basic relations that are used in the sequel.
First of all, we neglect elastic deformations. Since we consider relatively large plastic defor-
mations for metals, the elastic deformation can be ignored when we calculate the change in
geometry. Furthermore, for metals, the plastic deformation is related to slip on glide sys-
tems and therefore the total volume remains constant. In this chapter strains are always
logarithmic strains and therefore, the volume constraint requires

εx + εy + εz = 0 (9.1)

In a tensile test, we observe that after initial yielding of material, the tensile stress increases.
This is called work hardening and is often modelled with the Hollomon law:

σy = Cεn (Holloman law) (9.2)

where σy is the flow stress at a particular strain ε. For many steels this is a good approximation
for a large part of the hardening curve. To model the initial yield stress, an initial strain ε0

105
106 Tensile instability in sheet metal

(a) (b) (c)

Figure 9.1: Sheet forming instabilities: (a) uniform deformation, (b) diffuse necking and (c)
local necking.

is used in the Swift law:

σy = C(ε + ε0 )n (Swift law) (9.3)
The Hollomon and Swift relations directly apply to a uniaxial tensile test, but for a plane
stress state an equivalent stress must be defined that can be related to the flow stress σy . A
simple relation in principal stresses is given by the Von Mises criterion:
q
σy = σ12 + σ22 − σ1 σ2 (9.4)

Plastic deformation takes place if the equivalent stress reaches the value of the flow stress.
The strain ε in the Hollomon equation is then substituted by an equivalent strain that is
derived from the principle of energy equivalence:
Z t
εeq (t) = ε̇eq (τ ) dτ and σeq ε̇eq = σij ε̇ij (9.5)
0

In practice the Von Mises yield locus appears to be too rounded. To describe the higher cur-
vature that is often found near the uniaxial and equi-biaxial stress states, Hosford developed
an isotropic non-quadratic model in 1972 [4]. For the plane stress case the model can be
written as:
σeq = 21 (|σ1 − σ2 |m + |σ1 |m + |σ2 |m )1/m (9.6)
Experience shows good results for b.c.c. metals with m = 6 and for f.c.c. metals with m = 8.
The plane stress yield loci for m = 2 (Von Mises) and for m = 6 and m = 8 are presented in
9.2.
9.2 Diffuse necking 107

σ2

σ1

m=2
m=6
m=8

Figure 9.2: The Hosford yield locus for 3 different values of m.

Finally, the direction of plastic strain is perpendicular to the yield surface, or in other words
to the (hyper)plane of constant equivalent stress:
∂σeq
ε̇ij = λ (9.7)
∂σij

with the plastic multiplier λ.

9.2 Diffuse necking

For the analysis of diffuse necking, we consider a tensile specimen as in 9.1 with width w,
thickness t and gauge length l, the constant volume constraint requires that

A0 l0 = Al with A = wt (9.8)

An infinitesimal strain increment can then be written as:

dε1 = dl/l = − dA/A (9.9)

where the 1-direction indicates the tensile direction. The tensile load P on the specimen is

P = σ1 A (9.10)

or in differential form
dP = dσ1 A + σ1 dA (9.11)
At the load maximum we have dP = 0 and we obtain
dσ1 dA
=− = dε1 (9.12)
σ1 A
108 Tensile instability in sheet metal

or
dσ1
= σ1 (9.13)
dε1
Applying the Hollomon hardening function σ1 = Cεn1 the condition for diffuse necking be-
comes:
nCε1n−1 = Cεn1 ⇒ n = ε∗1 (9.14)

The value ε∗1 is called the uniform strain. For many steels the value of n is in the order of 0.2
and the preceding formula indicates that a uniaxial tensile specimen will yield in a uniform
way up to a natural strain of about 0.2. At that point diffuse necking will start as indicated
in 9.1.b. In the load–deflection curve this point coincides with the maximum load.

9.3 Localized necking

In a uniaxial tensile test, a relatively small strip is loaded in tension. Therefore the non-
necking part of the specimen does not have a large influence on the necking part. In normal
sheet forming processes this is not the case. Usually, if one region in the sheet is deforming
that would induce diffuse necking in an unconstrained situation, the rest of the sheets im-
poses constraints that precludes lateral deformation. The deformations in sheet forming must
usually be compatible with a rigid tool translation; therefore diffuse necking does not limit
sheet forming processes as in a tensile test.
In sheet forming, necking is possible if it is localized, so that it does not influence the global
strain distribution. These necks typically have a neck size in the order of the sheet thickness.
In thickness direction the neck is diffuse, but in the sheet plane the neck is localized.
The deformation states of material points can be represented in a so-called forming limit
diagram (FLD). In an FLD the major and minor in-plane strains (the principal strains in
the plane of the sheet) are plotted in a graph as in 9.3. The upper boundary of feasible
deformation points represents the deformation at which localization occurs. This boundary
is called the forming limit curve (FLC).
Experimental forming limit curves are derived from drawn specimens with a circle pattern
imprinted on the sheet. The circles deform to ellipses and the dimensions of the ellipse
along the major and minor axes, related to the original diameter determine the major and
minor principal strain. Formally, the forming limit diagrams are only derived for proportional
loading and non-proportional loading may yield different results.
FLDs are often determined by an experiment, devised by Nakazima et al (see 9.3). A hemi-
spherical punch deforms strips of fixed length but varying width. Strain ratios between −0.5
(uniaxial) and 1.0 (equi-biaxial) can be achieved in this way. FLDs can also be determined
theoretically. A distinction is made between the left part of the FLD (the drawing region)
and the part on the right (the stretching region). Both regions are treated next.

9.3.1 Localization in the drawing region

We first have a look at the localized necking criterion according to Hill [2]. For this purpose
we consider a sheet, biaxially loaded in a proportional process. As long as the deformations
are uniform, we have a constant ratio between the two in-plane principal stresses and another
constant ratio between the strain increments in the principal directions. Taking the principal
9.3 Localized necking 109

Figure 9.3: FLD determination using the strip technique of Nakazima et al [8].

stress and principal strain directions to be equal restricts the validity of the model to planar
isotropic material behavior.
σ2 = ασ1 , σ3 = 0 (9.15)
dε2 = β dε1 , dε3 = −(1 + β) dε1 (9.16)
The unit forces (per sheet length) are defined as

T1 = tσ1 , T2 = tσ2 (9.17)

Local necking is assumed to start if dT = 0. Due to the proportional loading, this will occur
for T1 and T2 simultaneously. We now have

dT1 = dtσ1 + t dσ1 = t(σ1 dε3 + dσ1 ) = t[−σ1 (1 + β) dε1 + dσ1 ] = 0 (9.18)

and the condition for localized necking becomes:

dσ1
= (1 + β)σ1 (9.19)
dε1

If the hardening of the material can be described by the Hollomon relation (9.2) it can be
derived that for a proportional deformation process the relation

σ1 = C ′ εn1 (9.20)
110 Tensile instability in sheet metal

holds, where C ′ = Cf (β, n). In this case the localized necking condition (9.19) becomes
nσ1
nC ′ ε1n−1 = = (1 + β)σ1 (9.21)
ε1

Therefore the strains at maximum T1 are:

ε∗1 = n/(1 + β) ε∗2 = βn/(1 + β) (9.22)

Taking the summation of ε∗1 and ε∗2 we see that

ε∗1 + ε∗2 = n (9.23)

In a forming limit diagram this condition can be represented by a straight line through the
point where the major strain is n and the minor strain is 0, parallel to the line ε1 = −ε2 .
For the Hollomon hardening function, the limit strains are independent of the yield function.
This is not the case for other hardening functions.

✻
σ2

θ t n 2 σ1
✛ ❑✯ ✻
✲1 ✲

Figure 9.4: Localization band in a uniformly strained region

If necking takes place in the described way, a localization band emerges at a specific angle to
the major strain direction (9.4). The straining in the band accelerates, while the strain outside
the band remains constant (neglecting elastic springback). In order to keep the displacements
inside and outside the band compatible, the strain in the band, in the direction along the band
must be zero. Therefore, since we have proportional loading up to the localization strain, the
band will develop in a direction where the strain is zero. The strain εtt in a direction with
angle θ to the major strain direction can be derived from the following equation (graphically
represented by the Mohr circle):

εtt = 21 (ε1 + ε2 ) + 12 (ε1 − ε2 ) cos 2θ (9.24)

The constraint εtt = 0 yields the condition

ε1 + ε2 1+β
cos 2θ = − =− (9.25)
ε1 − ε2 1−β
9.3 Localized necking 111

A
1 dεA

dεB
B 2 P
dεA dεB
Q dε2
A arctan α0
dε1A dε1B
φ=0
(a) Plate with imperfection. (b) Strain increments inside and outside the groove.

Figure 9.5: Principle of the Marciniak–Kuczynski analysis.

It can be concluded from this result that at values β > 0 no solution for θ exists. Therefore,
the localization criterion according to Hill is only valid for the extension–compression region
(the left side) of the FLD.
For uniaxial tensile stresses and isotropic material, the width and thickness strain are equal
and from the constant volume constraint we obtain

dε1 = − dε2 − dε3 = −2 dε2 (9.26)

and therefore
dε2 1
β= =− (9.27)
dε1 2

For uniaxial loading of an isotropic material we will get a localization band at an angle

1+β 1 1
arccos − 31 ≈ 55◦


cos 2θ = − =− ⇒ θ= 2 (9.28)
1−β 3

or 35◦ off from a neck perpendicular to the tensile direction.

9.3.2 Localisation in the stretching region

It was shown in the previous section that under proportional biaxial extension (ε1 > 0 and
ε2 > 0) localised necking within a sheet is impossible. However, due to inhomogeneities,
the local rate of deformation can rotate to a plane strain state, facilitating localised neck-
ing. Marciniak and Kuczynski introduced an analysis used for the determination of localised
necking in biaxial extension [7]. In their analysis a sheet is considered with an initial groove
with reduced thickness, as presented in 9.5(a). Quantities referring to the groove are given
an index B and quantities referring to the rest of the sheet are given an index A. The initial
thickness of the groove is tB,0 and outside the groove tA,0 . In practice, in the stretching
regime a neck is observed perpendicular to the direction of the major strain. Therefore, the
initial groove is aligned with the minor strain in the Marciniak–Kuczynski (M–K) analysis.
112 Tensile instability in sheet metal

Marciniak–Kuczynski Analysis Outside the groove a proportional deformation path is

assumed, described by

σ2A = α0 σ1A σ3A = 0 (9.29)

ε2A = β0 ε1A ε3A = −(1 + β0 )ε1A (9.30)

Compatibility between the uniform part A and the groove B requires that

dε2A = dε2B (9.31)

The force per unit sheet length in direction 1 must be transmitted through the groove, hence

T1 = σ1A tA = σ1B tB ⇒ σ1B = σ1A /f (9.32)

Where f = tB /tA is the current thickness ratio. As long as f ≈ 1, the stress ratio α0
approximately holds for both regions A and B. Since the stress σ1B in the groove is larger
than σ1A in the uniform part, the material in the groove reaches the yield surface first. In
9.5(b) this is approximately at position P . Because of the constraint equation (9.31) no
yielding takes place, since the uniform region is still fully elastic. The stress state in region
B will move along the yield locus to position Q, until also region A reaches the yield locus
at position P . This situation is depicted in 9.5(b). In that period σ1B has increased and σ2B
has decreased, hence the stress ratio α has decreased and proportional deformation is not
possible in the groove.
With the stress state in A and B at two different positions on the yield locus, the normal to
the yield surface is different and because of the constraint (9.31), the strain perpendicular to
the neck must be larger in the neck than in the uniform part. This is presented in the vector
plot in 9.5(b). As a consequence, the hardening in region B will be larger and the two stress
situations can no longer be projected on the same yield locus. With an equal strain increment
in the 2-direction and a larger strain in the 1-direction, the thickness strain in the neck will
be more negative than in the uniform part, decreasing the thickness ratio (f < f0 ).
The analysis of the deformation can be further developed numerically. The drawing region
can be included in the analysis by assuming an inclined groove. For planar isotropic material
the inclination angle can be determined beforehand by (9.25). Strain increments ∆εA are
prescribed on region A, respecting the proportionality ratio β0 . For every increment ∆εA a
strain increment ∆εB is calculated iteratively, such that the plate forces normal and tangential
to the groove in region A and B are in equilibrium. The sheet is considered to have failed if
the strain increment in the groove is larger than a prescribed multiple of the strain increment
in the uniform part. The strain in the uniform region at the time of failure is considered to
be the limiting strain for that particular strain ratio β0 .
In 9.6 the strain paths in the plate (region A) and in the groove (region B) are presented in
a forming limit diagram. An initial thickness ratio f0 = 0.99 was used. The failure criterion
stopped the increments when the strain increment in the groove was more than 10 times
the strain increment in the uniform region. A Von Mises yield locus was used with a Swift
hardening model with C = 446.3 MPa, n = 0.292 and ε0 = 0.0039.
In the drawing region of the FLD, the strains inside the groove in the directions 1 and 2
have the same ratio as in the uniform region. The limit strain follows a line at an angle of
45◦ with respect to the major strain axis, as predicted by the Hill analysis. The limit strain
for the plane strain path is lower than the corresponding n-value of 0.292 because of the
9.3 Localized necking 113

0.7

0.6

0.5
major strain

0.4

0.3

0.2

0.1
plate
groove
FLC
0
-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6
minor strain

Figure 9.6: The strain-paths inside and outside the groove in a Marciniak–Kuczynski analysis.

Figure 9.7: The Hill theory and the Marciniak–Kuczynski theory against experiment (◦) [8].

initial thickness ratio of 0.99. In the stretching region, the strains inside and outside the
groove are no longer in the same direction. It can be seen that the strain ratio ε1B /ε1A keeps
increasing until the rate of deformation in the groove approaches a plane strain situation. For
a fixed strain increment in part A, the strain increment in the groove becomes very large and
localisation is detected.
The results of the Marciniak–Kuczynski analysis are very sensitive to the initial thickness
ratio f0 (see 9.7). If experimental FLC values are known for the plane strain condition, this
is often used to calibrate f0 . The estimation at the equi-biaxial point is shown to be very
sensitive to the used yield function [9]. If a Von Mises yield function is used, the FLCs tend to
be much too high near the equi-biaxial direction (see 9.7). It can be appreciated from 9.5(b)
114 Tensile instability in sheet metal

0.7

0.6

0.5
major strain

0.4

0.3

0.2

m=2
0.1
m=6
m=8
0
-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
minor strain

Figure 9.8: Predicted FLCs using the Hosford yield function and f0 = 0.995.

that the local curvature of the yield locus is very important. Indeed, if the Hosford yield
function is used for representative values of m the shape of the predicted FLC becomes much
better, as can be seen in 9.8. Notice that—at least for the Swift hardening model—the left-
hand side of the FLC is not influenced by the shape of the yield function. This corresponds
to the results of the Hill local necking analysis.

-- END OF THEORY --
Appendix A

Principal stresses and invariants

In each point, at least three mutual perpendicular planes can be defined on which only a
normal stress acts. Then, the stress vector T has the direction of the normal vector:

T=σ·n (A.1)

The size of the principal stress is σ. The components of the stress vector can be written as
(see equation 4.3):

Tx = σ cos (α) = σxx cos (α) + σxy cos (β) + σxz cos (γ)
Ty = σ cos (β) = σxy cos (α) + σyy cos (β) + σyz cos (γ) (A.2)
Tz = σ cos (γ) = σxz cos (α) + σyz cos (β) + σzz cos (γ)

or in eigen-value representation:
     
σxx − σ σxy σxz 
 cos (α) 

 
0
 
   
 σxy
 σ yy − σ σ yz
 ·
  cos (β) = 0 (A.3)
 
   

σxz σyz σzz − σ  cos (γ)  0

⇐⇒ (σij − σδij )nj = 0 for i 6= j

(
0 for i 6= j
with δij =
1 for i = j
For the non-trivial solution, this eigen-value problem gives the following equation:
 
σxx − σ σxy σxz
 
det 
 σxy σyy − σ σyz  =0 (A.4)
σxz σyz σzz − σ

⇐⇒ |σij − σδij | = 0
Working out this equation yields a secular equation of the third degree, with three roots. The
three roots represent the principal stresses. By definition it is stated that: σ1 ≥ σ2 ≥ σ3 .
Elaboration of the determinant gives:

σ 3 − I1 σ 2 − I2 σ − I3 = 0 (A.5)

115
116 Principal stresses and invariants

The coefficients of this equation are the invariants of the stress tensor:

I1 = tr (σ) = σxx + σyy + σzz = σ1 + σ2 + σ3

I2 = 12 σ : σ − (tr (σ))2 ) = − (σxx σyy + σyy σzz + σzz σxx ) + (σxy
2 2 2


+ σyz + σzx )
(A.6)
= − (σ1 σ2 + σ2 σ3 + σ3 σ1 )
I3 = det(σ) = σ1 σ2 σ3

The sign ”:” means double contraction, comparable to the inner product of 2 vectors. In
index notation, A : B is written as Aij Bij

Deviatoric stress

The first invariant, the sum of the the principal stresses, equals the isotropic part of the
stress tensor:
1 1 1
3 I1 = 3 tr (σ) = 3 (σ1 + σ2 + σ3 ) = −p (A.7)
with p the hydrostatic pressure (average pressure). Next, σ can be split into an isotropic part
and a part of which the average normal stress is zero:

σ = s + 13 tr (σ)1 = s − p1 ⇔ σij = sij − pδij (A.8)

s is the deviatoric stress tensor, which can also be expressed as:

s = σ − 13 tr (σ)1 = 4 I − 13 11 : σ


(A.9)

with 4 I the fourth order identity tensor, defined by:

4 4 4
I : A = A (∀A) ⇔ Iijkl Akl = Aij ⇒ Iijkl = δik δjl (A.10)

The invariants of the deviatoric stress tensor read:

′
I1 = tr (s) = tr (σ) − 13 tr (σ)tr (1) = 0
′
I2 = 21 s : s = 12 (s1 s2 + s2 s3 + s3 s1 ) = I2 − 13 I12
′
2 3
I3 = det(s) = s1 s2 s3 = I3 − 13 I1 I2 + 27 I1

The secular equation for s now becomes:

′ ′
s 3 − I2 s − I3 = 0 (A.11)

Graphical representation, Mohr’s stress circles

Three equations form the foundation of the graphical representation of Mohr’s circles. If
the coordinate system is placed on the principal axis, then equation (A.2) degenerates to:

sx = σxx cos (α) = σ1 cos (α)

sy = σyy cos (β) = σ2 cos (β) (A.12)
sz = σzz cos (γ) = σ3 cos (γ)
 117

With σn = |s · n| = |sx cos (α) + sy cos (β) + sz cos (γ)| one can write:

σn = σ1 cos2 (α) + σ2 cos2 (β) + σ3 cos2 (γ) (A.13)

an expression for the shear stress components can be found by using the theorem of Pythagoras
and the symmetry condition of the stress tensor σij = σji :

τ 2 = |τ |2 = |s|2 − σn2
(A.14)
τ 2 = σ12 cos2 (α) + σ22 cos2 (β) + σ32 cos2 (γ) − σn2

The third equation defines the length of the unity normal:

cos2 (α) + cos2 (β) + cos2 (γ) = 1 (A.15)

From this system of equations, 2 unknowns can be eliminated: 1. Elimination of cos(β) and
cos(γ) gives:
(σn − A1 )2 + τ 2 = R12 (A.16)
This equation represents a circle with:

A1 = 12 (σ1 + σ2 ) (center)
q
R1 = (σ1 − σ2 )(σ1 − σ3 ) cos2 (α) + 41 (σ1 − σ2 )2 (radius, min. for cos(α) = 0)

2. Elimination of cos(α) and cos(β) gives:

(σn − A2 )2 + τ 2 = R22 (A.17)

with:

A2 = 21 (σ1 + σ2 ) (center)
q
R2 = (σ3 − σ1 )(σ3 − σ2 ) cos2 (γ) + 41 (σ1 − σ2 )2 (radius, min. for cos(γ) = 0)

3. Elimination of cos(α) and cos(γ) gives:

(σn − A3 )2 + τ 2 = R32 (A.18)

with:

A3 = 21 (σ1 + σ3 ) (center)
q
R3 = (σ2 − σ1 )(σ2 − σ3 ) cos2 (β) + 14 (σ1 − σ3 )2 (radius, min. for cos(β) = 0)

Three extremes are found for cos(α), cos(β) and cos(γ) = 0, respectively, i.e. Mohr’s stress
circles. The graphical representation of Mohr’s stress circles is depicted in figure A.1. The
stresses acting on an arbitrary plane are somewhere in the hatched areas.

Principal strains and invariants

118 Principal stresses and invariants

Figure A.1: Mohr’s stress circles

The determination of the principal strain is analogue to the determination of the principal
stresses:  
εxx − ε εxy εxz
 
det 
 εxy εyy − ε εyz  =0 (A.19)
εxz εyz εzz − ε
Elaboration of these equations ends up with a secular equation with three roots: i.e. the
principal strains ε1 ≥ ε2 ≥ ε3 . By definition it is stated that: ε1 ≥ ε2 ≥ ε3 . Elaboration of
the determinant gives:
ε3 − J1 ε2 − J2 ε − J3 = 0 (A.20)
The invariants of the strain tensor read:
J1 = tr (ε) = εxx + εyy + εzz = ε1 + ε2 + ε3
J2 = 12 (ε : ε − (tr (ε))2 ) = −εxx εyy + εyy εzz + εzz εxx + (ε2xy + ε2yz + ε2zx )
(A.21)
= −ε1 ε2 + ε2 ε3 + ε3 ε1
J3 = det(ε) = ε1 ε2 ε3

Deviatoric strain

The first invariant, the sum of the principal strains, equals the isotropic part of the strain
tensor. This invariant is related to the specific volume change for small strains:
∆V V − V0
= = (1 + εxx )(1 + εyy )(1 + εzz ) − 1 ≈ εxx + εyy + εzz = J1 (A.22)
V0 V0
Next, ε can be split into a part which involves volume changes (the isotropic part) and a part
that involves shape changes. The isotropic part of the latter is zero.

ε = e + 31 tr (ε)1 (A.23)

e is known as the deviatoric strain tensor and has the following invariants:
′
J1 = tr (e) = exx + eyy + ezz = e1 + e2 + e3
′
J2 = 21 (e : e − (tr (e))2 ) = −exx eyy + eyy ezz + ezz exx + (e2xy + e2yz + e2zx )
(A.24)
= −e1 e2 + e2 e3 + e3 e1
′
J3 = det(e) = e1 e2 e3
 119

Graphical representation, Mohr’s strain circles

The determination of Mohr’s strain circles is analogue to the determination of Mohr’s stress
circles, and the result is depicted in figure A.2, where γ represents the shear strain:

Figure A.2: Mohr’s strain circles

Bibliography

[1] Bonet, J., and Wood, R. D. Nonlinear continuum mechanics for finite element
analysis, 2nd ed. Cambridge University Press, Cambridge, 2008.

[2] Hill, R. On discontinuous plastic states, with special reference to localized necking in
thin sheets. Journal of the Mechanics and Physics of Solids 1 (1952), 19–30.

[3] Holzapfel, G. A. Nonlinear solid mechanics: a continuum approach for engineering.

Wiley, Chichester, 2000.

[4] Hosford, W. F. A generalized isotropic yield criterion. Journal of Applied Mechanics

39 (1972), 607–609.

[5] Kachanov, L. M. On the creep fracture time. (in russian). Izv. Acad. Nauk SSSR Otd.
Tekhn 8 (1958), 26–31.

[6] Krabiell, A., and Dahl, W. Zum Einfluss von Temperatur und Dehngeschwindigkeit
auf die Streckgrenze von Baustählen unterschiedlicher Festigkeit. Archiv für das
Eisenhüttenwesen 52 (1981), 429–436.

[7] Marciniak, Z., and Kuczynski, K. Limit strains in processing of stretch forming
sheet metal. Int. J. Mech. Sci. 9 (1967), 609–620.

[8] Pearce, R. Sheet Metal Forming. IOP Publishing Ltd, 1991.

[9] Vegter, H., An, Y., Pijlman, H. H., and Huétink, J. Different approaches to
describe the plastic material behaviour of steel and aluminium-alloys in sheet forming. In
Proceedings of the 2nd ESAFORM Conference on Material Forming (Guimarães, 1999),
J. A. Covas, Ed., pp. 127–132.

[10] Zener, C., and Hollomon, J. H. Effect of strain rate upon plastic flow of steel.
Journal of Applied Physics 15 (1944), 22.

121