Learn Modeling of Elastoplastic Materials

with Altair Compose

0
Table of Contents
1. Introduction to Modelling of Materials............................................................ 5

2. One Dimensional Models ................................................................................. 7

2.1. Hyperelastic Material ................................................................................... 7

2.1.1. Exercise 1: Modelling the Linear Elastic Material ............................. 10

2.2. Ideal Elastoplastic Material ........................................................................ 16

2.2.1. Algorithm for Ideal Elastoplasticity................................................... 17

2.2.2. Exercise 2: Ideal Elastoplasticity ....................................................... 21

2.3. Elastoplastic Material with Linear Isotropic Hardening ............................. 30

2.3.1. The Algorithm ................................................................................... 33

2.3.2. Exercise 3: Elastoplasticity with Isotropic Hardening ....................... 35

2.4. Elastoplastic Material with Linear Kinematic Hardening ........................... 40

2.5. Elastoplastic Material with Linear Isotropic and Kinematic Hardening ..... 41

2.5.1. The Algorithm ................................................................................... 43

2.5.2. Exercise 4: Elastoplasticity with Linear Isotropic and Kinematic

Hardening ......................................................................................... 45

2.6. Nonlinear Hardening – Ramberg-Osgood Model....................................... 51

2.6.1. The Algorithm ................................................................................... 52

2.6.2. Exercise 5: Elastoplasticity with Ramberg-Osgood Isotropic

Hardening ......................................................................................... 55

3. Three Dimensional Models ............................................................................ 63

3.1. Elastoplasticity with Linear Isotropic and Kinematic Hardening................ 65

3.1.1. Von Mises Yield Criterion.................................................................. 66

3.1.2. The Algorithm ................................................................................... 69

3.1.3. Exercise 6: 3D Elastoplasticity with Linear Isotropic Hardening ....... 72

1
3.1.4. Exercise 7: 3D Elastoplasticity with Linear Isotropic and Kinematic
Hardening ......................................................................................... 91

4. Finite Element Method Application ................................................................ 98

4.1. A Brief Overview of FEM Theory ................................................................ 98

4.1.1. Exercise 8: Elastoplasticity in 1D Truss System ............................... 102

5. Literature 116

2
About this Book
Altair Compose software is an environment for doing calculations, manipulating and
visualizing data (including from CAE simulations or test results), programming and
debugging scripts useful for repeated computations and process automation. This
guide has been prepared for beginners and experienced users to compute the stress
state of a material when given strain history data. By means of computation we will
try to produce accurate stress vs strain figures, that will sufficiently match the
generally observed characteristics during experiments. The first part of the guide will
present 1D material modelling with basics, but also a more advanced Ramberg-
Osgood hardening model. Next, we will create basic material models for 3D, where
tensor calculus will be applied. The capstone exercise of this book will be a simple
nonlinear quasi-static FEM example, where we will investigate elastoplasticity in 1D
two-truss element system.

We assume you have some knowledge of computer programming and understand

concepts like variables, constants, expressions, statements, etc. You can also refer the
book on compose on Altair University called Math, Scripting, Data Analysis &
Visualization with Altair ComposePlease note that a commercially released software
is a living “thing” and so at every release (major or point release) new methods, new
functions are added along with improvement to existing methods. This document is
written using Altair Compose 2017.3 version.

3
Acknowledgment
“If everyone is moving forward together, then success takes care
of itself”
Henry Ford (1863 -1947)
A very special Thank You goes to:

Gabriel Stankiewicz (Altair University Ambassador, email: [email protected]

),for creating all the chapters in this book,

Manuel Ramsaier (Altair University Ambassador, email:

[email protected] )

Rahul Ponginan for thoroughly reviewing the book, formatting and editing,

Priyanka Nagaraj for editing, formatting,

Markus Haller and Stefan Klug (Germany), Nelson Dias, Pavan Kumar CV, Mike
Heskitt and Sean Putman for their support,

Prakash Pagadala, Premanand Suryavanshi, Rahul Rajan, Pranav Harikrishnan

,Nimisha Srivastava and for their support and helpful inputs.

And now - enjoy this study guide and let us know whether it helped you to
successfully apply Altair Compose to your projects.

Best regards

Dr. Matthias Goelke

On behalf of “The Altair University Team”

Disclaimer

Every effort has been made to keep the book free from technical as well as other mistakes.
However, publishers and authors will not be responsible for loss, damage in any form and
consequences arising directly or indirectly from the use of this book.

© 2018 Altair Engineering, Inc. All rights reserved. No part of this publication may be
reproduced, transmitted, transcribed, or translated to another language without the written
permission of Altair Engineering, Inc. To obtain this permission, write to the attention Altair
Engineering legal department at: 1820 Big Beaver, Troy, Michigan, USA, or call +1-248-614-
2400.

4
1. Introduction to Modelling of
Materials
The idea behind this ebook is to compute the stress state of a material when given
strain history data. By means of computation we will try to produce accurate stress vs
strain figures, that will sufficiently match the generally observed characteristics during
experiments. The first part of the guide will present 1D material modelling with basics,
but also a more advanced Ramberg-Osgood hardening model. Next, we will create
basic material models for 3D, where tensor calculus will be applied. The capstone
exercise of this guide will be a simple nonlinear quasi-static FEM example, where we
will investigate elastoplasticity in 1D two-truss element system.

Figure 1.1 : Elastoplastic deformation in 1D.

5
6
2. One Dimensional Models
The one-dimensional models are used mainly in simple tension-compression
characteristics of the material. In 1D, we operate purely on scalars, it is useful for
understanding the basic concepts behind elastoplasticity and hardening. In the
following chapter we will consider different plastic deformation schemes, that are
observed mainly in metals and alloys. We will skip the viscoelastic and elasto-
viscoplastic material models (rate-dependent material models), which are on the
other hand used for polymers and organic materials.

2.1. Hyperelastic Material

Hyperelastic materials are mathematically the simplest ones to model. Although in
metals and alloys the stiffness of the material is constant (linear stress-strain relation),
the general case includes an arbitrary (nonlinear) relationship between stress and
strain (variable stiffness). Based on the general case, the constitutive equation for
hyperelastic materials is derived.

Figure 2.1: Strain energy density definition

7
We introduce the strain energy density, which is defined as an elastically stored
energy in a deformed body (see Figure 1.1):

𝑑𝑑𝜓𝜓 = 𝜎𝜎𝑑𝑑𝜀𝜀

Where:

𝑑𝑑𝜓𝜓 – inifinitesimal strain energy density,

𝜎𝜎 – stress,

𝑑𝑑𝜀𝜀 – infinitesimal strain.

In general, we obtain:

𝜀𝜀
𝜓𝜓 = � 𝑑𝑑𝜓𝜓
0

Since the main features of the hyperelastic material are:

• Reversibility – there are no plastic (irreversible) deformations,

• Path independency – the stress-strain path is not dependent on strain rate
(how fast the strain deformation is applied). Dependency on strain rate is
included in viscoelastic and elasto-viscoplastic material.

Fulfilment of these conditions means that loading and unloading always takes place
along the same path – there is no energy loss. To derive the constitutive equation, we
will consider the dissipation inequality (Clausius-Duhem inequality), which ensures
that there is always either an energy loss due to dissipation or no energy loss. The
energy of the system cannot be increased due to the law of thermodynamics. The
equation is presented as a difference between stress power and rate of change of
energy density:

Ɗ = 𝜎𝜎𝜀𝜀̇ − 𝜓𝜓̇ ≥ 0

Where:

Ɗ – dissipation power,

𝜀𝜀̇ – strain rate,

8
𝜓𝜓̇ – strain energy density rate,

𝜎𝜎𝜀𝜀̇ – stress power.

For the elastic case, this equation returns 0 due to lack of dissipation. By expanding
the energy density rate by means of the chain rule we obtain:

𝜕𝜕𝜓𝜓 𝜕𝜕𝜓𝜓 𝜕𝜕𝜕𝜕 𝜕𝜕𝜓𝜓

𝜓𝜓̇ = = = 𝜀𝜀̇
𝜕𝜕𝑡𝑡 𝜕𝜕𝜕𝜕 𝜕𝜕𝑡𝑡 𝜕𝜕𝜕𝜕

By substituting this form of energy density rate into the dissipation inequality, we
finally get the constitutive relation:

𝜕𝜕𝜕𝜕
𝜎𝜎 =
𝜕𝜕𝜕𝜕

For the linear case, where the stiffness is constant, the energy function is modelled
as:

1
𝜓𝜓(𝜀𝜀) = 𝛦𝛦𝜀𝜀 2
2

Where:

𝛦𝛦 – Young’s modulus.

The stress equation follows then as:

𝜎𝜎 = 𝛦𝛦𝛦𝛦

The above relation can be represented with a single spring in terms of the rheological
model:

Figure 2.2: Rheological model for a linear elastic material.

9
2.1.1. Exercise 1: Modelling the Linear Elastic Material

All the numerical examples in this book will be based on known strain data. The strain
data will not be given as a function of time (experimental or simulation results), but
as an array of values discretized in time. That means a suitable timestep should be
chosen to obtain sufficient accuracy. In the case of linear elastic material, calculation
is very straight forward, since the only equation used is the σ = Εε relation.

1. Open Altair Compose and in the editor window type in:

to reset the variables and close all the plots. Additionally, for completeness,
provide a comment for the units used in the exercise. As you may notice, this
is the unit set used most frequently in HyperWorks.

2. Save your file as Elastic.oml. Remember to create a working directory in which

you will be saving all your files. You can name it as “Material_Modelling”.

3. Provide all relevant parameters for the calculation. This includes time
discretization (timestep), material data (in this case only Young Modulus E).

4. The time step (tstep) and number of timesteps (n_tstep) should be defined.
Given these two, the end_time can be calculated. We also need to provide an
array of all the time points t, since we will need it for plotting the stress and
strain evolution in time, it should start with t = 0 and end with t = end_time.
The size of t is equal to n_tstep+1.

10
5. For this exercise we will consider three different strain histories:
a) Linearly increasing
b) Nonlinearly increasing
c) Sinusoidal
We let the user choose which loadcase should be considered. For that, we ask
for a number between 1 and 3 and store the input in the loadcase variable.

6. The strain data (input) will be generated inside our code. First of all, create an
array of zeroes and store it in eps variable.

The first entry in the brackets corresponds to number of rows and the second
to number of columns.

7. Now the strain data array will be updated according to the number chosen by
the user. A switch-case conditional statement will be applied.

For more information regarding the switch-case statement, please refer to

the online help documentation (press F1 when in Altair Compose).

11
a) The first case is defined by using the linspace command. First value is 0, the
last (maximum) 0.002, which corresponds to 0,2% of strain. The number of
entries is n_tstep+1.
b) To obtain the nonlinearly increasing strain, we use a for loop to input a
quadratic function. We update all the values in eps except the 1 – we leave 0
as it is. Therefore the loop starts with i = 2. The value (i-1)/(n_tstep+1) is
always between 0 and 1. It is then squared and multiplied by 0.001 so that
the maximum strain is 0,1 %.
c) The third case is analogous to the second one. Inside the sin function an
interval between 0 and 4*pi is considered to model two full cycles.

8. Now we are ready to compute the stress values. For this we simply need to
multiply each entry in the strain array eps with Young`s Modulus E. Use the

command .* to multiply element-wise. We store the stress array in the sig

variable.
9. Now it is time for the postprocessing part. We will plot three figures:
a) Strain vs Time
b) Stress vs Time
c) Stress vs Strain

12
We use subplot command to obtain three figures in a row in single window.

10. When done, save your file and run it. In command window, you will be asked
to input the loadcase number.

13
Type 1 and press enter. Three figures will pop up. Expand the plotting window
by clicking on the marked icon.

Figure 2.3: Exercise 1 results

After expansion you will be able to see the plots with their titles.

Figure 2.4: Exercise 1 Results

Notice that the first plot is presenting the array we defined with linspace
command. On the second figure the stress value is presented, which is a result

14
of multiplying the strain array with E = 210000. We can check the correctness
of this plot by taking the last strain value:
𝛦𝛦 ∗ 𝜀𝜀 = 210000 ∗ 0.002 = 420 𝛭𝛭𝑃𝑃𝑃𝑃
The last plot is a linear relation between stress and strain – due to constant
thickness.
11. Try out two remaining load cases by rerunning the script. You will notice that
although the strain and stress is not linear in time, the stress vs strain figure
will always be the same. This is due to hyperelastic material’s property of
reversibility and path-independency.

Figure 2.5 : Exercise 1 results for 2 loadcase

Figure 2.6 : Exercise 1 results for 3 loadcase

15
2.2. Ideal Elastoplastic Material
The main characteristic of an elastoplastic material is its two-stage stress-strain
relationship. In addition to the linear elastic material from the previous chapter we
introduce a yield stress σy , which defines a threshold value for plastic deformation.
When the stress state is not exceeding the yield stress, then we speak about an elastic
behaviour. However, when the yield stress is reached (exceeded) plastic (irreversible)
deformation is observed. In a general elastoplastic case, the stress value can go above

Figure 2.7: Ideal elastoplastic material stress-strain relation.

yield stress and a certain stiffness in the plastic deformation state can be observed.
However, right now let’s concentrate on the ideal elastoplasticity, that means the
stiffness in the plastic state is 0 and a stress value is always equal to σy (see figure
below).

As you may notice from the figure above, after the plastic deformation has occurred
(stage 2), the unloading path (stage 3) is not the same as before (stage 1). It means
that some amount of strain energy was dissipated (this amount is equal to the area of
the parallelogram enclosed by stages 1,2 and 3). The dissipation inequality will return
a positive value in that case. It can be obtained by splitting the total strain into the
elastic εe and plastic εp part and considering the fact that strain energy depends only
on the elastic part of the strain.

16
𝜕𝜕𝜕𝜕 𝑒𝑒
𝐷𝐷 = 𝜎𝜎�𝜀𝜀 𝑒𝑒̇ + 𝜀𝜀 𝑝𝑝̇ � − 𝜀𝜀 ̇ ≥ 0
𝜕𝜕𝜀𝜀 𝑒𝑒

Where:

𝜀𝜀 𝑒𝑒̇ – elastic strain rate,

𝜀𝜀 𝑝𝑝̇ – plastic strain rate.

Since:

𝜕𝜕𝜕𝜕
𝜎𝜎 =
𝜕𝜕𝜀𝜀 𝑒𝑒

We obtain a reduced dissipation inequality:

𝐷𝐷 𝑟𝑟𝑟𝑟𝑟𝑟 = 𝜎𝜎𝜀𝜀 𝑝𝑝̇ ≥ 0

The rheological model for ideal elastoplasticity is composed of a spring and frictional
device, which depicts the plastic deformation.

Figure 2.8: Rheological model for ideal elastoplasticity.

2.2.1. Algorithm for Ideal Elastoplasticity

To derive an algorithm for the ideal elastoplasticity, an evolution equation describing

the plastic deformation is needed. We will use the reduced dissipation inequality
introduced before and apply the postulate of maximum dissipation.

17
max�𝐷𝐷 𝑟𝑟𝑟𝑟𝑟𝑟 �

Additionally, we need to introduce a so called yield function φ(σ), which is responsible

for checking if the current stress value is admissible, i.e. if the elastic part of
deformation is taking place. If the elastic part is linear, it has a form:

𝛷𝛷(𝜎𝜎) = |𝜎𝜎| − 𝜎𝜎𝑦𝑦 ≤ 0

Where:

𝜎𝜎𝑦𝑦 – yield stress.

The way to maximize the reduced dissipation function is to apply a mathematical

optimization method - the Lagrange multiplier method. We will not go into detail with
the calculation process, but present a final evolution equation:

𝜕𝜕𝜕𝜕
𝜀𝜀 𝑝𝑝̇ = 𝜆𝜆̇ = 𝜆𝜆̇𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠(𝜎𝜎)
𝜕𝜕𝜕𝜕

Where:

𝜆𝜆̇ – plastic multiplier (Lagrange inequality multiplier).

The evolution equation is accompanied by following conditions:

𝜆𝜆̇ ≥ 0, 𝛷𝛷 ≤ 0, 𝜆𝜆̇𝛷𝛷 = 0

The plastic multiplier λ is basically an absolute value of plastic deformation. The

inequality conditions are so-called KKT conditions and they indicate whether elastic
deformation takes place (𝜆𝜆̇ > 0 and 𝛷𝛷 = 0) or plastic deformation (𝜆𝜆̇ = 0 and 𝛷𝛷 <
0).

The evolution equation for plastic deformation rate is a basis for numerical algorithm.
This is a first order differential equation. Normally, when hand calculations are done,
a definite integral is applied to obtain the plastic strain value. However, in our case,
we want to apply the method into a computer and that means we need to use a
numerical integration scheme. The numerical integration schemes use time
discretized function data and function steepness information to obtain the solution

18
at the next time point. Our choice is the implicit Euler backward scheme, which is
based on steepness information at the n+1 point:

𝑎𝑎𝑛𝑛+1 = 𝑎𝑎𝑛𝑛 + ∆𝑡𝑡 𝑎𝑎̇ (𝑡𝑡𝑛𝑛+1 , 𝑎𝑎𝑛𝑛+1 )

On the figure below, one step of Euler backward method is depicted. As you may
notice, each step introduces an error when compared to the continuous solution. Of
course, the smaller the timestep (finer discretization), the more accurate the solution.

Figure 2.9: Implicit Euler Backward method.

In the ideal elastoplastic application, the following scheme is used:

𝑝𝑝 𝑝𝑝 𝑝𝑝 𝑝𝑝 𝑝𝑝
𝜀𝜀𝑛𝑛+1 = 𝜀𝜀𝑛𝑛 + ∆𝑡𝑡𝜀𝜀̇𝑛𝑛+1 = 𝜀𝜀𝑛𝑛 + ∆𝑡𝑡𝜆𝜆̇𝑛𝑛+1 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠(𝜎𝜎𝑛𝑛+1 ) = 𝜀𝜀𝑛𝑛 + ∆𝜆𝜆𝑛𝑛+1 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠(𝜎𝜎𝑛𝑛+1 )

Where:

𝑝𝑝
𝜀𝜀𝑛𝑛 – known plastic deformation from the n-th timestep (previous timestep),

𝑝𝑝
𝜀𝜀𝑛𝑛+1 – unknown plastic deformation, to be calculated in the current timestep,

𝜆𝜆̇𝑛𝑛+1 – plastic multiplier at the current timestep.

Since at n+1 point we know the total strain rate 𝜀𝜀̇𝑛𝑛+1, it can be split into the plastic
and elastic component. The elastic component is calculated according to the scheme:

𝑒𝑒 𝑝𝑝 𝑝𝑝
𝜀𝜀𝑛𝑛+1 = 𝜀𝜀𝑛𝑛+1 − 𝜀𝜀𝑛𝑛+1 = 𝜀𝜀𝑛𝑛+1 − 𝜀𝜀𝑛𝑛 − ∆𝜆𝜆𝑛𝑛+1 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠(𝜎𝜎𝑛𝑛+1 )

19
Where:

𝑡𝑡 𝑝𝑝
𝜀𝜀𝑛𝑛+1 = 𝜀𝜀𝑛𝑛+1 − 𝜀𝜀𝑛𝑛 , the trial strain value.

The trial value is computed according to known total strain (input) and plastic strain
𝑡𝑡
from the previous step. The following approach is to compute trial stress 𝜎𝜎𝑛𝑛+1 and
𝑡𝑡
trial yield function 𝛷𝛷𝑛𝑛+1 values. The remaining part of the equation ∆𝜆𝜆𝑛𝑛+1 𝑠𝑠𝑖𝑖𝑖𝑖𝑖𝑖(𝜎𝜎𝑛𝑛+1 )
is called the correction part and is only needed when plastic deformation appears.

The trial yield function is used to check if plastic deformation appears.

If:

𝑡𝑡 𝑡𝑡 |
𝛷𝛷𝑛𝑛+1 = |𝜎𝜎𝑛𝑛+1 − 𝜎𝜎𝑦𝑦 > 0

Then, the (actual) yield function must be equal to 0 according to KKT conditions, so:

𝑡𝑡
𝛷𝛷𝑛𝑛+1 = 𝛷𝛷𝑛𝑛+1 − 𝐸𝐸∆𝜆𝜆𝑛𝑛+1 = 0

So, in case of plastic deformation the correction parameter is calculated as

𝑡𝑡
𝛷𝛷𝑛𝑛+1
∆𝜆𝜆𝑛𝑛+1 =
𝐸𝐸

If the trial yield function return value is smaller or equal to 0, then the trial values are
the ready solution for n+1 point.

This algorithm then gets the following form:

%Known:
𝑝𝑝
𝜀𝜀𝑛𝑛+1 , 𝜀𝜀𝑛𝑛

%Compute trial values (elastic prediction step):

𝑡𝑡 𝑝𝑝
𝜎𝜎𝑛𝑛+1 = 𝐸𝐸(𝜀𝜀𝑛𝑛+1 − 𝜀𝜀𝑛𝑛 )

𝑡𝑡 𝑡𝑡 |
𝛷𝛷𝑛𝑛+1 = |𝜎𝜎𝑛𝑛+1 − 𝜎𝜎𝑦𝑦

%Check if yielding occurs:

𝑡𝑡
If 𝛷𝛷𝑛𝑛+1 ≤0

𝑡𝑡
𝜎𝜎𝑛𝑛+1 = 𝜎𝜎𝑛𝑛+1

20
𝑝𝑝 𝑝𝑝
𝜀𝜀𝑛𝑛+1 = 𝜀𝜀𝑛𝑛

Else

𝑡𝑡
𝛷𝛷𝑛𝑛+1
∆𝜆𝜆𝑛𝑛+1 =
𝐸𝐸

𝑡𝑡
𝜎𝜎𝑛𝑛+1 = 𝜎𝜎𝑛𝑛+1 − 𝐸𝐸∆𝜆𝜆𝑛𝑛+1 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠(𝜎𝜎𝑛𝑛+1 )

𝑝𝑝 𝑝𝑝
𝜀𝜀𝑛𝑛+1 = 𝜀𝜀𝑛𝑛 + ∆𝜆𝜆𝑛𝑛+1 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠(𝜎𝜎𝑛𝑛+1 )

The algorithm can be visualised in the following manner:

Figure 2.10: Algorithm for the ideal elastoplastic material.

2.2.2. Exercise 2: Ideal Elastoplasticity

For the following exercise a successful completion of Exercise 1 is necessary.

The file Elastic.oml should be made available in the “Material_Modelling”
working directory.

1. Open Altair Compose and make sure you are in your „Material Modelling”
working directory. Since some of the parts of the following code will be similar
to the linear elastic example, we can organize them in separate files, called

21
functions. Let’s open Elastic.oml file by dragging it from the File Browser into
the Editor Window.

2. Make sure the loadcase definition section is defined as follows:

Now create a new file and name it loadcase.oml. This file will be used as an
external function that can be called from other .oml files.

3. First line in newly created file should start with function definition:

It always must start with a keyword function. The next position defines the
return variables, in our case it is a strain data array. After the equality sign we
have a function name and input arguments – load case number n and number
of timesteps n_tstep.

4. Next cut and paste the loadcase definition section into the function. Change
the switch loadcase line to switch n. Now we can introduce a fourth load

22
case to our function, which will consist of 4 linear parts. In the end your
function should look like the following:

Let’s create a new file for ideal elastoplasticity. Type in the new editor
window, You can save the file as IdealElastoplasticity.oml

5. Now that we already have 3 files, it will be convenient to organize them into
a project.

6. In the Project Browser, right click on the <Project: Untitled> position.

7. In the pop up menu, choose “Save Project”.

8. Save the project file as “MaterialModelling.aprj”

23
9. Now we can add our files to the project. Make sure all three files: Elastic.oml
loadcase.oml and IdealElastoplastic.oml are opened in the EditorWindow.

10. Right click on the bookmark of each file and in the pop up menu choose “Add
to project”

24
Repeat the procedure for the remaining two files, Notice that the
loadcase.oml file is expandable, it is because it contains a “loadcase” function.
A currently open file will be highlighted in the Project Browser.

11. Let’s continue with the IdealElastoplastic.oml file. We need to provide initial
parameters for our model. First we define time discretization parameters and
the material parameters.
We choose a steel material with yield strength of 180 MPa.

12. The loadcase choice will be made by user within I/O interaction.

13. Now the loadcase function needs to be called. It will return an array of strain
data based on the passed loadcase number and number of time steps to
define the length of the array:

Initialize the plastic strain array and stress array with zeroes:

25
14. Having all the necessary inputs setup, we can provide the algorithm code.
Let’s recall the main steps:
1. Compute the trial values (elastic prediction step)
2. Check if for given trial stress we exceed the yield value
2.1. If not: Assign trial values to the actual solution
2.2. If yes: Calculate the correction parameter lambda and apply the
correction to the trial values to obtain the solution.

15. Plot the results. In addition to the figures from Exercise 1 we also want to plot
plastic strain vs time( this time with subplot 2x2 window ). Additionally, strain
values are in this case represented as percentage, therefore eps array is
scaled by 100 and [%] unit information is provided.

26
16. Save your file. Now we can execute the model. First, the function file needs
to be run in order to be usable. View the loadcase.oml file and press “Run”
button or ctrl+E. After that, switch to the IdealElastoplastic.oml file and
execute it. Try out all the load cases.

27
Loadcase 1:

In the first loadcase, the strain is increasing linearly. However the stress value
stops increasing when it reaches the yield value. At the same time point
plastic strain starts increasing which indicates that the plastic deformation
phase began.

Loadcase 2:
In the second loadcase the situation is analogous in comparison to the first
one.

28
Loadcase 3:

The third loadcase is cyclic loading with both tension and compression
behaviour. As the strain vs time curve is a pure sinusoid, the stress vs time
curve is sliced as the values exceed the yield value in both tension and
compression. Now let’s look directly at the fourth figure: Plastic Strain vs
Time. You can notice that the plastic strain is undergoing change only during
the time periods when stress value is equal to yield. For the rest of the time
it remains constant and is not going back to zero (irreversibility). The typical
path-dependency is depicted on the third figure, where you can clearly notice
that the elastic unloading does not follow the same path as loading with
plastic deformation. If the sinusoidal load would be continued for a longer
period of time, the stress-strain characteristic will not change but continue
following the parallelogram path.

29
Loadcase 4:

Loadcase 4 presents what happens when the loading path is less regular
(unsymetric with respect to time axis). Although the strain never reaches
negative value, we actually notice compressive values of stress. This is
because the amplitude between peaks in the first figure is big enough to
decrease the stress value by more than 180 MPa (yield value). The Stress-
Strain relationship also becomes irregular and each elastic deformation
follows a different path.

2.3. Elastoplastic Material with Linear Isotropic

Hardening
The isotropic hardening behaviour is observed as an expansion of elastic region while
plastic deformation happens. It means that the yield stress value is increased in both

30
tension and compression symmetrically (see fig. below), when a body is undergoing
plastic deformation with increasing stress value, as opposed to ideal elastoplasticity.

Figure 2.11: Isotropic hardening in elastoplastic material.

Look at the figure above. When the first stage of deformation takes place, plastic
yielding begins when the 𝜎𝜎𝑦𝑦 is reached. The plastic deformation causes the stress to
increase further (we do not discuss about the ideal plastic deformation any further).
We stop deforming in the tensile direction at some point. The stress value is now
equal to 𝜎𝜎𝑦𝑦 yield stress + some value 𝑅𝑅1 . When we apply deformation in the
compressive direction, we notice that in the compression state, yield occurs only after
the −(𝜎𝜎𝑦𝑦 + 𝑅𝑅1 ) stress value is reached. The elastic deformation region has now
symmetrically expanded due to plastic deformation. This process of increasing the
elastic region will continue until the part will fail.

Isotropic hardening is included in the strain energy density definition by introducing

an additional variable 𝛼𝛼, which we can say is responsible for hardening strain.

𝜕𝜕𝜕𝜕 𝑒𝑒 𝜕𝜕𝜕𝜕
𝐷𝐷 = 𝜎𝜎𝜀𝜀̇ − 𝜀𝜀̇ − 𝛼𝛼̇ ≥ 0
𝜕𝜕𝜀𝜀 𝑒𝑒 𝜕𝜕𝜕𝜕

31
Where:

𝛼𝛼 – hardening strain-like variable.

Assuming that isotropic hardening is linear and defines the K stiffness during plastic
loading, the strain Energy density is expressed as a function of elastic strain and alpha:

1 1
𝜓𝜓 = 𝜓𝜓(𝜀𝜀 𝑒𝑒 , 𝛼𝛼) = 𝐸𝐸(𝜀𝜀 𝑒𝑒 )2 + 𝐾𝐾𝛼𝛼 2
2 2

Where:

𝐾𝐾 – the hardening stiffness.

When we consider the elastic deformation, then the dissipation is equal to zero. We
can formulate the reduced dissipation inequality (reduced by the elastic components)
which introduces the hardening stress 𝑅𝑅, that was presented on the previous figure.

𝐷𝐷 𝑟𝑟𝑟𝑟𝑟𝑟 = 𝜎𝜎𝜀𝜀̇ 𝑝𝑝 + 𝑅𝑅𝛼𝛼̇ ≥ 0

Where 𝑅𝑅 = −𝜕𝜕𝜕𝜕/𝜕𝜕𝜕𝜕 is introduced as a negative value. The hardening stress R now

needs to be included in the yield function definition:

𝛷𝛷 = |𝜎𝜎| − 𝐺𝐺(𝛼𝛼) = |𝜎𝜎| − (𝜎𝜎𝑦𝑦 − 𝑅𝑅) ≤ 0

Where:

𝐺𝐺(𝛼𝛼) – is a function containing the isotropic hardening behaviour (for example: linear
hardening: 𝐺𝐺(𝛼𝛼) = 𝜎𝜎𝑦𝑦 − 𝑅𝑅).

And after taking the partial derivative of the energy density function with respect to
α we get:

𝑅𝑅 = −𝐾𝐾𝛼𝛼

32
2.3.1. The Algorithm

The implicit Euler backward step uses the following evolution equations computed
with the Lagrangian optimization methods:

𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕
𝜀𝜀̇ 𝑝𝑝 = 𝜆𝜆̇ , 𝛼𝛼̇ = 𝜆𝜆̇
𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕

With the KKT conditions as usual for elastoplasticity:

𝜆𝜆̇ ≥ 0, 𝛷𝛷 ≤ 0, 𝜆𝜆̇𝛷𝛷 = 0

We can expand the algorithm formulae with the alpha strain variable:

𝑝𝑝 𝑝𝑝
𝜀𝜀𝑛𝑛+1 = 𝜀𝜀𝑛𝑛 + ∆𝜆𝜆𝑛𝑛+1 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠(𝜎𝜎𝑛𝑛+1 )

𝛼𝛼𝑛𝑛+1 = 𝛼𝛼𝑛𝑛 + ∆𝜆𝜆𝑛𝑛+1

By considering the trial function:

𝑡𝑡 𝑡𝑡 | 𝑡𝑡 |
𝛷𝛷𝑛𝑛+1 = |𝜎𝜎𝑛𝑛+1 − 𝐺𝐺(𝛼𝛼𝑛𝑛 ) = |𝜎𝜎𝑛𝑛+1 − �𝜎𝜎𝑦𝑦 − 𝑅𝑅𝑛𝑛 �

Where:

𝐺𝐺(𝛼𝛼𝑛𝑛+1 ) – is a function defining the type of isotropic hardening (in that case linear)

By considering the yield function:

𝛷𝛷𝑛𝑛+1 = |𝜎𝜎𝑛𝑛+1 | − (𝜎𝜎𝑦𝑦 − 𝑅𝑅𝑛𝑛+1 )

Where we substitute:

𝑅𝑅𝑛𝑛+1 = 𝑅𝑅𝑛𝑛 + 𝐾𝐾∆𝜆𝜆𝑛𝑛+1 = 𝐾𝐾(𝛼𝛼𝑛𝑛 + 𝜆𝜆𝑛𝑛+1 )

𝑡𝑡 |
|𝜎𝜎𝑛𝑛+1 | = |𝜎𝜎𝑛𝑛+1 − 𝐸𝐸𝛥𝛥𝛥𝛥𝑛𝑛+1

The yield function is obtained as:

𝑡𝑡
𝛷𝛷𝑛𝑛+1 = |𝛷𝛷𝑛𝑛+1 | − ∆𝜆𝜆𝑛𝑛+1 (𝐸𝐸 + 𝐾𝐾)

We derive the correction variable according to KKT conditions – when the plastic
loading takes place we get ∆𝜆𝜆𝑛𝑛+1 > 0 and 𝛷𝛷𝑛𝑛+1 = 0, so:

33
𝑡𝑡
𝛷𝛷𝑛𝑛+1
∆𝜆𝜆𝑛𝑛+1 =
𝐸𝐸 + 𝐾𝐾

Now we can present the full algorithm:

%Known:
𝑝𝑝
𝜀𝜀𝑛𝑛+1 , 𝜀𝜀𝑛𝑛 , 𝜶𝜶𝒏𝒏

%Compute trial values (elastic prediction step):

𝑡𝑡 𝑝𝑝
𝜎𝜎𝑛𝑛+1 = 𝐸𝐸(𝜀𝜀𝑛𝑛+1 − 𝜀𝜀𝑛𝑛 )

𝑡𝑡 𝑡𝑡 |
𝛷𝛷𝑛𝑛+1 = |𝜎𝜎𝑛𝑛+1 − (𝜎𝜎𝑦𝑦 + 𝑲𝑲𝜶𝜶𝒏𝒏 )

%Check if yielding occurs:

𝑡𝑡
If 𝛷𝛷𝑛𝑛+1 ≤0

𝑡𝑡
𝜎𝜎𝑛𝑛+1 = 𝜎𝜎𝑛𝑛+1

𝑝𝑝 𝑝𝑝
𝜀𝜀𝑛𝑛+1 = 𝜀𝜀𝑛𝑛

𝜶𝜶𝒏𝒏+𝟏𝟏 = 𝜶𝜶𝒏𝒏

Else

𝑡𝑡
𝛷𝛷𝑛𝑛+1
∆𝜆𝜆𝑛𝑛+1 =
𝐸𝐸 + 𝑲𝑲

𝑡𝑡
𝜎𝜎𝑛𝑛+1 = 𝜎𝜎𝑛𝑛+1 − 𝐸𝐸∆𝜆𝜆𝑛𝑛+1 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠(𝜎𝜎𝑛𝑛+1 )

𝑝𝑝 𝑝𝑝
𝜀𝜀𝑛𝑛+1 = 𝜀𝜀𝑛𝑛 + ∆𝜆𝜆𝑛𝑛+1 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠(𝜎𝜎𝑛𝑛+1 )

𝜶𝜶𝒏𝒏+𝟏𝟏 = 𝜶𝜶𝒏𝒏 + ∆𝝀𝝀𝒏𝒏+𝟏𝟏

The additional parts with respect to ideal elastoplasticity is marked in bold text. You
can notice that this is mainly keeping track of the hardening strain variable 𝛼𝛼 and using
it to expand the elastic region by means of the yield function.

34
Figure 2.12: Algorithm step for elastoplasticity with isotropic hardening.

The reason why, on the figure above, the trial yield function exceeds the actual stress-
strain curve (blue dashed line) is because it considers alpha strain from the previous
time step. However, the delta lambda correction parameter is divided by a sum of E
and K, so that the actual strain difference between fully elastic strain component and
total strain is the distance from the blue dashed line to the trial step line.

2.3.2. Exercise 3: Elastoplasticity with Isotropic Hardening

To develop this exercise we will use the code from the Exercise 2, since a case with
isotropic hardening is basically an extension to ideal elastoplasticity code.
1. Open the MaterialModelling.aprj project file. Click on the yellow folder icon
in the top-left corner. In the pop up window, change the File type to .aprj:

35
Then select your MaterialModelling.aprj file and open it. You will notice that
all the three files, that were created until now, pop up in the Project Browser
together with the project name.

2. Create a new file. Copy-paste the whole code from the file
IdealElastoplasticity to your newly created file by using shortcuts (windows):
ctrl+A to select, ctrl+C to copy and in the new file window: ctrl+V to paste.
3. Save your new file as EPIsoHardening.oml.
4. Add the new file to our project: Right click on the EPIsoHardening.oml
bookmark and select “Add to Project”.

36
5. Now it is time to extend the code. First provide an additional parameter for
the hardening stiffness – K:

6. Initialize the alpha array with zeros:

7. Extend the algorithm with alpha update and K stiffness:

Now you can see that the yield function (line 31) is extended to check the
yielding with respect to the hardening stress R additionally, which is stored in

37
alpha array components multiplied by K stiffness: R = -Kα. When the yielding
is not occuring (line 40) then alpha component remain the same just like the
plastic strain. The correction parameter (line 45) is also updated with respect
to E + K additionally, so the correction is different from the ideal
elastoplasticity case.

8. Now let’s save the file and add another load case to the loadcase.oml
function. Add case 5 by copying the case 4 code and modifying the vertices
magnitudes:

9. Add also case 6 by copying the case 3 code and changing the amplitude of the
oscillation to 0.002:

10. Save the loadcase.oml file. Update the sentence in the input definition to let
the user choose number between 1 and 6:

11. Our last step would be to add a new plot to show the alpha parameter
evolution. In order to do that expand the subplot domain in all the figures to
2 x 3:

Remember to extend this domain for all plots!

38
12. Now add a fifth plot to show alpha evolution in time:

13. Save the file and execute the file – press ctrl+E. Choose the load case nr 5. The
resultant plots should look as follows:

Let’s first take a look at the first graph: This is our strain evolution defined in loadcase
function. It has an increasing amplitude in order to let the elastic region expand in
each cycle. Looking at the Stress vs Time graph let us notice that each of the four
loading stages cause plastic deformation – we can see it by clear change in the stress
increase rate. The Stress vs Strain graph is probably the most suitable to notice the
isotropic hardening phenomenon. You can see that in each loading stage bigger stress
value was needed to cause plastic deformation and that this value expanded
symmetrically. Let’s also take a look at the last figure – Alpha vs Time. You can see
that the alpha variable is cummulative – it increases it’s value whenever plastic
deformation happens and consequently isotropic hardening takes place. It’s final
value is equal to sum of all plastic strain changes plotted on the Figure 2.3.

39
Let’s take a look at loadcase nr 6:

This time the load amplitude does not expand within the next cycle. Therefore on the
Stress vs Strain curve you may notice that in the compression state the first cycle has
larger plastic deformation, but less elastic strain. Whereas the next cycle presents
greater elastic deformation, so the plastic component is smaller.

2.4. Elastoplastic Material with Linear

Kinematic Hardening
Kinematic hardening is a second possibility of how the elastic region can change
during plastic deformation. Instead of expanding, a switch in position of elastic region
is observed. Let’s look at the figure below:

40
Figure 2.53: Elastoplastic material with kinematic hardening.

In the initial loading (starting from the centre) the plastic deformation occurs after the
yield stress 𝜎𝜎𝑦𝑦 is reached. During that, the stress value reaches 𝜎𝜎𝑦𝑦 + 𝐵𝐵, where B is
called the Back stress and its value is equal to the displacement of elastic region due
to plastic deformation. The elastic stress range is still 2𝜎𝜎𝑦𝑦 , but the position has
displaced towards tensile stress by B value.

2.5. Elastoplastic Material with Linear Isotropic

and Kinematic Hardening
In experiments, most of the elastoplastic materials, mainly metals and alloys, present
hardening behaviour which can be classified somewhere between isotropic and
kinematic hardening. We then talk about the percentage of participation of each
hardening model in the material characteristic.

For kinematic hardening, like isotropic, additional variable 𝛽𝛽 with units of strain is
introduced, which is responsible for the displacement of elastic region. In the 1D case
this is a scalar value, in 3D space however, it is a tensor and will be discussed in further
chapters.

41
Our strain energy density function is further expanded to include 𝛽𝛽:

𝜓𝜓 = 𝜓𝜓(𝜀𝜀 𝑒𝑒 , 𝛼𝛼, 𝛽𝛽)

and the dissipation inequality obtains additional term of partial derivation of strain
energy density with respect to 𝛽𝛽:

𝜕𝜕𝜕𝜕 𝑒𝑒 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕

𝐷𝐷 = 𝜎𝜎𝜀𝜀̇ − 𝑒𝑒
𝜀𝜀̇ − 𝛼𝛼̇ − 𝛽𝛽̇ ≥ 0
𝜕𝜕𝜀𝜀 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕

Getting rid of the elastic deformation part from the dissipation inequality (during
elastic loading there is no energy loss – no dissipation), the reduced form of the
equation is as follows:

𝐷𝐷 𝑟𝑟𝑟𝑟𝑟𝑟 = 𝜎𝜎𝜀𝜀̇ 𝑝𝑝 + 𝑅𝑅𝛼𝛼̇ + 𝐵𝐵𝛽𝛽̇ ≥ 0

The Back stress 𝐵𝐵 = −𝜕𝜕𝜓𝜓/𝜕𝜕𝜕𝜕 is introduced, which is equal to the vertical switch of
the elastic region on the stress-strain figure. In the 3D case, the back stress is a tensor.
Now we need to include the Back stress into our yield function:

𝛷𝛷 = |𝜎𝜎 + 𝐵𝐵| − (𝜎𝜎𝑦𝑦 − 𝑅𝑅) ≤ 0

The hardening stress R is combined with yield stress to count in the expansion of
elastic region and the Back Stress B is combined with the current stress to enable the
unsymmetric stress region (different tensile and compressive stress limits). The 𝜎𝜎 +
𝐵𝐵 sum will be later referred to as the relative stress ξ.

Since we need to take care of both hardening methods, the strain energy density
function for the linear isotropic and linear kinematic hardening looks as follows:

1 1 1
𝜓𝜓 = 𝜓𝜓(𝜀𝜀 𝑒𝑒 , 𝛼𝛼, 𝛽𝛽) = 𝐸𝐸(𝜀𝜀 𝑒𝑒 )2 + 𝑟𝑟𝑟𝑟𝛼𝛼 2 + (1 − 𝑟𝑟)𝐾𝐾𝛽𝛽 2
2 2 2

Where 𝑟𝑟 ∶= [0,1] scalar is introduced, to define the percentage of participation of

isotropic hardening in the material model.

The hardening and Back stress are getting following formulas in the end:

𝑅𝑅 = −𝑟𝑟𝑟𝑟𝛼𝛼, 𝐵𝐵 = −(1 − 𝑟𝑟)𝐾𝐾𝛽𝛽

42
2.5.1. The Algorithm

The implicit Euler backward step uses the following evolution equations computed
with the Lagrangian optimization methods:

𝜕𝜕𝛷𝛷 𝜕𝜕𝛷𝛷 𝜕𝜕𝛷𝛷

𝜀𝜀̇ 𝑝𝑝 = 𝜆𝜆̇ , 𝛼𝛼̇ = 𝜆𝜆̇ , 𝛽𝛽̇ = 𝜆𝜆̇
𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕

With the KKT conditions as usual for elastoplasticity:

𝜆𝜆̇ ≥ 0, 𝛷𝛷 ≤ 0, 𝜆𝜆̇𝛷𝛷 = 0

The equations are then as follows:

𝑝𝑝 𝑝𝑝
𝜀𝜀𝑛𝑛+1 = 𝜀𝜀𝑛𝑛 + ∆𝜆𝜆𝑛𝑛+1 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠(𝜎𝜎𝑛𝑛+1 + 𝐵𝐵𝑛𝑛+1 )

𝛼𝛼𝑛𝑛+1 = 𝛼𝛼𝑛𝑛 + ∆𝜆𝜆𝑛𝑛+1

𝛽𝛽𝑛𝑛+1 = 𝛽𝛽𝑛𝑛 + 𝛥𝛥𝛥𝛥𝑛𝑛+1 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠(𝜎𝜎𝑛𝑛+1 + 𝐵𝐵𝑛𝑛+1 )

Now we need to introduce the afore mentioned relative stress, which is included in
the sign function in the previous equations:

𝜉𝜉𝑛𝑛+1 = 𝜎𝜎𝑛𝑛+1 + 𝐵𝐵𝑛𝑛+1

The relative stress does not have a critical meaning for the algorithm, it is just included
to clear up the algorithm a little bit by omitting using the sum of current and Back
stress multiple times.

The computation of the trial stress and strain is the same as in the pure isotropic
hardening case, only the trial yield function now includes the relative trial stress:

𝑡𝑡 𝑡𝑡
𝜉𝜉𝑛𝑛+1 = 𝜎𝜎𝑛𝑛+1 + 𝐵𝐵𝑛𝑛+1

𝑡𝑡 𝑡𝑡 |
𝛷𝛷𝑛𝑛+1 = |𝜉𝜉𝑛𝑛+1 − (𝜎𝜎𝑦𝑦 − 𝑅𝑅𝑛𝑛 )

The rest is calculated analogously, the final algorithm is:

%Known:
𝑝𝑝
𝜀𝜀𝑛𝑛+1 , 𝜀𝜀𝑛𝑛 , 𝛼𝛼𝑛𝑛 , 𝜷𝜷𝒏𝒏

%Compute trial values (elastic prediction step):

43
𝑡𝑡 𝑝𝑝
𝜎𝜎𝑛𝑛+1 = 𝐸𝐸(𝜀𝜀𝑛𝑛+1 − 𝜀𝜀𝑛𝑛 )

𝝃𝝃𝒕𝒕𝒏𝒏+𝟏𝟏 = 𝝈𝝈𝒕𝒕𝒏𝒏+𝟏𝟏 − (𝟏𝟏 − 𝒓𝒓)𝑲𝑲𝜷𝜷𝒏𝒏

𝜱𝜱𝒕𝒕𝒏𝒏+𝟏𝟏 = |𝝃𝝃𝒕𝒕𝒏𝒏+𝟏𝟏 | − �𝝈𝝈𝒚𝒚 + 𝒓𝒓𝒓𝒓𝜶𝜶𝒏𝒏 �

%Check if yielding occurs:

𝑡𝑡
If 𝛷𝛷𝑛𝑛+1 ≤0

𝑡𝑡
𝜎𝜎𝑛𝑛+1 = 𝜎𝜎𝑛𝑛+1

𝑝𝑝 𝑝𝑝
𝜀𝜀𝑛𝑛+1 = 𝜀𝜀𝑛𝑛

𝛼𝛼𝑛𝑛+1 = 𝛼𝛼𝑛𝑛

𝜷𝜷𝒏𝒏+𝟏𝟏 = 𝜷𝜷𝒏𝒏

Else

𝑡𝑡
𝛷𝛷𝑛𝑛+1
∆𝜆𝜆𝑛𝑛+1 =
𝐸𝐸 + 𝐾𝐾

𝑡𝑡
𝜎𝜎𝑛𝑛+1 = 𝜎𝜎𝑛𝑛+1 − 𝐸𝐸∆𝜆𝜆𝑛𝑛+1 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠(𝝃𝝃𝒕𝒕𝒏𝒏+𝟏𝟏 )

𝑝𝑝 𝑝𝑝
𝜀𝜀𝑛𝑛+1 = 𝜀𝜀𝑛𝑛 + ∆𝜆𝜆𝑛𝑛+1 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠(𝝃𝝃𝒕𝒕𝒏𝒏+𝟏𝟏 )

𝛼𝛼𝑛𝑛+1 = 𝛼𝛼𝑛𝑛 + ∆𝜆𝜆𝑛𝑛+1

𝜷𝜷𝒏𝒏+𝟏𝟏 = 𝜷𝜷𝒏𝒏 + ∆𝝀𝝀𝒏𝒏+𝟏𝟏 𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔(𝝃𝝃𝒕𝒕𝒏𝒏+𝟏𝟏 )

Again, the equations in bold are the newly introduced ones. We compute the relative
stress according to the kinematic hardening percentage and the signum function also
uses only the relative stress to consider the asymmetry of yielding in tension and
compression.

44
2.5.2. Exercise 4: Elastoplasticity with Linear Isotropic and
Kinematic Hardening

In this exercise we will further extend the code for isotropic linear hardening.

1. Open Altair Compose and open the MaterialModelling.aprj project. Open the
EPIsoHardening.oml file and copy the content to the new window in the
Editor. Save the new file as EPIsoKinHardening.oml.

2. Add the new file to our project by right clicking on the file bookmark. In the
project browser you should notice your newly created file:

Now let’s modify the code. We will not add any further loadcases in this exercise,
but will just work with the new EPIsoKinHardening.oml file.

3. Before the line with load_n input, add the r variable for defining the
isotropic/kinematic hardening percentage.
We let the user to choose a real number between 0 and 1:

45
4. Now we need to initialize the beta strain array. Since the “beta” name is not
available, let’s shorten it and use a name “bet”:

5. Expand the algorithm by adding following formulae:

Now our algorithm allows to freely choose the amount of isotropic and kinematic
hardening in the material model. The ksi_t relative stress is used in each sign
function and the beta variables are updated analogously to alpha and plastic
strain.

46
6. The last thing in the code is to add another plot to monitor the beta strain
evolution:

7. Now it is time to run the algorithm and visualize the results. Execute the file
with ctrl+E and try out following parameters:
1: Kinematic hardening, loadcase 5

47
2: 50% Isotropic hardening, loadcase 5

3: 100% Isotropic hardening, loadcase 5

Notice how in the first case on the Figure 2.2 the plastic loading takes the
same path in each cycle. Lack of isotropic hardening does not allow the elastic
region to expand, but just displaces the yielding value so that it always follows
this path. In the second case adding 50 % of isotropic hardening allows
expansion of elastic region, but still the compression and tension yield stress

48
values are not the same at the beginning and the end of one elastic
deformation phase. The last case is 100% isotropic hardening, that means the
elastic deformation phases are symmetric considering the yield stress values.
Notice that more the isotropic hardening that takes place, less are the plastic
deformation and higher the stress values (Figures 2.1,2.3). Now let’s consider
loadcase 6.

4: Kinematic hardening, loadcase 6:

49
5: 50% Isotropic hardening, loadcase 6:

6: 100% Isotropic hardening, loadcase 6:

50
2.6. Nonlinear Hardening – Ramberg-Osgood
Model
Until now we have always used the assumption that the stress-strain evolution in
plasticity is linear. In addition, we assumed that the transition from elastic to plastic
state does not have any additional features. In this chapter we will introduce a more
advanced way of modelling the plasticity in 1D. We will consider Ramberg-Osgood
model, which is frequently used and is an accurate method for many metals and
alloys.

Figure 2.64: Nonlinear hardening models. [1]

The linearity of the previously presented models was contained in the yield function
criterion:

𝛷𝛷(𝜎𝜎, 𝑅𝑅) = |𝜎𝜎| − �𝜎𝜎𝑦𝑦 − 𝑅𝑅�

Where:

𝐺𝐺(𝛼𝛼) = 𝜎𝜎𝑦𝑦 − 𝑅𝑅 = 𝜎𝜎𝑦𝑦 + 𝐾𝐾𝛼𝛼

51
Now we linearly introduce a 𝐺𝐺(𝛼𝛼) function for the Ramberg-Osgood model, proposed
by [2]:

𝐺𝐺(𝛼𝛼) = 𝜎𝜎𝑦𝑦 + 𝐶𝐶𝛼𝛼 𝑚𝑚

Where:

C, m – are constant parameters determined to fit the experimental data.

Originally the equation modelling the stress-strain relation has a following form:

1/𝑛𝑛
𝜎𝜎 𝜎𝜎
𝜀𝜀 = + 0,002 � �
𝛦𝛦 𝜎𝜎𝑦𝑦

Where the 0,002 multiplier taken from a 0,2% yield strain. The n variable is a strain
hardening exponent. The first part of the equation 𝜎𝜎/𝛦𝛦 is the elastic part of strain and
0,002(𝜎𝜎/𝜎𝜎𝑦𝑦 )1/𝑛𝑛 is the plastic part of the strain.

2.6.1. The Algorithm

Now we will derive the new formula for the quadratic isotropic hardening. To be more
precise, we will derive a new formula for the correction parameter 𝛥𝛥𝛥𝛥𝑛𝑛+1 .

Given the quadratic hardening function:

𝑚𝑚
𝐺𝐺(𝛼𝛼𝑛𝑛+1 ) = 𝜎𝜎𝑦𝑦 + 𝐶𝐶𝛼𝛼𝑛𝑛+1

We get the yield function criterion for plastic state as:

0 = 𝛷𝛷𝑛𝑛+1 = |𝜎𝜎𝑛𝑛+1 | − 𝐺𝐺(𝛼𝛼𝑛𝑛+1 ) = |𝜎𝜎𝑛𝑛+1 | − 𝐺𝐺(𝛼𝛼𝑛𝑛+1 ) − 𝐺𝐺(𝛼𝛼𝑛𝑛 ) + 𝐺𝐺(𝛼𝛼𝑛𝑛 )

But considering the stress equation:

𝑡𝑡 |
|𝜎𝜎𝑛𝑛+1 | = |𝜎𝜎𝑛𝑛+1 − 𝐸𝐸𝛥𝛥𝛥𝛥𝑛𝑛+1

And trial yield function:

𝑡𝑡 𝑡𝑡 |
𝛷𝛷𝑛𝑛+1 = |𝜎𝜎𝑛𝑛+1 − 𝐺𝐺(𝛼𝛼𝑛𝑛 )

We obtain the final equation:

52
𝑡𝑡
0 = 𝛷𝛷𝑛𝑛+1 = 𝛷𝛷𝑛𝑛+1 − 𝐸𝐸𝛥𝛥𝛥𝛥𝑛𝑛+1 − 𝐺𝐺(𝛼𝛼𝑛𝑛+1 ) + 𝐺𝐺(𝛼𝛼𝑛𝑛 )

After a substitution of G function and the relation 𝛼𝛼𝑛𝑛+1 = 𝛼𝛼𝑛𝑛 + 𝛥𝛥𝛥𝛥𝑛𝑛+1 , we obtain:

𝑡𝑡
0 = 𝛷𝛷𝑛𝑛+1 − 𝐸𝐸𝛥𝛥𝛥𝛥𝑛𝑛+1 − 𝐶𝐶((𝛼𝛼𝑛𝑛 + 𝛥𝛥𝛥𝛥𝑛𝑛+1 )𝑚𝑚 − 𝛼𝛼𝑛𝑛𝑚𝑚 )

The equation given above cannot be solved directly in the implicit Euler Backward
algorithm. Instead an iterative Newton-Raphson method is used to approximate the
solution for 𝛥𝛥𝛥𝛥𝑛𝑛+1 . The Newton-Raphson algorithm is additionally extended by 𝜔𝜔 ∈
(0,1) parameter to enable convergence. This approach is called Successive Over-
Relaxation (SOR) and helps converge for diverging functions [3], i.e.:

𝑓𝑓(𝑥𝑥) = |𝑥𝑥|𝑚𝑚 , 0 < 𝑚𝑚 < 0.5

updated by following procedure:

𝑥𝑥𝑛𝑛+1 = 𝑓𝑓𝑢𝑢𝑢𝑢𝑢𝑢𝑢𝑢𝑢𝑢𝑢𝑢 (𝑥𝑥𝑛𝑛 )

in the following manner:

𝑥𝑥𝑛𝑛+1 = (1 − 𝜔𝜔)𝑥𝑥𝑛𝑛 + 𝜔𝜔𝑓𝑓𝑢𝑢𝑢𝑢𝑢𝑢𝑢𝑢𝑢𝑢𝑢𝑢 (𝑥𝑥𝑛𝑛 )

The final algorithm obtains a following form:

%Known:
𝑝𝑝
𝜀𝜀𝑛𝑛+1 , 𝜀𝜀𝑛𝑛 , 𝜶𝜶𝒏𝒏

%Compute the hardening function value:

𝐺𝐺(𝛼𝛼𝑛𝑛 ) = 𝜎𝜎𝑦𝑦 + 𝐶𝐶𝛼𝛼𝑛𝑛𝑚𝑚

%Compute trial values (elastic prediction step):

𝑡𝑡 𝑝𝑝
𝜎𝜎𝑛𝑛+1 = 𝐸𝐸(𝜀𝜀𝑛𝑛+1 − 𝜀𝜀𝑛𝑛 )

𝑡𝑡 𝑡𝑡 |
𝛷𝛷𝑛𝑛+1 = |𝜎𝜎𝑛𝑛+1 − 𝐺𝐺(𝛼𝛼𝑛𝑛 )

%Check if yielding occurs:

𝑡𝑡
If 𝛷𝛷𝑛𝑛+1 ≤0

𝑡𝑡
𝜎𝜎𝑛𝑛+1 = 𝜎𝜎𝑛𝑛+1

𝑝𝑝 𝑝𝑝
𝜀𝜀𝑛𝑛+1 = 𝜀𝜀𝑛𝑛

53
𝜶𝜶𝒏𝒏+𝟏𝟏 = 𝜶𝜶𝒏𝒏

Else

%Solve the equation d_lambda using Newton-Rhapson iterations

𝑡𝑡
∆𝜆𝜆𝑛𝑛+1  0 = 𝛷𝛷𝑛𝑛+1 − 𝐸𝐸𝛥𝛥𝛥𝛥𝑛𝑛+1 − 𝐶𝐶((𝛼𝛼𝑛𝑛 + 𝛥𝛥𝛥𝛥𝑛𝑛+1 )𝑚𝑚 − 𝛼𝛼𝑛𝑛𝑚𝑚 )

𝑡𝑡
𝜎𝜎𝑛𝑛+1 = 𝜎𝜎𝑛𝑛+1 − 𝐸𝐸∆𝜆𝜆𝑛𝑛+1 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠(𝜎𝜎𝑛𝑛+1 )

𝑝𝑝 𝑝𝑝
𝜀𝜀𝑛𝑛+1 = 𝜀𝜀𝑛𝑛 + ∆𝜆𝜆𝑛𝑛+1 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠(𝜎𝜎𝑛𝑛+1 )

𝜶𝜶𝒏𝒏+𝟏𝟏 = 𝜶𝜶𝒏𝒏 + ∆𝝀𝝀𝒏𝒏+𝟏𝟏

%NewtonRhapson algorithm

𝑡𝑡
%Given 𝛼𝛼𝑛𝑛 , 𝛷𝛷𝑛𝑛+1 , 𝐸𝐸, 𝐶𝐶, 𝑚𝑚

%Initialize iteration parameters

𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚 = 30 %max number of iterations

𝑖𝑖 = 0 %current iteration counter

𝑡𝑡𝑡𝑡𝑡𝑡 = 10𝑒𝑒 − 5 %threshold value for convergence – R residual need to

decrease below tol

𝑑𝑑 = 0 %increment value for d_lambda

𝜔𝜔 = 0.2 %SOR parameter

%Initialize the d_lambda value to 0

∆𝜆𝜆𝑛𝑛+1 = 0

%Calculate the residual function (which need to tend to 0):

𝑡𝑡
𝑅𝑅 = 𝛷𝛷𝑛𝑛+1 − 𝐸𝐸𝛥𝛥𝛥𝛥𝑛𝑛+1 − 𝐶𝐶((𝛼𝛼𝑛𝑛 + 𝛥𝛥𝛥𝛥𝑛𝑛+1 )𝑚𝑚 − 𝛼𝛼𝑛𝑛𝑚𝑚 )

%Run the iteration loop until convergence criterion or maxiter has been reached

while |𝑅𝑅| > 𝑡𝑡𝑡𝑡𝑡𝑡 and 𝑖𝑖 < 𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚

𝑑𝑑𝑑𝑑
= −𝐸𝐸 − 𝑚𝑚𝑚𝑚(𝛼𝛼𝑛𝑛 + 𝛥𝛥𝛥𝛥𝑛𝑛+1 )𝑚𝑚−1
𝑑𝑑𝛥𝛥𝛥𝛥𝑛𝑛+1

−1
𝑑𝑑𝑑𝑑
𝑑𝑑 = − � � 𝑅𝑅
𝑑𝑑𝛥𝛥𝛥𝛥𝑛𝑛+1

54
𝑓𝑓 = 𝛥𝛥𝛥𝛥𝑛𝑛+1 + 𝑑𝑑 %update function

𝛥𝛥𝛥𝛥𝑛𝑛+1 = (1 − 𝜔𝜔)𝛥𝛥𝛥𝛥𝑛𝑛+1 + 𝜔𝜔𝑓𝑓 %update with SOR

𝑡𝑡
𝑅𝑅 = 𝛷𝛷𝑛𝑛+1 − 𝐸𝐸𝛥𝛥𝛥𝛥𝑛𝑛+1 − 𝐶𝐶((𝛼𝛼𝑛𝑛 + 𝛥𝛥𝛥𝛥𝑛𝑛+1 )𝑚𝑚 − 𝛼𝛼𝑛𝑛𝑚𝑚 ) %recalculate R

𝑖𝑖 = 𝑖𝑖 + 1 %increment the iteration counter

end

2.6.2. Exercise 5: Elastoplasticity with Ramberg-Osgood

Isotropic Hardening

This exercise requires successful completion of the Exercise 3: Elastoplasticity with

linear isotropic hardening.

For the purpose of this exercise we will create two .oml files. The first one will be an
extension to linear isotropic model and the second one will be a function file, which
will compute the d_lambda variable by means of Newton-Rhapson iterations.

1. Open the project file MaterialModelling.aprj. Create a new file and copy the
content of EPIsoHardening.oml file to its editor window.
2. Save the file and name it EPIsoHardeningRO.oml, where RO stands for
Ramberg-Osgood.
3. Let’s start with updating the Material parameters:

We update the value for K stiffness and introduce two new parameters:
C – the RO stiffness constant, which is included in the RO material model.

m – the RO exponent, which is applied to accumulated plastic strain α.

55
4. Next update will concern the algorithm itself, as we leave out the loadcase
choice and variable initialization as it was before. At the beginning of the for

loop, we compute the G function value at the current α:

We use the .^ operator to apply the power operation.

5. Now we need to update the computation of trial values, since the trial yield
function now uses the G function to compare with the trial stress.

6. In the if-else scope, we leave the elastic update part as it was before and
update the plastic correction step accordingly:

We need to obtain current alpha value as a scalar in order to pass it to the

newtonRhapson function (which later in this tutorial will be defined). The
newtonRhapson function needs to take all the parameters used for iterations,
that is: Trial yield function, current alpha, Young modulus and RO constants.
The newtonRhapson function returns the d_lambda value, which is next used

56
to update the stress, plastic strain and alpha. We leave the rest of the code as
it is.

7. Save your file and add it to the project using a right mouse click on the file
bookmark and choosing “Add to Project” option.

8. Now it is time to create newtonRhapson.oml function. Create a new file and

name it as before.
9. Add a function definition:

Hint: Remember to save the file with the same name as function name !

10. Inside the function scope, we first need to specify the RO method parameters:

57
The maxiter is set to 50, since most of the calls need between 40 and 50
iterations to obtain the convergence. Such a slow convergence rate is due to
a small value of relaxation parameter – the smaller its value, the slower the
convergence. However it helps in obtaining convergence for such functions
like the RO model. The iteration counter i and the increment variable d are
initialized to 0. The tolerance for the error is chosen as 10e-5. The relaxation
parameter w is set to 0.2, which means that only 20% of the increment will
come from the update function and 80% from previous value of d_lambda.

Before the iterations start, we need to initialize our target variable d_lambda
to starting value of 0 and the initial value for R residual, which is basically our
measure of error of the d_lambda value estimation.

11. Now we can start the while loop with two termination conditions:

In order for the while loop to continue running, both conditions need to be
satisfied. The residual (error) needs to be greater than the tolerance and the
number of iterations cannot exceed the limit value we set in the beginning. If
one of the conditions is not satisfied – either convergence has been reached
or the number of iterations are exceeded, the while loop will end and the
current d_lambda value will be returned.

12. At the beginning of the while loop we need to calculate the derivative of the
residual function.

58
𝑑𝑑𝑑𝑑
= −𝐸𝐸 − 𝑚𝑚𝑚𝑚(𝛼𝛼𝑛𝑛 + 𝛥𝛥𝛥𝛥𝑛𝑛+1 )𝑚𝑚−1
𝑑𝑑𝛥𝛥𝛥𝛥𝑛𝑛+1

There is however a small issue at the beginning of computation. For the first
call of the function (first occurrence of plastic strain) the αn is still equal to 0
and in the first iteration the d_lambda is equal to 0, the computation of the
following expression: 𝑚𝑚𝑚𝑚(𝛼𝛼𝑛𝑛 + 𝛥𝛥𝛥𝛥𝑛𝑛+1 )𝑚𝑚−1 will return −∞ for 𝑚𝑚 = 0.2.

Therefore it is necessary to include an if-else condition to take care of the

case with −∞:

The single expression (𝛼𝛼𝑛𝑛 + 𝛥𝛥𝛥𝛥𝑛𝑛+1 )𝑚𝑚−1 is initialized to 0, if the sum of α and
Δλ returns 0, if it doesn’t, the expression is computed in the normal way.

13. The rest of the algorithm is defined in the following way:

The derivative of the residual is stored in dR variable, the increment d is

computed and the update function f is obtained. However the actual
increment of the d_lambda variable is done with respect to the SOR method
– with use of ω relaxation parameters to get convergence.
In the end the residual R is recomputed and iteration counter i incremented.
Remember to include the end statements to close the appropiate scopes.

59
14. Save the file. Now we are ready to run the EPIsoHardeningRO.oml file. Make
sure your loadcase.oml, newtonRhapson.oml file are available in the
directory.

15. The results for the chosen loadcases are presented below:

Loadcase 1:

The simple loadcase with linearly increasing strain best depicts the idea
behind the RO method. The yielding is much more smooth, there is no vertex
at the transition between elastic and plastic deformation, which for sure
presents a more realistic figure of elastoplasticity.

60
Loadcase 4:

The figures are very similar to those in the linear isotropic hardening example,
except for the smooth yielding.

Loadcase 6:

In the 6th loadcase, the expansion of the elastic region due to isotropic
hardening is not that strong in previous exercises, since the stiffness K has
been significantly reduced.

61
62
3. Three Dimensional Models
In the one dimensional case, we could only speak about one direction of loading,
stresses and strains. There was no need to count in the different direction of the
variables, therefore we assumed all the variables to be scalars. In the 3D case, we
need to consider also vectors and tensors. The stress tensor is 2nd order (two
dimensional), consists of 9 components and is symmetric, therefore we speak about
6 independent variables:

𝜎𝜎11 𝜎𝜎12 𝜎𝜎13

𝜎𝜎
𝝈𝝈 = � 21 𝜎𝜎22 𝜎𝜎23 �, where: 𝜎𝜎12 = 𝜎𝜎21 , 𝜎𝜎13 = 𝜎𝜎31 , 𝜎𝜎23 = 𝜎𝜎32
𝜎𝜎31 𝜎𝜎32 𝜎𝜎33

A tensor or a vector will be written in bold. Whereas components of vectors and

tensors are not and have one (vector) or two (2nd order tensor) subscripts.

Displacement is defined as a vector:

𝑢𝑢1
𝒖𝒖 = �𝑢𝑢2 �
𝑢𝑢3

And the strain tensor is a symmetric gradient of the displacement field. The gradient
of a vector is a 2nd order tensor where all combinations of partial derivatives are
considered:

𝑢𝑢1,1 𝑢𝑢1,2 𝑢𝑢1,3

𝑢𝑢
∇𝒖𝒖 = � 2,1 𝑢𝑢2,2 𝑢𝑢2,3 �
𝑢𝑢3,1 𝑢𝑢3,2 𝑢𝑢3,3

𝜕𝜕𝑢𝑢1
Where 𝑢𝑢1,2 = 𝜕𝜕𝑥𝑥2
, is a derivative of the first component of the displacement vector

with respect to a second coordinate.

𝜀𝜀11 𝜀𝜀12 𝜀𝜀13

1
𝜺𝜺 = (∇𝒖𝒖 + ∇𝑇𝑇 𝒖𝒖) = �𝜀𝜀21 𝜀𝜀22 𝜀𝜀23 �
2 𝜀𝜀 𝜀𝜀32 𝜀𝜀33
31

In the algorithms (and in the continuum mechanics in general) 2nd order tensors like
stress or strain are often decomposed in the volumetric and deviatoric part. When we

63
consider a strain tensor of a cubic finite element, the volumetric part describes the
change in volume without shape change and the deviatoric part describes the change
in shape preserving the volume. Consider an arbitrary 2nd order tensor:

1
𝑨𝑨 = 𝑨𝑨𝑣𝑣𝑣𝑣𝑣𝑣 + 𝑨𝑨𝑑𝑑𝑑𝑑𝑑𝑑 , where 𝑨𝑨𝑣𝑣𝑣𝑣𝑣𝑣 = 𝑡𝑡𝑡𝑡(𝑨𝑨)𝑰𝑰 and 𝑨𝑨𝑑𝑑𝑑𝑑𝑑𝑑 = 𝑨𝑨 − 𝑨𝑨𝑣𝑣𝑣𝑣𝑣𝑣
3

The trace operator tr() returns a sum of diagonal elements of a tensor (scalar):

𝑡𝑡𝑡𝑡(𝑨𝑨) = 𝐴𝐴11 + 𝐴𝐴22 + 𝐴𝐴33

And I is a unit tensor:

1 0 0
𝑰𝑰 = �0 1 0�
0 0 1

As you may notice a volumetric part of a tensor returns a diagonal matrix with equal
components, we can simply imagine it as a scaling operator.

Besides that we introduce the following operators:

• Double contraction:

3 3

𝑨𝑨 ∶ 𝑩𝑩 = � � 𝐴𝐴𝑖𝑖𝑖𝑖 𝐵𝐵𝑖𝑖𝑖𝑖 = 𝐴𝐴𝑖𝑖𝑖𝑖 𝐵𝐵𝑖𝑖𝑖𝑖

𝑖𝑖=1 𝑗𝑗=1

In the last step we skipped the sum symbols, this is due to the Einstein summation
convention – when the same index appears twice, we sum up the operation for all
indices from 1 to 3.

Example for two vectors:

𝑎𝑎𝑖𝑖 𝑏𝑏𝑖𝑖 = 𝑎𝑎1 𝑏𝑏1 + 𝑎𝑎2 𝑏𝑏2 + 𝑎𝑎3 𝑏𝑏3

64
3.1. Elastoplasticity with Linear Isotropic and
Kinematic Hardening
Let’s begin with the elastic part of the strain energy density function, which is here
defined by means of Lame constants:

𝜆𝜆
𝜓𝜓 𝑒𝑒 (𝜺𝜺) = (𝜺𝜺 ∶ 𝑰𝑰)2 + 𝜇𝜇(𝜺𝜺2 ∶ 𝑰𝑰)
2

Where:
λ – First Lame constant
μ – Second Lame constant

The strain energy density function can also be expressed by means of volumetric and
deviatoric components of strain:

1
𝜓𝜓 𝑒𝑒 (𝜺𝜺) = 𝜅𝜅(𝜺𝜺𝑣𝑣𝑣𝑣𝑣𝑣 ∶ 𝑰𝑰)2 + 𝜇𝜇(𝜺𝜺𝑑𝑑𝑑𝑑𝑑𝑑 2 ∶ 𝑰𝑰)
2

Where:

2
κ – Bulk modulus: 𝜅𝜅 = 𝜆𝜆 + 𝜇𝜇
3

Extension to isotropic and kinematic hardening needs an introduction of plastic

component:

𝜓𝜓(𝜺𝜺) = 𝜓𝜓 𝑒𝑒 (𝜺𝜺) + 𝜓𝜓 𝑝𝑝 (𝜺𝜺)

Where:

1 1
𝜓𝜓 𝑝𝑝 (𝜺𝜺) = 𝐾𝐾𝑟𝑟𝛼𝛼 2 + 𝐾𝐾(1 − 𝑟𝑟)‖𝜷𝜷‖2
2 2

The isotropic hardening strain remains a scalar value, as it indicates only an expansion
of elastic region, the β strain is a tensor, since it explicitly defines a displacement of
elastic region in 3D space. Therefore, we need to consider a norm of the tensor to

65
compute a scalar quantity representing the tensor. The norm here is defined as
Euclidean norm (or p=2 norm).

The hardening stress 𝑅𝑅 (scalar) and Back stress 𝑩𝑩 (which in 3D case is a tensor) are
analogous to the 1D case obtained as partial derivatives of strain energy density
function:

𝜕𝜕𝜓𝜓 𝜕𝜕𝜕𝜕
𝑅𝑅 = − = −𝐾𝐾𝐾𝐾𝛼𝛼, 𝜝𝜝 = − = −𝐾𝐾(1 − 𝑟𝑟)𝜷𝜷
𝜕𝜕𝜕𝜕 𝜕𝜕𝜷𝜷

Yield stress scalar is not anymore, a sustainable criterion for plasticity in 3D, since the
stress state is described by 6 independent components. To define a limit criterion for
elastic loading, we can express a stress state by means of stress invariants, so called
principal stresses 𝜎𝜎1 , 𝜎𝜎2 , 𝜎𝜎3 :

𝜎𝜎1 0 0
𝝈𝝈 = � 0 𝜎𝜎2 0�
0 0 𝜎𝜎3

This form is obtained by choosing such a coordinate system, that shear stresses
disappear, and we only consider normal stresses. Such a form is obtained with help of
eigenvalue and eigenvector calculation, where eigenvalues are the magnitudes of
principal stresses and eigenvectors the directions of the new coordinate system.

3.1.1. Von Mises Yield Criterion

Stress tensor expressed by means of principal directions (only 3 components) enables

a 3D visualisation of a stress state. There are many yield stress criteria for 3D, however
we will only concentrate on the VonMises yield criterion. The VonMises criterion says,
that the yielding only depends on the deviatoric part of stress tensor. Volumetric
stress does not have an influence. The elastic region space can be visualised by means
of a cylinder:

66
Figure 3.1:Stress space for Von Mises yield criterion.

When the stress tensor “reaches” the boundaries of a cylinder, plastic deformation
occurs. On the figure you can see a green vector representing a certain stress state. It
can be decomposed into volumetric and deviatoric parts. The stress can be
represented as a vector, since we have only three components (principal stresses).
1
The hydrostatic axis has a direction along the vector: 𝑛𝑛 = [1; 1; 1]𝑇𝑇 . The hydrostatic
√3

axis defines the stress state where all the principal stresses are equal 𝜎𝜎1 = 𝜎𝜎2 = 𝜎𝜎3 .
The VonMises yield criterion is applied to develop the yield function:

2
𝛷𝛷(𝝈𝝈) = � �𝝈𝝈𝑑𝑑𝑑𝑑𝑑𝑑 � − 𝜎𝜎𝑦𝑦 ≤ 0
3

Where:

𝝈𝝈𝑑𝑑𝑑𝑑𝑑𝑑 – deviatoric stress tensor – tensors and vectors are marked as bolded,

𝜎𝜎𝑦𝑦 – scalar yield stress value.

Where 𝜎𝜎𝑦𝑦 is known from 1D as yield stress and remains a scalar value. Considering
the hardening components R and B we obtain the following form of yield function:

67
3
𝛷𝛷(𝝈𝝈, 𝑅𝑅, 𝑩𝑩) = �𝝈𝝈𝑑𝑑𝑑𝑑𝑑𝑑 + 𝑩𝑩� − � (𝜎𝜎𝑦𝑦 − 𝑅𝑅) ≤ 0
2

Where:

𝑩𝑩 – Back stress tensor,

𝑅𝑅 – hardening stress scalar.

Hardening can be visualised by observing VonMises stress cylinder in the deviatoric

1
plane (plane with a normal vector: 𝑛𝑛 = [1; 1; 1]𝑇𝑇 ):
√3

Figure 3.2: Hardening laws in 3D.

From the figure above, we can deduce that hardening stress R is only a scalar, because
it defines the increase in the cylinder radius, whereas Back Stress B must be a tensor
to fully define coordinates of a new centre of elastic region.

The dissipation inequality for 3D case is a scalar quantity, therefore stress tensor
needs to be double contracted with the strain rate:

𝐷𝐷 = 𝝈𝝈 ∶ 𝜺𝜺̇ − 𝜓𝜓̇ ≥ 0

After some manipulations and elimination of elastic components a reduced

dissipation inequality is obtained:

68
𝐷𝐷 𝑟𝑟𝑟𝑟𝑟𝑟 = 𝝈𝝈 ∶ 𝜺𝜺̇ 𝑝𝑝 + 𝑅𝑅𝛼𝛼̇ + 𝑩𝑩 ∶ 𝜷𝜷̇ ≥ 0

3.1.2. The Algorithm

As in one dimensional case, the basis for the algorithm will be the evolution equations.
Their form is also derived as a product of plastic multiplier and a partial derivative of
the yield function:

𝜕𝜕𝜕𝜕 𝝈𝝈𝑑𝑑𝑑𝑑𝑑𝑑 + 𝑩𝑩
𝜺𝜺̇ 𝑝𝑝 = 𝜆𝜆̇ = 𝜆𝜆̇ 𝑑𝑑𝑑𝑑𝑑𝑑
𝜕𝜕𝝈𝝈 ‖𝝈𝝈 + 𝑩𝑩‖

Where the second component of the product is in fact a unit (normalized) direction
of a deviatoric stress relative to the Back stress correction:

𝝈𝝈𝑑𝑑𝑑𝑑𝑑𝑑 + 𝑩𝑩
= 𝒏𝒏
‖𝝈𝝈𝑑𝑑𝑑𝑑𝑑𝑑 + 𝑩𝑩‖

The evolution equation for the hardening strain variables 𝛼𝛼 and 𝜷𝜷 are as follows:

2
𝛼𝛼̇ = � 𝜆𝜆̇, 𝜷𝜷̇ = 𝜆𝜆̇𝒏𝒏
3

By applying the Euler Backward algorithm to the strain variables, we obtain the basic
formulae in time discretization:

𝑝𝑝 𝑝𝑝 𝑝𝑝 𝑝𝑝
𝜺𝜺𝑛𝑛+1 = 𝜺𝜺𝑛𝑛 + ∆𝑡𝑡𝜺𝜺̇ 𝑛𝑛+1 = 𝜺𝜺𝑛𝑛 + ∆𝜆𝜆𝑛𝑛+1 𝒏𝒏𝑛𝑛+1

2
𝛼𝛼𝑛𝑛+1 = 𝛼𝛼𝑛𝑛 + ∆𝑡𝑡𝛼𝛼̇ 𝑛𝑛+1 = 𝛼𝛼𝑛𝑛 + ∆𝜆𝜆𝑛𝑛+1 �
3

𝜷𝜷𝑛𝑛+1 = 𝜷𝜷𝑛𝑛 + ∆𝑡𝑡𝜷𝜷̇𝑛𝑛+1 = 𝜷𝜷𝑛𝑛 + ∆𝜆𝜆𝑛𝑛+1 𝒏𝒏𝑛𝑛+1

By following similar procedures as in 1D case to develop the trial deviatoric stress and
trial yield function we obtain final algorithms for linear isotropic hardening:

%Known:

𝑝𝑝
𝜺𝜺𝑛𝑛+1 , 𝜺𝜺𝑛𝑛 , 𝛼𝛼𝑛𝑛

69
%Elastic prediction – compute trial values

𝑣𝑣𝑣𝑣𝑣𝑣
𝝈𝝈𝑛𝑛+1 = 𝜅𝜅 ∗ 𝑡𝑡𝑡𝑡(𝜺𝜺𝒏𝒏+𝟏𝟏 ) ∗ 𝑰𝑰

𝒑𝒑
𝝈𝝈𝒅𝒅𝒅𝒅𝒅𝒅,𝒕𝒕 𝒅𝒅𝒅𝒅𝒅𝒅
𝒏𝒏+𝟏𝟏 = 2𝜇𝜇(𝜺𝜺𝒏𝒏+𝟏𝟏 − 𝜺𝜺𝒏𝒏 )

𝑡𝑡
𝑅𝑅𝑛𝑛+1 = −𝐾𝐾𝛼𝛼𝑛𝑛

2
𝑡𝑡
𝛷𝛷𝑛𝑛+1 = �𝝈𝝈𝒅𝒅𝒅𝒅𝒅𝒅,𝒕𝒕
𝒏𝒏+𝟏𝟏 � −
𝑡𝑡
� (𝜎𝜎𝑦𝑦 − 𝑅𝑅𝑛𝑛+1 )
3

%Check if plastic deformation occurs

𝑡𝑡
if 𝛷𝛷𝑛𝑛+1 ≤0

𝛼𝛼𝑛𝑛+1 = 𝛼𝛼𝑛𝑛

𝑝𝑝 𝑝𝑝
𝜺𝜺𝑛𝑛+1 = 𝜺𝜺𝑛𝑛

𝒅𝒅𝒅𝒅𝒅𝒅,𝒕𝒕
𝝈𝝈𝒅𝒅𝒅𝒅𝒅𝒅
𝒏𝒏+𝟏𝟏 = 𝝈𝝈𝒏𝒏+𝟏𝟏

𝑡𝑡
𝑅𝑅𝑛𝑛+1 = 𝑅𝑅𝑛𝑛+1

else

𝒏𝒏𝑛𝑛+1 = 𝝈𝝈𝒅𝒅𝒅𝒅𝒅𝒅,𝒕𝒕 𝒅𝒅𝒅𝒅𝒅𝒅,𝒕𝒕

𝒏𝒏+𝟏𝟏 /�𝝈𝝈𝒏𝒏+𝟏𝟏 �

𝑡𝑡 2
∆𝜆𝜆𝑛𝑛+1 = 𝛷𝛷𝑛𝑛+1 /(2𝜇𝜇 + ∗ 𝛫𝛫)
3

2
𝛼𝛼𝑛𝑛+1 = 𝛼𝛼𝑛𝑛 + ∆𝜆𝜆𝑛𝑛+1 �
3

𝑝𝑝 𝑝𝑝
𝜺𝜺𝑛𝑛+1 = 𝜺𝜺𝑛𝑛 + ∆𝜆𝜆𝑛𝑛+1 𝒏𝒏𝑛𝑛+1

𝑑𝑑𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑𝑑𝑑,𝑡𝑡
𝝈𝝈𝑛𝑛+1 = 𝝈𝝈𝑛𝑛+1 − 2𝜇𝜇∆𝜆𝜆𝑛𝑛+1 𝒏𝒏𝑛𝑛+1

𝑣𝑣𝑣𝑣𝑣𝑣 𝑑𝑑𝑑𝑑𝑑𝑑
𝜎𝜎𝑛𝑛+1 = 𝜎𝜎𝑛𝑛+1 + 𝜎𝜎𝑛𝑛+1

𝑅𝑅𝑛𝑛+1 = −𝐾𝐾𝐾𝐾𝑛𝑛+1

And for linear isotropic and kinematic hardening:

%Known:

𝑝𝑝
𝜺𝜺𝑛𝑛+1 , 𝜺𝜺𝑛𝑛 , 𝛼𝛼𝑛𝑛 , 𝜷𝜷𝒏𝒏

70
%Elastic prediction – compute trial values

𝑣𝑣𝑣𝑣𝑣𝑣
𝝈𝝈𝑛𝑛+1 = 𝜅𝜅 ∗ 𝑡𝑡𝑡𝑡(𝜺𝜺𝒏𝒏+𝟏𝟏 ) ∗ 𝑰𝑰

𝒑𝒑
𝝈𝝈𝒅𝒅𝒅𝒅𝒅𝒅,𝒕𝒕 𝒅𝒅𝒅𝒅𝒅𝒅
𝒏𝒏+𝟏𝟏 = 2𝜇𝜇(𝜺𝜺𝒏𝒏+𝟏𝟏 − 𝜺𝜺𝒏𝒏 )

𝑡𝑡
𝑅𝑅𝑛𝑛+1 = −𝑟𝑟𝑟𝑟𝛼𝛼𝑛𝑛

𝑡𝑡
𝑩𝑩𝑛𝑛+1 = −(1 − 𝑟𝑟)𝐾𝐾𝜷𝜷𝑛𝑛

𝝃𝝃𝑡𝑡𝑛𝑛+1 = 𝝈𝝈𝑑𝑑𝑑𝑑𝑑𝑑,𝑡𝑡 𝑡𝑡
𝑛𝑛+1 + 𝑩𝑩𝑛𝑛+1

2
𝑡𝑡
𝛷𝛷𝑛𝑛+1 = ‖𝝃𝝃𝒕𝒕𝒏𝒏+𝟏𝟏 ‖ − � (𝜎𝜎𝑦𝑦 − 𝑅𝑅𝑛𝑛+1
𝑡𝑡
)
3

%Check if plastic deformation occurs

𝑡𝑡
if 𝛷𝛷𝑛𝑛+1 ≤0

𝛼𝛼𝑛𝑛+1 = 𝛼𝛼𝑛𝑛

𝑝𝑝 𝑝𝑝
𝜺𝜺𝑛𝑛+1 = 𝜺𝜺𝑛𝑛

𝒅𝒅𝒅𝒅𝒅𝒅,𝒕𝒕
𝝈𝝈𝒅𝒅𝒅𝒅𝒅𝒅
𝒏𝒏+𝟏𝟏 = 𝝈𝝈𝒏𝒏+𝟏𝟏

𝑡𝑡
𝑅𝑅𝑛𝑛+1 = 𝑅𝑅𝑛𝑛+1

𝜷𝜷𝑛𝑛+1 = 𝜷𝜷𝑛𝑛

𝑡𝑡
𝜝𝜝𝑛𝑛+1 = 𝑩𝑩𝑛𝑛+1

else

𝒏𝒏𝑛𝑛+1 = 𝝃𝝃𝒕𝒕𝒏𝒏+𝟏𝟏 /‖𝝃𝝃𝒕𝒕𝒏𝒏+𝟏𝟏 ‖

𝑡𝑡 2
∆𝜆𝜆𝑛𝑛+1 = 𝛷𝛷𝑛𝑛+1 /(2𝜇𝜇 + ∗ 𝛫𝛫)
3

2
𝛼𝛼𝑛𝑛+1 = 𝛼𝛼𝑛𝑛 + ∆𝜆𝜆𝑛𝑛+1 �
3

𝜷𝜷𝑛𝑛+1 = 𝜷𝜷𝑛𝑛 + ∆𝜆𝜆𝑛𝑛+1 𝒏𝒏𝑛𝑛+1

𝑝𝑝 𝑝𝑝
𝜺𝜺𝑛𝑛+1 = 𝜺𝜺𝑛𝑛 + ∆𝜆𝜆𝑛𝑛+1 𝒏𝒏𝑛𝑛+1

𝑑𝑑𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑𝑑𝑑,𝑡𝑡
𝝈𝝈𝑛𝑛+1 = 𝝈𝝈𝑛𝑛+1 − 2𝜇𝜇∆𝜆𝜆𝑛𝑛+1 𝒏𝒏𝑛𝑛+1

𝑣𝑣𝑣𝑣𝑣𝑣 𝑑𝑑𝑑𝑑𝑑𝑑
𝜎𝜎𝑛𝑛+1 = 𝜎𝜎𝑛𝑛+1 + 𝜎𝜎𝑛𝑛+1

71
𝑅𝑅𝑛𝑛+1 = −𝑟𝑟𝐾𝐾𝐾𝐾𝑛𝑛+1

𝜝𝜝𝑛𝑛+1 = −(1 − 𝑟𝑟)𝐾𝐾𝜷𝜷𝑛𝑛+1

3.1.3. Exercise 6: 3D Elastoplasticity with Linear Isotropic

Hardening

In this exercise we need to introduce tensor values as 3x3 matrices. The strain
history array will actually be an array of matrices, hence, a 3D array, which
can be visualized as follows:

Figure 3.3: 3D Array for Tensors visualization.

Let’s then start with defining a new loadcase generating function loadcase3D.oml

1. Open the MaterialModelling.aprj project file if it is not yet open. Create a new
file named loadcase3D.oml and add it to the project. Function definition
should look as follows:

2. Then we initialize the 3D array of strains in the following manner:

The first and second dimension is corresponding to tensor dimesions, the
third dimension is the direction of the array.

72
3. Now it’s time to define 5 loadcases. First 3 will contain only normal
components (no shear strain), they will depict a “linear in time” behaviour.
The fourth one is a simple shear strain and last one sinusoidally pulsating in
normal components.

a. Pure volumetric expansion

Only the diagonal components are updated and in each timestep they
have equal values. There is no deviatoric strain, in case a decomposition
into volumetric and deviatoric would be applied. There is a linear
evolution from 0 to 0,2% of strain.

b. Tension with different normal components

In this case the 2nd normal component evolves 3 times faster than the
rest. That means deviatoric behaviour is already present. There is no
shear strain.

c. Tension and compression with different normal components

73
As you can notice, the strain values for 1st and 3rd direction have been
reversed to cause compression. In that case there is more deviatoric
strain observed, however still no shear strain is applied.

d. Simple shear

For simple shear, there is only one shear strain variable subjected to
change. Since strain tensor is symmetric, the corresponding two
components need to be equal all over the time evolution. In that case it
is a shear in 1,2-plane. Therefore 12 and 21 components are updated.

e. Cyclic sinusoidal tension/compression with different cycle time.

Sinusoidal pulsation applied to each volumetric component with varying

amplitudes (0,1%, 0,2% and 0,15%) and varying cycle times (4,8 and 6
used as multiplicators accordingly).

4. Save your loadcase3D.oml file.

74
5. Now it is time to create a code for elastoplasticity with isotropic hardening.
Create a new file named: EPIsoHardening3D.oml, save it and add to our
project

6. Introduce time discretization parameters as usual.

7. When it comes to material parameters, we additionally need to define

Poisson’s ratio ν. Given both E and ν, full 3D isotropic material stiffness is
defined.
They can be converted to obtain Kirchhoff modulus, Lame constants as well
as bulk modulus. For our code, we need to obtain bulk modulus and second
Lame constant, as defined in our strain energy density function.

8. Convert the E and ν to bulk modulus κ and second Lame constant μ:

75
9. In the tensor calculus, we need quite frequently to use the unit tensor.
Fortunately to obtain it, we simply need to call the eye() function. To get 3x3
unit tensor, we need to input 3 as an argument.

10. Now it is time to select a loadcase. We let the user choose it during run time:

Given a loadcase number, we can call the loadcase3D function and provide a
loadcase number together with number of timesteps.

11. Before the algorithm starts, we need to initialize all the variables, that need
to occur in form of arrays:

The plastic strain array, the volumetric and deviatoric components of the
input array (strain array), the stress array together with its volumetric and

76
deviatoric components are all initialized as a tensor array. The hardening
stress and hardening strain variable alpha still remain scalar, since they only
count for expansion of the VonMises “circle”.

12. Additionally, we can in advantage compute the volumetric and deviatoric

strains, so that we clear the algorithm a little bit from this calculation.

We use the trace() function to obtain the trace of the strain tensor at the
position i in the array. Reminder: Id is a unit tensor.

13. Now it is time to provide the algorithm

We use the norm() function to obtain an Euclidean norm of a tensor.

77
For reference, check for the algorithm structure provided before the exercise
chapter. The idea is the same as in the 1D case, except for tensor update
instead of scalar update for stresses and alpha. n_dir variable is the unit
direction of deviatoric stress. The stress tensor is updated as a sum of
deviatoric and volumetric components, therefore can be updated in the
common scope outside if-else statements.

14. Next step is to define proper plots. Since now we need to investigate single
components of tensors, the ideal way would be to use 3x3 subplots for each
tensor and plot them within different figures. Plotting all the components of
the tensors would require a lot of code normally, therefore it will be easier to
use for loops to traverse through tensor components.

Figure 1.1 wil depict the strain evolution for each of its components. We use
two for loops to traverse in row and column direction. Each subplot is
activated as a combination of i,j indices. Since the eps data is a 3D array, we
cannot use it as an argument for the plot() function. Instead, we need to

78
create temporary arrays, which will store time history of a current strain
tensor component, indicated by i,j indices. Therefore, we need to insert a
third for loop to traverse through all time steps and assign the ij component
value. We plot both total and plastic strain on the same plots (see figure
below):

Next plot will contain results for total, volumetric and deviatoric stress
tensors. We follow the same procedure as for plotting strain data.
Last figure will provide us with alpha and hardening stress R results:

79
15. Save your file. Now we are ready to run the script. In the command window
you will be asked for the loadcase number. Try the script for all the loadcases.

80
The figures may appear too small on laptop monitors, therefore it is more
convenient to detach the plot window and expand it freely:

Loadcase 1:

81
The strain increases linearly, there is no plastic strain, since the strain is purely
volumetric.
The total and volumetric stress is equal, no deviatoric stress occurred. Shear
stress is equal to 0.

Loadcase 2:

82
Since in loadcase 2, the normal strain components are different in values,
the deviatoric part of strains is nonzero this time. You can see that plastic
strain for the component 11 and 33 is negative, since the tension in 22
direction is much greater than those in 11 and 33 causing the deviatoric
strains to obtain negative values in these directions.
We can check how it happens by doing a simple calculation. Let’s assume a
strain state to be as follows:
1 0 0
[𝜀𝜀𝑖𝑖𝑖𝑖 ] = �0 3 0�
0 0 1
Then we compute the volumetric strain using the formulae given before:
1
𝜺𝜺𝑣𝑣𝑣𝑣𝑣𝑣 = 𝑡𝑡𝑡𝑡(𝜺𝜺)𝐼𝐼
3
Which gives:
5/3 0 0
𝑣𝑣𝑣𝑣𝑣𝑣
[𝜀𝜀𝑖𝑖𝑖𝑖 ]=� 0 5/3 0 �
0 0 5/3
The deviatoric part is computed as a difference between the total and
volumetric part:

83
5 2
⎡1 − 0 0 ⎤ ⎡− 0 0 ⎤
⎢ 3 ⎥ ⎢ 3 ⎥
5 4
𝑑𝑑𝑑𝑑𝑑𝑑
�𝜀𝜀𝑖𝑖𝑖𝑖 �=⎢ 0 3− 0 ⎥=⎢ 0 0 ⎥
⎢ 3 ⎥ ⎢ 3 ⎥
⎢ 5⎥ ⎢ 2⎥
⎣ 0 0 1− ⎦ ⎣ 0 0 − ⎦
3 3
Which proves, that although the body was subjected to only tensile strains,
the deviatoric strain tensor results in negative components due to magnitude
difference in the total strain.

84
The stress figures look as follows,

here you can also notice that the deviatoric part of the stresses is in the
negative range for the normal stress components except for the direction 22.
The total and volumetric stress has similar values, but for the 11 and 33 cases
the volumetric stress is even greater than the total.

This time hardening also occurred:

We can see that the elastic region has explanded by 12MPa in this loadcase.

85
Loadcase 3 looks similar to the 2nd one, except that a different scale is
observed:

The stress values:

86
And the hardening variables:

Loadcase 4:

This time we consider shear stresses. As you can see, the plastic deformation
also occurs in the same manner as for the previous cases.

87
Important thing to notice here is that there is no volumetric stress in the case
of simple shear. The volume is preserved and only shape deformation is
observed, therefore the total and deviatoric stress coincides in the figure
above for 12 (and 21) direction.

Hardening variables:

88
Loadcase 5:

The strain history for each normal component varies in both amplitude and
cycle time, therefore the combined deformation is quite complicated, which
results in very irregular plastic strain evolution.

89
Stress values:

For the stress the results happen to be even more complicated. For counting
in the plastic behaviour we are mostly interrested in deviatoric stress, when
the values of deviatoric stress exceed a certain limit, plastic strain is observed.

Hardening variables:

Each occurrence of plastic strain results in an increase of alpha and hardening

stress, which are of course cummulative parameters.

90
3.1.4. Exercise 7: 3D Elastoplasticity with Linear Isotropic and
Kinematic Hardening

For this exercise, a successful completion of Exercise 6 is necessary, as this

exercise is in fact an extension to the previously written code. The aim is to include
the kinematic part of hardening into the exercise.

1. Create a new file and copy-paste the content of the EPIsoHardening3D.oml

file. Save it as EPIsoKinHardening3D.oml and add to the MaterialModelling

project.
2. We should live the part defining the discretization parameters and start
extending the code from requesting the percentage of isotropic hardening in
the model within I/O. For that, include a line of following code just before the
line for reading loadcase number:

Save the value in the r variable.

91
3. Next extension will be to add initialization of β and B variables:

Since these are tensors, we need to initialize then an array of 3x3 matrices.
Now it is time to move into the algorithm. As you may remember from the 1D
isotropic-kinematic model, the variables α and β are dependent on the r
percentage of isotropic hardening, therefore this variable needs to be added
to the lines 54 and 55. We add also the ξ variable to count in the shift in elastic
region center coming from kinematic hardening. Now the ξ variable is used in
the trial yield function and later in defining the direction of radial return.

4. Expand the elastic update part for kinematic hardening variables:

92
5. A few changes have to be made for Radial Return part. Since now we use the
correction direction n_dir by the Back stress value:

6. Instead of adding more plots for kinematic hardening variables, we are more
interested in seeing how the elastic region has changed during simulation.
Since it only depends on deviatoric stress, we will create a 2D plot in the
1
deviatoric plane, that is a plane with a normal vector of: 𝒏𝒏 = [1; 1; 1]𝑇𝑇 . In
√3

order to visualise it, we need to use hardening stress R to count in expansion

of a circle and additional Back stress B to model displacement of the center.
While using hardening stress is trivial, since this is a scalar, visualization of the
Back stress is more complicated. The procedure is to extract the Eigenvalues
of B so that we obtain the values for principal directions, only then can we
transfer the tensor to 3D plot. However in our case we create a 2D plot with
a given orientation before. Therefore we need to project the 3 components
into a 2D system. All the operations are done within visualize.oml function,
which you need to create, copy the following code:

function [a] = visualize(sig_y,R,B)

b = eig(B);

x = 1/sqrt(3) * (b(1)*cos(pi/6) - b(2)*cos(pi/6));

y = 1/sqrt(3) * (b(1)*cos(pi/3) + b(2)*cos(pi/3) - b(3));

f4 = figure(4);

th = 0:pi/50:2*pi;

93
x1 = sig_y * cos(th);

y1 = sig_y * sin(th);

plot(x1,y1);

title('Deviatoric stress plane');

hold on;

x2 = (-R+sig_y) * cos(th) + x;

y2 = (-R+sig_y) * sin(th) + y;

plot(x2, y2);

legend('Initial','Deformed');

hold off;

end

7. Remember to save the function in the same directory and add it to the
project. Now we need to call the function at the end of the
EPIsoKinHardening3D.oml file:

94
We need to extract a 2D matrix out of 3D Back stress array for the last iteration.
The same for hardening stress – only last iteration value is passed. In the end we
pass 2 scalars: yield and hardening stress and one tensor – Back stress. Save all
the files and execute it.

8. Let’s try to compare different hardening behaviours for loadcase 3:

100% isotropic hardening:

95
100% kinematic hardening:

96
50 % isotropic and 50% kinematic hardening:

97
4. Finite Element Method
Application
In this chapter the aim will be to implement elastoplasticity in a simple one
dimensional FEM model. The model will contain the simplest of the 1D elements –
trusses, which support only axial loading. Thanks to that, it would be relatively easy
to implement the solution for 1D elastoplasticity with Isotropic hardening presented
in the previous chapters.

4.1. A Brief Overview of FEM Theory

Please note that the following material is assuming, that the user is familiar with
theoretical background of Finite Element Method and hence, only an overview of
basic FEM theory will be presented. Since a simple FEM example for two truss system
is already discussed in Compose E-Book (https://altairuniversity.com/free-ebook-
math-scripting-compose/), we will mainly focus on presenting theoretical basics for
nonlinear FEM analysis.

But first let’s collect basic equations necessary to build linear FEM analysis. In the
Compose eBook following equation was derived by means of local and global
formulation:

𝑲𝑲𝑲𝑲 = 𝒇𝒇

K – global stiffness matrix,

u – global displacement vector,
f – global external force vector.

The global stiffness matrix K is built up out of local, elemental stiffness matrices. The
dimension of stiffness matrix is equal to total number of nodes in the model. When

98
two elements share one node, the corresponding entry in global stiffness matrix is
combining two entries from elemental stiffness matrices.
For 1D Truss element:
𝐷𝐷𝐷𝐷 1 1
𝐾𝐾𝑒𝑒𝑒𝑒 = � �
𝑙𝑙 −1 1
D – tangent modulus,
A – cross section of an element,
l – length of an element.

For linear elastoplasticity with isotropic hardening:

Elastic state:
𝐷𝐷 = 𝐸𝐸
Elastoplastic state:
𝐸𝐸𝐸𝐸
𝐷𝐷 =
𝐸𝐸 + 𝐻𝐻
E – Young’s modulus,
H – Hardening stiffness.

Figure 4.1 : Elastoplasticity with hardening

99
The stiffness matrix ”K”is defined as:

𝑲𝑲 = � 𝑩𝑩𝑻𝑻 𝑫𝑫 𝑩𝑩 𝑑𝑑𝛺𝛺
𝛺𝛺

Where “B” matrix is used for mapping between elemental quantities and nodal
quantities. For one truss element (elemental “B” matrix) is defined as:
1 −1
𝑩𝑩𝒆𝒆𝒆𝒆 = � �
𝑙𝑙 1
For two truss elements we get a global “B” matrix:
−1/𝑙𝑙1 0
𝑩𝑩 = � 1/𝑙𝑙1 −1/𝑙𝑙2 �
0 1/𝑙𝑙2
Since in the nonlinear FEM we will be operating on load and displacement increments,
we use the following relation:
∆𝜺𝜺 = 𝑩𝑩Δ𝒖𝒖
And for stress-strain relation we assume the following approximation:
∆𝝈𝝈 = 𝑫𝑫∆𝜺𝜺
Nonlinear analysis will be conducted using incremental-iterative procedure. The
external load vector will be divided into a certain number of increments and applied
one by another. After each load increment, the equilibrium of external and internal
forces will be seeked using Newton-Raphson iterations.
The equilibrium is achieved when the residual approaches zero:
𝒓𝒓 = 𝒇𝒇 − 𝒒𝒒
The internal force vector “q” must be a nodal quantity, therefore the stresses are
mapped using the “B” matrix:

𝒒𝒒 = � 𝑩𝑩𝜯𝜯 𝝈𝝈 𝑑𝑑𝛺𝛺
𝛺𝛺

The complete algorithm looks as follows [4]:

User input: Material constants, Boundary conditions, geometry

Initialize: 𝒖𝒖 = 0, 𝒇𝒇 = 0, 𝝈𝝈 = 0

Compute: 𝑩𝑩, ∆𝒇𝒇

100
Load increments: for 𝑘𝑘 = 1: 𝑛𝑛_𝑖𝑖𝑖𝑖𝑖𝑖

𝑲𝑲 = 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠(… );

Split stiffness matrix, force, displacement: 𝑲𝑲𝒖𝒖𝒖𝒖 , 𝑲𝑲𝒖𝒖𝒖𝒖 , …

Displacement predictor: ∆𝒖𝒖 = 𝜥𝜥−𝟏𝟏 (∆𝒇𝒇𝒖𝒖 − 𝑲𝑲𝒖𝒖𝒖𝒖 𝒖𝒖𝒑𝒑 );

Reaction forces: 𝒇𝒇𝒓𝒓𝒓𝒓𝒓𝒓,𝒑𝒑 = 𝒇𝒇𝒓𝒓𝒓𝒓𝒓𝒓,𝒑𝒑 + 𝑲𝑲𝒑𝒑𝒑𝒑 ∆𝒖𝒖𝒑𝒑 + 𝑲𝑲𝒑𝒑𝒑𝒑 ∆𝒖𝒖𝒖𝒖 − ∆𝒇𝒇𝒑𝒑

Newton-Rhapson: for j = 1:max_iter

Strain increment: ∆𝜺𝜺 = 𝑩𝑩∆𝒖𝒖

Stress: 𝝈𝝈𝒏𝒏𝒏𝒏𝒏𝒏 = 𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟(𝝈𝝈𝒄𝒄𝒄𝒄𝒄𝒄𝒄𝒄𝒄𝒄𝒄𝒄𝒄𝒄 , ∆𝜺𝜺)

Internal force: 𝒒𝒒 = 𝑩𝑩𝑻𝑻 (𝝈𝝈.∗ 𝑨𝑨.∗ 𝒍𝒍)

Residual: 𝒓𝒓 = (𝒇𝒇 + ∆𝒇𝒇 + 𝒇𝒇𝒓𝒓𝒓𝒓𝒓𝒓 ) − 𝒒𝒒

if 𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛(𝒓𝒓) < 𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛(𝒇𝒇 + ∆𝒇𝒇 + 𝒇𝒇𝒓𝒓𝒓𝒓𝒓𝒓 )

Converged: break

else

𝑲𝑲𝒑𝒑𝒑𝒑 −𝟏𝟏
Update Δu: 𝜹𝜹𝒖𝒖𝒖𝒖 = � � 𝒓𝒓
𝑲𝑲𝒖𝒖𝒖𝒖

∆𝒖𝒖 = ∆𝒖𝒖 + 𝜹𝜹𝜹𝜹

end: if-else

end: Newton-Raphson

Update u: 𝒖𝒖 = 𝒖𝒖 + ∆𝒖𝒖

Update f: 𝒇𝒇 = 𝒇𝒇 + ∆𝒇𝒇

end: Load increments

Print: 𝒖𝒖, 𝝈𝝈

101
4.1.1. Exercise 8: Elastoplasticity in 1D Truss System

In the following exercise we will model a very simple two element truss system,
where one of the elements is meant to undergo plastic deformation. The system
will be encastred at the one end and will be subjected to two external loads: an

evenly distributed force “q” and a concentrated force “F” applied to the middle
node.
The values for the forces will be selected such that the left element will deform
plastically and the right one will remain in the elastic state. Since we stay with one
dimensional model, all the quantities involved in the calculations remain scalars.
Since the element type is truss, only axial loading and deformation is supported.
Create a new file, type in:

And save the file as TrussSystem.oml.

1. Model Setup
The first part of the code will contain all the neccesary data to create a model
setup and to choose the criteria for the simulation. In comparison to the typical
HyperWorks analysis workflow, we can refer to this step as preprocessing.

102
1.1. Material Parameters
We will start with providing the necessary material data. Since the material
model is elastoplastic (with isotropic hardening), we need three parameters:
Young`s modulus E, Hardening stiffness H, and Yield stress sig_y:

1.2. Geometry Approximation

As mentioned before, the geometry is approximated by two, of equal length
and of equal cross section truss elements. We will frequently use variables
storing number of elements and nodes, therefore we need to define them.

The dimensions of trusses are provided as column arrays:

1.3. Boundary Conditions

For fixed displacement, the strategy will be to create a boolean array

indicating which of the nodes are belonging to Dirichlet Boundary Conditions
(DBC, constraints on displacements). The reason for this is that later in the
code we will need to split the stiffness matrix and external load vector into
parts containing DBC nodes and the rest.
The external concentrated forces are provided as array, where each entry
coresponds to one node.

103
The external distributed force array however contains one entry per element.

1.4. Convergence Criteria

The last part of the user input is to set up control parameters for the
incremental-iterative algorithm. As introduced before, the external load will
be applied in the so-called increments. The number of increments chosen by
the user can be critical when it comes to convergence of the calculations.
Obviously the smaller the increments, the more likely the simulation will
converge. Therefore, we choose to divide our load in 100 increments.

Analogous situation is with the limit of Newton-Rhapson iterations for each

load increment. However, very big numbers of max. iterations are usually not
helpful when the problem with convergence appears to be caused by large
load increments. Properly set up models should converge within a small
number of iterations.

The tolerance value is responsible for the accuracy of the final results.
Besides errors due to geometry approximation, material model simplification
and BC approximation, the convergence tolerance is an additional factor that
may introduce unwanted error in the nonlinear analysis. Here, the tolerance
value is a relative quantity used to compare the norm of residual array
(vector) with the norm of external force array (vector).

104
2. Initialization of Variables
The next section is devoted to initialize all the necessary variables before the
actual algorithm is started. We want to compute as much of the quantities as
possible outside algorithm so that the simulation is more robust.

2.1. Global Variables

We start with initialization to zeroes of all the variables that need to remain
in the global scope of the code. These are: total displacement, total stress,
total external force, total reaction force, force increment (force increment is
constant) and the B matrix (constant). All these variable arrays are
cumulative (updated in each loop run with certain increment value) except
for the force increment and B matrix.

2.2. The “B” Matrix

Definition of B matrix is not necessary, but much helpful in keeping the code
short and clean. It is a matrix that maps strain data array into displacement
array according to the geometry information (since displacement is related
to strain by the length of elements in case of trusses).

Notice that the B matrix is not a square matrix. It has n_el number of rows
and n_nod number of columns, since strain is an elemental quantity and
displacement a nodal quantity.

105
2.3. Dirichlet Boundary Indices
As mentioned before, the information about which nodes belong to DBC and
which do not is relevant for division of stiffness matrix and force vector. That
operation is trivial for two element system however for more complex
systems a general algorithm might be useful. Therefore, for the sake of
practice, we will apply a general procedure:

We initialize the DBC and nDBC as empty arrays since we don’t know what their
final size will be. Next we use a for loop and append the indices to
correspoding arrays.

2.4. External Load Vector Setup

Since both the external concentrated load (Neumann boundary) as well as
distributed load remain constant and, most importantly, the load increment
is decided to be constant (the same) for each step, we can compute the total
force vector outside the load increment loop.
Our aim is to compute the total external force increment df, therefore we
first divide our arrays of Neumann force df_n and distributed force dq_o by
number of increments n_inc.

106
Since the external load vector (increment) df is a nodal array, we need to
convert the distributed load array dq_o into equivalent nodal forces. For that,

we compute the equivalent nodal forces for each element (df_el) and
append them to our global load increment (df) within a for loop. When we
are done with the distributed load, we simply add the external concentrated
forces (df_n).

3. The Algorithm
The main strategy for the simulation is to apply the force increments and calculate
the displacement increment, obtain an equilibrium in the new loading conditions
and repeat this procedure until the final force values are reached. We use a for
loop to go through all the load increments.

3.1. Current Stiffness matrix “K”

The stiffness matrix, as opposed to force vector, cannot be computed
globally before the algorithm, since the stiffness varies depending on
whether the elements are in plastic or elastic deformation state. Therefore
it needs to be recalculated in each load increment.
To keep the main code short and tidy, we create an external .oml function to
return the global stiffness matrix.

107
For the function, we need to pass the current stress state, geometry
information and material parameters. All these parameters allow to decide
what is the current stiffness for each element. The function is defined

separately, but located in the same directory, file named stiffnessMatrix.oml.

The tangent modulus D is a slope of the stress-strain curve. For plastic state,
the slope is defined as a combination of Young’s modulus and hardening
stiffness. The global stiffness matrix K is a square matrix of n_nod dimesion, it
is a setup of elemental stiffness matrices for all elements (K_el, 2x2 size - two
nodes for each element).

3.2. Stiffness Matrix Split

Due to presence of Dirichlet Boundary Conditions (DBC), the stiffness matrix
must be split to allow separate calculation of unknown displacements and
reaction forces. The total stiffness equation has a following form:
𝐾𝐾𝑝𝑝𝑝𝑝 𝐾𝐾𝑝𝑝𝑝𝑝 𝑢𝑢𝑝𝑝 𝑓𝑓𝑝𝑝 𝑓𝑓
� � �𝑢𝑢 � = � � + � 𝑟𝑟𝑟𝑟𝑟𝑟 �
𝐾𝐾𝑢𝑢𝑢𝑢 𝐾𝐾𝑢𝑢𝑢𝑢 𝑢𝑢 𝑓𝑓𝑢𝑢 0

108
3.3. Displacement Increment Predictor Step
The calculation of displacement increment predictor du is nothing else but
the assumption that the deformation is purely elastic. The previously
performed stiffness matrix division allows the extraction of equation for
unknown displacements:
𝐾𝐾𝑢𝑢𝑢𝑢 𝑢𝑢𝑢𝑢 = 𝑓𝑓𝑢𝑢 − 𝐾𝐾𝑢𝑢𝑢𝑢 𝑢𝑢𝑝𝑝
Where the last part of the RHS equation 𝐾𝐾𝑢𝑢𝑢𝑢 𝑢𝑢𝑝𝑝 = 0, when DBC
displacements are fixed (there can be enforced displacement also). For the
completeness, this term is included in the code. In our case we compute the
increment of displacement du:

3.4. Reaction Forces

When it comes to reaction forces, we need their total values for the residual
evaluation. Therefore the approach is to increment reaction forces in each
loop run, so that it balances the external loads accurately (sum of external
loads and reaction forces must give 0 in the end). Since the force increment
is constant and calculated before the algorithm, the reaction forces are
calculated by means of du obtained in the predictor step, this gives proper
values of f_rec. The equation for reaction forces is also extracted thanks to
stiffness matrix split:
𝑓𝑓𝑟𝑟𝑟𝑟𝑟𝑟 = 𝐾𝐾𝑝𝑝𝑝𝑝 𝑢𝑢𝑝𝑝 + 𝐾𝐾𝑝𝑝𝑝𝑝 𝑢𝑢𝑢𝑢 − 𝑓𝑓𝑝𝑝
This time the term 𝐾𝐾𝑝𝑝𝑝𝑝 𝑢𝑢𝑝𝑝 = 0, but it is still included in the code for the sake
of completeness.

Notice that all the terms updating the reaction forces are based on increments
and not total values of displacement and force.

109
3.5. Modified Newton-Rhapson Iterations
Given the displacement increment from the predictor step, we are ready to
perform Newton-Rhapson iterations to obtain the true displacement
increments, which will count in the eventual plasticity. The aim is to provide
a stress state that will balance out the external forces – external<->internal
force equilibrium. Since the stress state is based indirectly on displacements
(strains to be precise), the displacement obtained in the predictor step
should be updated until the residual will aproach zero:
𝑟𝑟 = 𝑓𝑓 𝑒𝑒𝑥𝑥𝑥𝑥 − 𝑓𝑓 𝑖𝑖𝑖𝑖𝑖𝑖 = (𝑓𝑓 + 𝑑𝑑𝑑𝑑 + 𝑓𝑓𝑟𝑟𝑟𝑟𝑟𝑟 ) − 𝑞𝑞 → 0
We need to create another for loop for iterations. First we evaluate the
strain increment deps based on current displacement increment prediction
du, later on we use the known from previous chapters radial return algorithm

to compute current stress state prediction.

Next we compute internal force vector according to:
2𝑥𝑥 𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 2𝑙𝑙
𝒒𝒒 = � 𝜝𝜝𝛵𝛵 𝝈𝝈 𝑑𝑑𝛺𝛺 �⎯⎯⎯⎯⎯� � 𝜝𝜝𝛵𝛵 𝝈𝝈 𝐴𝐴(𝑥𝑥) 𝑑𝑑𝑑𝑑
𝑑𝑑𝛺𝛺 0

For two truss elements numerical calculation of the given integral gives exact
results for the following formula:
𝑞𝑞1 𝐵𝐵
𝒒𝒒𝒆𝒆𝒆𝒆 = �𝑞𝑞 � = � 1 � (𝜎𝜎𝜎𝜎𝑙𝑙) = 𝑩𝑩𝑻𝑻𝒆𝒆𝒆𝒆 (𝜎𝜎𝑒𝑒𝑒𝑒 𝛢𝛢𝑙𝑙)
2 𝐵𝐵2

110
If the norm of the residual norm(r) is fulfilling the convergence criterion, then
the iterations are stopped. If however it is not fulfilled, the displacement
increment du is updated by update_du value computed using stiffness matrix
and residual vector.

Additionally, we want to make sure that our iterations always converge for
each load increment. Hence, in case diverging step occurs, we want to print
out a Warning message informing us in which load increment the modified
Newton-Rhapson iterations didn’t converge and what is the residual error.

3.6. Radial Return Function

For the stress state update we use the well-known radial return algorithm.
The only modification in comparison to previous exercises is that we do not
have the alpha strain keeping track of the hardening. Instead we obtain
hardening stress directly by subtracting yield stress from the current stress
state. For more explanation please refer to 1D Isotropic hardening chapter.

111
3.7. Variables Update
Outside the scope of Modified Newton-Rhapson iterations, but still in the
scope of incremental algorithm, we want to update the total force f and total
displacement u with the known force increment df and computed
displacement increment du.

At the end of the simulation, we want to print the final displacement array
and stress aray to the command window.

4. Evaluate the Code

Now it is time to run our program. With the data given in this exercise the results
should be as follows:

112
Don’t worry if your results are not exactly the same, but differ a bit. This is due to
round-off error and happens when the order of calculations is different. These
results have been compared to OptiStruct analysis. A 1D linear static analysis was
done in order to compare the stresses (in 1D there is no stress concentrations,
but on the other hand it does not support nonlinear analysis). The displacements
were compared with 3D cylinder analysis with Nonlinear quasi-static loadstep.

1D Linear Static load step to compare the stresses (the visible rod diameter is not
a true one, rather a random one, since it is not possible to visualise exact cross
sections in HyperView):

113
3D Nonlinear Quasi-Static analysis to compare the displacement (exact geometry
visible):

As you can see it is fair to say that the stress values were computed accurately.
Regarding the displacement, considering the fact that the analysis in OptiStruct
was 3D, it is understandable that the results will vary a little bit more. Still, the
difference is of the 1-2% range, hence we may say that the Compose code is
providing reasonable results.

114
115
5. Literature
[1] L. Yaw, Louie, “Nonlinear Static - 1D Plasticity - Various Forms of Isotropic
Hardening”, Walla Walla University, January 25, 2012.

[2] M. Kojic and K. J. Bathe. Inelastic Analysis of Solids and Structures. Springer-Verlag,

New York, 2004.

[3] https://en.wikipedia.org/wiki/Successive_over-
relaxation#Other_applications_of_the_method

[4] K. Krabbenhoft, “Basic Computational Plasticity”, Department of Civil Engineering,

Technical University of Denmark, June 2002.

116