Computational Solid Mechanics
Computational Plasticity
Chapter 2. 1D Plasticity Algorithms
C. Agelet de Saracibar
ETS Ingenieros de Caminos, Canales y Puertos, Universidad Politcnica de Catalua (UPC), Barcelona, Spain
International Center for Numerical Methods in Engineering (CIMNE), Barcelona, Spain
1D Plasticity Algorithms > Contents
Contents
Contents
1. Introduction
2. 1D Rate independent plasticity models
1.
2.
3.
4.
Return mapping algorithm
Consistent elastoplastic tangent modulus
Step by step algorithm
Nonlinear isotropic hardening
3. 1D Rate dependent plasticity models
1.
2.
3.
4.
Return mapping algorithm
Consistent elastoplastic tangent modulus
Step by step algorithm
Nonlinear isotropic hardening
4. 1D Computational plasticity assignment
April 1, 2015
Carlos Agelet de Saracibar
1D Plasticity Algorithms > Contents
Contents
Contents
1. Introduction
2. 1D Rate independent plasticity models
1.
2.
3.
4.
Return mapping algorithm
Consistent elastoplastic tangent modulus
Step by step algorithm
Nonlinear isotropic hardening
3. 1D Rate dependent plasticity models
1.
2.
3.
4.
Return mapping algorithm
Consistent elastoplastic tangent modulus
Step by step algorithm
Nonlinear isotropic hardening
4. 1D Computational plasticity assignment
April 1, 2015
Carlos Agelet de Saracibar
1D Plasticity Algorithms > Introduction
Time integration algorithm
En+1
Enp
April 1, 2015
Time integration
algorithm
Carlos Agelet de Saracibar
p
En+
1
1D Plasticity Algorithms > Contents
Contents
Contents
1. Introduction
2. 1D Rate independent plasticity models
1.
2.
3.
4.
Return mapping algorithm
Consistent elastoplastic tangent modulus
Step by step algorithm
Nonlinear isotropic hardening
3. 1D Rate dependent plasticity models
1.
2.
3.
4.
Return mapping algorithm
Consistent elastoplastic tangent modulus
Step by step algorithm
Nonlinear isotropic hardening
4. 1D Computational plasticity assignment
April 1, 2015
Carlos Agelet de Saracibar
1D Plasticity Algorithms > Rate Independent Plasticity Models
1D Rate independent plasticity model
1. Additive split of strains
E := Ee + E p, E :=
{ , 0, 0} ,
E p :=
p
e
,
E
{
} :=
{ , , }
2. Constitutive equations
S=
Ee (Ee ) =
C Ee , S :=
{ , q, q }, C := diag {E , K , H }
3. Associative plastic flow rule
p= f (S )
E
S
4. Yield function
f (S ) = q Y + q
5. Kuhn-Tucker loading/unloading conditions
0,
April 1, 2015
f (S) 0, f (S) =
0
Carlos Agelet de Saracibar
1D Plasticity Algorithms > Rate Independent Plasticity Models
Return mapping algorithm
Associative plastic flow rule: plastic strains at time n+1
p= f (S )
E
S
Using a Backward-Euler (BE) time integration scheme yields,
Enp+1 =Enp + n +1 Sn+1 f (Sn +1 )
np+1 =
np + n +1 sgn ( n +1 qn +1 )
n + n +1
n +=
1
n n +1 sgn ( n +1 qn +1 )
n +1 =
April 1, 2015
Carlos Agelet de Saracibar
1D Plasticity Algorithms > Rate Independent Plasticity Models
Return mapping algorithm
Constitutive equations: stress state at time n+1
Ee (Ee ) =
S=
C Ee =
C (E E p ) , C := diag { E , K , H }
The time-discrete constituve equation at time n+1 takes the
form,
=
Sn +1 C=
Een +1 C (En +1 Enp+1 )
Substituting the plastic strain variables at time n+1 yields,
(
C (E
p
S=
C
E
E
n +1
n +1
n n +1 Sn+1 f (Sn +1 )
April 1, 2015
E
n +1
n ) n +1C Sn+1 f (Sn +1 )
Carlos Agelet de Saracibar
1D Plasticity Algorithms > Rate Independent Plasticity Models
Return mapping algorithm
Trial state at time n+1
The trial state at time n+1 is defined by freezing the plastic
behaviour at the time step, yielding
p
:
Enp+,trial
=
E
1
n
trial
p ,trial
:
Een,+=
E
E
1
n +1
n +1
e ,trial
p ,trial
p
:
Strial
=
C
E
=
C
E
E
=
C
E
E
( n+1 n+1 ) ( n+1 n )
n +1
n +1
trial
f ntrial
f
:
=
S
( n+1 )
+1
April 1, 2015
Carlos Agelet de Saracibar
1D Plasticity Algorithms > Rate Independent Plasticity Models
Return mapping algorithm
Return mapping algorithm
The return mapping algorithm takes the form,
Sn +1 =Strial
n +1 n +1 C Sn+1 f (Sn +1 )
n +1 =
ntrial
+1 n +1 E sgn ( n +1 qn +1 )
trial
=
q
q
:
n +1
n +1 n +1 K
trial
=
q
q
:
n +1 + n +1 H sgn ( n +1 qn +1 )
n +1
April 1, 2015
Carlos Agelet de Saracibar
10
1D Plasticity Algorithms > Rate Independent Plasticity Models
Return mapping algorithm
1. Additive split of strains at time n+1
e
p
E=
:
E
+
E
n +1
n +1
n +1
2. Stresses at time n+1. Return mapping algorithm
Sn +1 = Strial
n +1 n +1C Sn+1 f (Sn +1 )
3. Plastic internal variables at time n+1
Enp+1 =Enp + n +1 Sn+1 f (Sn +1 )
4. Yield function at time n+1
f n +1 := f (Sn +1 ) = n +1 qn +1 Y + qn +1
5. Kuhn-Tucker loading/unloading conditions at time n+1
n +1 0,
April 1, 2015
f n +1 0, n +1 f n +1 =
0
Carlos Agelet de Saracibar
11
1D Plasticity Algorithms > Rate Independent Plasticity Models
Return mapping algorithm
Theorem 1. Elastic step/plastic step
If the yield function is convex and the constitutive matrix is
definite-positive, the following condition holds,
f (Strial
n +1 ) f (Sn +1 )
and Kuhn-Tucker loading/unloading conditions can be decided
just in terms of the trial yield function according to,
f (Strial
0 Elastic step
n +1 =
n +1 ) < 0
trial
f
S
n +1 > 0 Plastic step
(
n +1 ) > 0
April 1, 2015
Carlos Agelet de Saracibar
12
1D Plasticity Algorithms > Rate Independent Plasticity Models
Return mapping algorithm
Proof 1. Convexity of the yield function yields,
trial
f (Strial
f
S
S
( n+1 ) ( n+1 Sn+1 ) Sn+1 f (Sn+1 )
n +1 )
Using the return mapping equation, yields,
Sn +1 = Strial
n +1 n +1C Sn+1 f (Sn +1 )
f (Strial
n +1 ) f (Sn +1 ) n +1Sn+1 f (Sn +1 ) C Sn+1 f (Sn +1 )
Definite-positiveness of the constitutive matrix yields,
f (Strial
n +1 ) f (Sn +1 ) n +1Sn+1 f (Sn +1 ) C Sn+1 f (Sn +1 ) 0
f (Strial
n +1 ) f (Sn +1 )
April 1, 2015
Carlos Agelet de Saracibar
13
1D Plasticity Algorithms > Rate Independent Plasticity Models
Return mapping algorithm
Proof 2. The trial yield function at time n+1 determines the
discrete plastic loading/elastic unloading conditions according to,
if
then
f (Strial
n +1 ) < 0
f (Sn +1 ) f (Strial
n +1 ) < 0
0 n +1 =
0 Elastic step
n +1 f (Sn +1 ) =
else if
then
f (Strial
n +1 ) 0
Sn +1 = Strial
n +1 n +1C Sn+1 f (Sn +1 )
0 f (Sn +1 ) =
0 Plastic step
n +1 > 0, n +1 f (Sn +1 ) =
end if
April 1, 2015
Carlos Agelet de Saracibar
14
1D Plasticity Algorithms > Rate Independent Plasticity Models
Return mapping algorithm
Theorem 2. Closest-point-projection
The stress state at time n+1 is the closest-point-projection (cpp)
of the trial stress state at n+1 onto the space of admissible
stresses, measured in the complementary energy norm,
trial
S
=
arg
min
S
( n+1 S) S E
n +1
where the complementary energy is given by,
(S
trial
n +1
S
=
)
1
2
trial
n +1
2
C 1
1
trial
1
trial
=
S
S
C
S
) ( n+1 S)
n +1
2(
April 1, 2015
Carlos Agelet de Saracibar
15
1D Plasticity Algorithms > Rate Independent Plasticity Models
Return mapping algorithm
Proof 1. The stress state at time n+1, the closest-point-projection
of the trial stress state at n+1 onto the space of admissible
stresses, is the solution of the following constrained
minimization problem,
trial
S
=
arg
min
S
( n+1 S) S E
n +1
Using the Lagrange multipliers method, the constrained
minimization problem can be transformed into an unconstrained
minimization problem defined as,
Sn +1
L
April 1, 2015
arg min L
trial
S
( n+1 ,S; n+1 ) S
trial
trial
S
,
S;
:
=
S
( n+1
) ( n+1 S) + f (S)
Carlos Agelet de Saracibar
16
1D Plasticity Algorithms > Rate Independent Plasticity Models
Return mapping algorithm
The optimality conditions of the unconstrained minimization
problem read,
S L
n +1
(S
trial
n +1
S Strial
Sn +1 + n +1 S f (Sn +1 )
, Sn +1 ; n +1 :=
n +1
n +1
n +1
C 1 Strial
Sn +1 + n +1 S f (Sn +1 ) =
:=
0
n +1
n +1 0,
n +1
f (Sn +1 ) 0, n +1 f (Sn +1 ) =
0
The return mapping algorithm arises as the solution of the
unconstrained minimization problem yielding,
Sn +1 =Strial
n +1 n +1 C Sn+1 f (Sn +1 )
April 1, 2015
Carlos Agelet de Saracibar
17
1D Plasticity Algorithms > Rate Independent Plasticity Models
Return mapping algorithm
Return mapping algorithm
The return mapping algorithm takes the form,
Sn +1 =Strial
n +1 n +1 C Sn+1 f (Sn +1 )
n +1 =
ntrial
+1 n +1 E sgn ( n +1 qn +1 )
trial
=
q
q
:
n +1
n +1 n +1 K
trial
=
q
q
:
n +1 + n +1 H sgn ( n +1 qn +1 )
n +1
April 1, 2015
Carlos Agelet de Saracibar
18
1D Plasticity Algorithms > Rate Independent Plasticity Models
Return mapping algorithm
The solution for the return mapping algorithm yields,
trial
n +1 qn +1 = ntrial
q
+1
n +1 n +1 ( E + H ) sgn ( n +1 qn +1 )
trial
trial
trial
sgn
n +1 qn +1 sgn ( n +1 qn +1 ) = ntrial
q
( n+1 n+1 )
n +1
+1
n +1 ( E + H ) sgn ( n +1 qn +1 )
n +1
qn +1 + n +1 ( E + H ) ) sgn ( n +1 qn +1 ) =
trial
trial
trial
=
ntrial
q
sgn
q
(
+1
n +1
n +1
n +1 )
trial
n +1 qn +1 + n +1 ( E + H ) = ntrial
q
+1
n +1
trial
trial
sgn ( n +1 qn=
sgn
q
( n+1 n+1 )
+1 )
April 1, 2015
Carlos Agelet de Saracibar
19
1D Plasticity Algorithms > Rate Independent Plasticity Models
Return mapping algorithm
For the non-trivial case (plastic loading), using the discrete KuhnTucker loading/unloading conditions, the discrete plastic
multiplier (or discrete plastic consistency parameter) reads,
f ( n +1 , qn +1 , qn +1 ) 0, n +1 f ( n +1 , qn +1 , qn +1 ) =
0
n +1 0,
if
n +1 > 0 then
f ( n +1 , qn +1 , qn +1 ) =
0
f ( n +1 , qn +1 , qn +1 ) = n +1 qn +1 Y + qn +1
trial
trial
= ntrial
+
+
+
q
E
K
H
q
) Y n+1
+1
n +1
n +1 (
trial
trial
= f ( ntrial
,
q
,
q
0
+1
+1 ( E + K + H )
n +1
n +1 ) n=
trial
trial
n +1 = ( E + K + H ) f ( ntrial
,
q
,
q
+1
n +1
n +1 ) > 0
1
April 1, 2015
Carlos Agelet de Saracibar
20
1D Plasticity Algorithms > Rate Independent Plasticity Models
Return mapping algorithm
The return mapping algorithm takes the form,
n +1 =
ntrial
+1 n +1 E sgn ( n +1 qn +1 )
trial
=
q
q
:
n +1
n +1 n +1 K
trial
=
q
q
:
n +1 + n +1 H sgn ( n +1 qn +1 )
n +1
1
trial
trial
trial
trial
trial
n +1= ntrial
,
,
sgn
E
K
H
f
q
q
E
+
+
(
)
(
)
(
+1
n +1
n +1
n +1
n +1
n +1 )
trial
trial
trial
trial
:
,
,
q
q
E
K
H
f
q
q
=
+
+
(
)
( n+1 n+1 n+1 ) K
n +1
n +1
1
trial
trial
trial
trial
trial
trial
:
,
,
sgn
q
q
E
K
H
f
q
q
H
q
=
+
+
+
(
) ( n+1 n+1 n+1 )
n +1
( n+1 n+1 )
n +1
April 1, 2015
Carlos Agelet de Saracibar
21
1D Plasticity Algorithms > Rate Independent Plasticity Models
Return mapping algorithm
The update of the plastic internal variables takes the form,
np+1 =
np + n +1 sgn ( n +1 qn +1 )
n + n +1
n +=
1
n n +1 sgn ( n +1 qn +1 )
n +1 =
trial
trial
trial
trial
np+1 = np + ( E + K + H )1 f ( ntrial
q
q
q
,
,
sgn
)
(
+1
n +1
n +1
n +1
n +1 )
trial
trial
trial
=
+
+
+
E
K
H
f
q
q
,
,
(
)
( n+1 n+1 n+1 )
n +1
n
1
trial
trial
trial
trial
trial
=
+
+
E
K
H
f
q
q
q
,
,
sgn
(
) ( n+1 n+1 n+1 ) ( n+1 n+1 )
n +1
n
April 1, 2015
Carlos Agelet de Saracibar
22
1D Plasticity Algorithms > Rate Independent Plasticity Models
Consistent elastoplastic tangent modulus
Consistent elastoplastic tangent modulus
The consistent elastoplastic tangent modulus is computed taking
the variation of the stress at time n+1, yielding,
n +1=
trial
n +1
d n=
d
+1
trial
trial
( E + K + H ) f ntrial
q
sgn
( n+1 n+1 )
+1
trial
n +1
trial
trial
( E + K + H ) df ntrial
E
sgn
q
( n+1 n+1 )
+1
1
where
e ,trial
p
d ntrial
Ed
E
d
E d n +1
=
=
=
(
+1
n +1
n +1
n )
trial
trial
trial
trial
trial
trial
sgn
df ntrial
d
q
d
q
=
( n+1 n+1 ) ( n+1 n+1 )
+1
n +1
n +1
trial
trial
trial
trial
sgn ( ntrial
sgn
q
d
q
=
( n+1 n+1 ) E d n+1
+1
n +1 )
n +1
April 1, 2015
Carlos Agelet de Saracibar
23
1D Plasticity Algorithms > Rate Independent Plasticity Models
Consistent elastoplastic tangent modulus
The consistent elastoplastic tangent modulus takes the form,
d n=
d
+1
trial
n +1
trial
trial
( E + K + H ) df ntrial
q
sgn
( n+1 n+1 )
+1
1
d n +1 = E d n +1
trial
trial
trial
sgn
( E + K + H ) sgn ( ntrial
q
E
d
q
( n+1 n+1 )
+1
n +1 )
n +1
1
d n +1 = E 1 E ( E + K + H )
) d
n +1
d n +1 = E d n +1 , E := E 1 E ( E + K + H )
ep
April 1, 2015
ep
Carlos Agelet de Saracibar
)
24
1D Plasticity Algorithms > Rate Independent Plasticity Models
1D Plasticity algorithm
1D Plasticity algorithm
Step 1. Given the strain at time n+1 (strain driven problem), and
the stored plastic internal variables at time n (plastic internal
variables database)
Step 2. Compute the trial state at time n+1
p
:
=
Enp+,trial
E
1
n
trial
p ,trial
:
Een,+=
E
E
1
n +1
n +1
e ,trial
p ,trial
p
:
Strial
=
C
E
=
C
E
E
=
C
E
E
( n+1 n+1 ) ( n+1 n )
n +1
n +1
trial
trial
trial
f ntrial
q
:
=
+
+1
n +1
n +1
Y
n +1
April 1, 2015
Carlos Agelet de Saracibar
25
1D Plasticity Algorithms > Rate Independent Plasticity Models
1D Plasticity algorithm
Step 3. Check the trial yield function at time n+1
if f
trial
n +1
0 then set
( )n+1 =
( )n+1 ,
trial
E ep =
E and exit
Step 4. Compute discrete plastic multiplier at time n+1
n +1 = ( E + K + H ) f ntrial
+1
1
April 1, 2015
Carlos Agelet de Saracibar
26
1D Plasticity Algorithms > Rate Independent Plasticity Models
1D Plasticity algorithm
Step 5. Return mapping algorithm (closest-point-projection)
trial
Sn +1 =Strial
f
( n+1 )
n +1
n +1
Strial
n+1
1
trial
trial
trial
n +1= ntrial
+
+
E
K
H
f
E
sgn
q
(
)
(
+1
n +1
n +1
n +1 )
trial
trial
q
=
q
E
+
K
+
H
f
:
(
)
n +1
n +1
n +1 K
1
trial
trial
trial
trial
q
=
q
+
E
+
K
+
H
f
H
q
:
sgn
(
) n+1
( n+1 n+1 )
n +1
n +1
April 1, 2015
Carlos Agelet de Saracibar
27
1D Plasticity Algorithms > Rate Independent Plasticity Models
1D Plasticity algorithm
Step 6. Update plastic internal variables database at time n+1
Enp+1 =Enp + n +1 Strial f (Strial
n +1 )
n+1
trial
trial
np+1 = np + ( E + K + H )1 f ntrial
q
sgn
(
+1
n +1
n +1 )
trial
=
+
+
+
E
K
H
f
(
) n+1
n +1
n
1
trial
trial
trial
=
+
+
E
K
H
f
q
sgn
)
(
( n+1 n+1 )
n +1
n
n +1
Step 7. Compute the consistent elastoplastic tangent modulus
E := E 1 E ( E + K + H )
ep
April 1, 2015
Carlos Agelet de Saracibar
28
1D Plasticity Algorithms > Rate Independent Plasticity Models
Nonlinear isotropic hardening
Nonlinear isotropic hardening
q := = ( ) = ( )
Exponential saturation law + linear hardening
( ) = ( Y ) (1 exp ( ) ) + K
q := := ( Y ) (1 exp ( ) ) K
April 1, 2015
Carlos Agelet de Saracibar
29
1D Plasticity Algorithms > Rate Independent Plasticity Models
Nonlinear isotropic hardening
Time discrete nonlinear isotropic hardening
qn +1 := ( n +1 ) = ( n + n +1 )
trial
qntrial
=
:
( n+1 ) = (n ) = qn
+1
trial
q=
q
:
n +1
n +1 ( n + n +1 ) + ( n )
April 1, 2015
Carlos Agelet de Saracibar
30
1D Plasticity Algorithms > Rate Independent Plasticity Models
Nonlinear isotropic hardening
Plastic loading: Yield function at time n+1
n )) 0
=
f n +1 f ntrial
+1 n +1 ( E + H ) ( ( n + n +1 ) ( =
Nonlinear residual scalar equation on the plastic multiplier at
time n+1
g=
f=
0
n +1 : g ( n=
+1 )
n +1
trial
n )) 0
g=
:
f
n +1
n +1 n +1 ( E + H ) ( ( n + n +1 ) ( =
April 1, 2015
Carlos Agelet de Saracibar
31
1D Plasticity Algorithms > Rate Independent Plasticity Models
Nonlinear isotropic hardening
Newton-Raphson iterative solution algorithm
Step 1. Initialize iteration counter and plastic multiplier
=
k 0,=
nk+1 0
Step 2. Compute the residual g at time n+1, iteration k
Step 3. While the absolute value of the current residual at time
n+1, iteration k, is greater than a tolerance
Step 4. Solve the linarized equation
g nk+1 + Dg nk+1 nk+1 =
0
Step 5. Update the plastic multiplier at time n+1, iteration k+1
+1
k
k
nk=
+1
n +1
n +1
Step 6. Compute the residual g at time n+1, iteration k+1
Step 7. Increment iteration counter k=k+1 and go to Step 3
April 1, 2015
Carlos Agelet de Saracibar
32
1D Plasticity Algorithms > Rate Independent Plasticity Models
Nonlinear isotropic hardening
Newton-Raphson iterative solution algorithm
g nk+1 + Dg nk+1 nk+1 =
0
k
trial
k
k
g=
:
f
E
H
)
( n n+1 ) (n )
n +1
n +1
n +1 (
Dg nk+1 nk+1 := ( E + H ) nk+1 ( n + nk+1 ) nk+1
:= E + ( n + nk+1 ) + H nk+1
nk+1 := nk++11 nk+1
April 1, 2015
Carlos Agelet de Saracibar
33
1D Plasticity Algorithms > Rate Independent Plasticity Models
Nonlinear isotropic hardening
Consistent elastoplastic tangent modulus
The consistent elastoplastic tangent modulus is computed taking
the variation of the stress at time n+1, yielding,
trial
trial
n +1 =
ntrial
E
sgn
q
( n+1 n+1 )
+1
n +1
trial
trial
d n +1 =
d ntrial
q
sgn
( n+1 n+1 )
+1
n +1
where the variations of the trial stress tensor, and plastic
multiplier, at time n+1, have to be computed.
April 1, 2015
Carlos Agelet de Saracibar
34
1D Plasticity Algorithms > Rate Independent Plasticity Models
Nonlinear isotropic hardening
The variation of the plastic multiplier at time n+1 is computed
setting the variation of the residual at time n+1 equal to zero,
trial
n) 0
g=
f
:
n +1
n +1 n +1 ( E + H ) ( n + n +1 ) ( =
trial
=
dg
df
:
n +1
n +1 d n +1 ( E + H ) ( n + n +1 ) d n +1
=: df ntrial
H 0
+1 d n +1 E + ( n + n +1 ) +=
d n +=
1
April 1, 2015
E + ( n + n +1 ) + H
df ntrial
+1
Carlos Agelet de Saracibar
35
1D Plasticity Algorithms > Rate Independent Plasticity Models
Nonlinear isotropic hardening
Substituting the variations shown before, the following discrete
tangent constitutive equation is obtained,
d n +1 = Enep+1d n +1
where the consistent elastoplastic tangent modulus at time n+1
is given by
E= E 1 E E + ( n + n +1 ) + H
ep
n +1
April 1, 2015
Carlos Agelet de Saracibar
36
1D Plasticity Algorithms > Contents
Contents
Contents
1. Introduction
2. 1D Rate independent plasticity models
1.
2.
3.
4.
Return mapping algorithm
Consistent elastoplastic tangent modulus
Step by step algorithm
Nonlinear isotropic hardening
3. 1D Rate dependent plasticity models
1.
2.
3.
4.
Return mapping algorithm
Consistent elastoplastic tangent modulus
Step by step algorithm
Nonlinear isotropic hardening
4. 1D Computational plasticity assignment
April 1, 2015
Carlos Agelet de Saracibar
37
1D Plasticity Algorithms > Rate Dependent Plasticity Models
1D Rate dependent plasticity model
1. Additive split of strains
E := Ee + E p, E :=
{ , 0, 0} ,
E p :=
p
e
,
E
{
} :=
{ , , }
2. Constitutive equations
S=
Ee (Ee ) =
C Ee , S :=
{ , q, q }, C := diag {E , K , H }
3. Associative plastic flow rule
p= f (S )
E
S
4. Yield function
f (S ) = q Y + q
5. Plastic multiplier
April 1, 2015
f (S ) 0
Carlos Agelet de Saracibar
38
1D Plasticity Algorithms > Rate Dependent Plasticity Models
Return mapping algorithm
Associative plastic flow rule: plastic strains at time n+1
p= f (S )
E
S
Using a Backward-Euler (BE) time integration scheme yields,
p
Enp+=
E
1
n + n +1 t Sn+1 f (Sn +1 )
April 1, 2015
Carlos Agelet de Saracibar
39
1D Plasticity Algorithms > Rate Dependent Plasticity Models
Return mapping algorithm
Constitutive equations: stress state at time n+1
Ee (Ee ) =
S=
C Ee =
C (E E p ) , C := diag { E , K , H }
The time-discrete constituve equation at time n+1 takes the
form,
=
Sn +1 C=
Een +1 C (En +1 Enp+1 )
Substituting the plastic strains at time n+1 yields,
(
C (E
=
Sn +1 C En +1 Enp n +1 t Sn+1 f (Sn +1 )
=
April 1, 2015
E
n +1
n ) n +1t C Sn+1 f (Sn +1 )
Carlos Agelet de Saracibar
40
1D Plasticity Algorithms > Rate Dependent Plasticity Models
Return mapping algorithm
Trial state at time n+1
The trial state at time n+1 is defined by freezing the plastic
behaviour at the time step, yielding
p
Enp+,trial
:
=
E
1
n
trial
p ,trial
Een,+=
:
E
E
1
n +1
n +1
e ,trial
p ,trial
p
Strial
:
=
C
E
=
C
E
E
=
C
E
E
( n+1 n+1 ) ( n+1 n )
n +1
n +1
trial
f nrial
f
:
=
S
( n+1 )
+1
April 1, 2015
Carlos Agelet de Saracibar
41
1D Plasticity Algorithms > Rate Dependent Plasticity Models
Return mapping algorithm
Return mapping algorithm
The return mapping algorithm takes the form,
trial
Sn +=
S
1
n +1 n +1 t C Sn+1 f (Sn +1 )
n +1 = ntrial
+1 n +1 t E sgn ( n +1 qn +1 )
trial
q
:
q
=
n +1
n +1 n +1 t K
trial
q
:
q
=
n +1 + n +1 t H sgn ( n +1 qn +1 )
n +1
April 1, 2015
Carlos Agelet de Saracibar
42
1D Plasticity Algorithms > Rate Dependent Plasticity Models
Return mapping algorithm
1. Additive split of strains at time n+1
e
p
E=
:
E
+
E
n +1
n +1
n +1
2. Stresses at time n+1. Return mapping algorithm
trial
Sn +=
S
1
n +1 n +1t C Sn+1 f (Sn +1 )
3. Plastic internal variables at time n+1
p
Enp+=
E
1
n + n +1 t Sn+1 f (Sn +1 )
4. Yield function at time n+1
f n +1 := f (Sn +1 ) = n +1 qn +1 Y + qn +1
5. Plastic multiplier at time n+1
=
n +1
April 1, 2015
f (Sn +1 ) 0
Carlos Agelet de Saracibar
43
1D Plasticity Algorithms > Rate Dependent Plasticity Models
Return mapping algorithm
Return mapping algorithm
The return mapping algorithm takes the form,
trial
Sn +=
S
1
n +1 n +1 t C Sn+1 f (Sn +1 )
n +1 = ntrial
+1 n +1 t E sgn ( n +1 qn +1 )
trial
q
:
q
=
n +1
n +1 n +1 t K
trial
q
:
q
=
n +1 + n +1 t H sgn ( n +1 qn +1 )
n +1
April 1, 2015
Carlos Agelet de Saracibar
44
1D Plasticity Algorithms > Rate Dependent Plasticity Models
Return mapping algorithm
The solution for the return mapping algorithm yields,
trial
n +1 qn +1 = ntrial
q
+1
n +1 n +1t ( E + H ) sgn ( n +1 qn +1 )
trial
trial
trial
n +1 qn +1 sgn ( n +1 qn +1 ) = ntrial
q
sgn
q
( n+1 n+1 )
+1
n +1
n +1t ( E + H ) sgn ( n +1 qn +1 )
n +1
qn +1 + n +1t ( E + H ) ) sgn ( n +1 qn +1 ) =
trial
trial
trial
=
ntrial
q
sgn
q
(
+1
n +1
n +1
n +1 )
trial
n +1 qn +1 + n +1t ( E + H )= ntrial
q
+1
n +1
trial
trial
sgn ( n +1 qn=
sgn
q
( n+1 n+1 )
+1 )
April 1, 2015
Carlos Agelet de Saracibar
45
1D Plasticity Algorithms > Rate Dependent Plasticity Models
Return mapping algorithm
For the non-trivial case (plastic loading), the yield function and
the discrete plastic multiplier are greater than zero, and the
following expressions hold,
if
n +1 > 0 then
f ( n +1 , qn +1 , qn +1 ) =
n +1 > 0
f ( n +1 , qn +1 , qn +1 ) = n +1 qn +1 Y + qn +1
trial
trial
q
t
E
K
H
q
= ntrial
+
+
+
(
)
n +1
n +1
Y
n +1
+1
trial
trial
,
,
q
q
n +1
= f ( ntrial
n +1
n +1 ) =
n +1t ( E+K +H )
+1
trial
trial
,
,
n +1t= E + K + H + f ( ntrial
q
q
n +1
n +1 ) > 0
+1
t
April 1, 2015
Carlos Agelet de Saracibar
46
1D Plasticity Algorithms > Rate Dependent Plasticity Models
Return mapping algorithm
The return mapping algorithm takes the form,
n +1 = ntrial
+1 n +1 t E sgn ( n +1 qn +1 )
trial
:
q
=
q
n +1
n +1 n +1 t K
trial
:
q
q
=
n +1 + n +1 t H sgn ( n +1 qn +1 )
n +1
1
trial
trial
trial
n +1= ntrial
sgn
E
K
H
t
f
E
+
+
+
(
)
(
n +1
n +1
n +1 )
+1
trial
trial
q
:
q
E
K
H
t
f
=
+
+
+
(
) n+1 K
n +1
n +1
1
trial
trial
trial
trial
q
:
q
E
K
H
t
f
H
sgn
q
=
+
+
+
+
(
)
( n+1 n+1 )
n +1
n +1
n +1
April 1, 2015
Carlos Agelet de Saracibar
47
1D Plasticity Algorithms > Rate Dependent Plasticity Models
Return mapping algorithm
The update of the plastic internal variables takes the form,
np+1 = np + n +1t sgn ( n +1 qn +1 )
n +1 := n + n +1t
n +1 := n n +1t sgn ( n +1 qn +1 )
trial
trial
np+1 = np + ( E + K + H + t )1 f ntrial
sgn
q
(
+1
n +1
n +1 )
trial
=
+
+
+
+
E
K
H
t
f
:
(
)
n +1
n
n +1
1
trial
trial
trial
=
+
+
+
:
E
K
H
t
f
sgn
q
(
)
(
n +1
n
n +1
n +1
n +1 )
April 1, 2015
Carlos Agelet de Saracibar
48
1D Plasticity Algorithms > Rate Dependent Plasticity Models
Consistent elastoplastic tangent modulus
Consistent elastoplastic tangent modulus
The consistent elastoplastic tangent modulus is computed taking
the variation of the stress at time n+1, yielding,
trial
trial
trial
+
+
+
n +1= ntrial
E
K
H
t
f
E
sgn
q
(
)
(
+1
n +1
n +1
n +1 )
1
d n=
d
+1
trial
n +1
trial
trial
( E + K + H + t ) df ntrial
E
sgn
q
(
+1
n +1
n +1 )
1
where
e ,trial
p
d ntrial
=
Ed
=
E
d
=
( n+1 n ) E d n+1
+1
n +1
trial
trial
trial
trial
trial
trial
sgn
df ntrial
=
d
q
=
q
d
q
( n+1 n+1 ) ( n+1 n+1 )
+1
n +1
n +1
trial
trial
trial
trial
sgn ( ntrial
sgn
=
q
d
q
)
(
+1
n +1
n +1
n +1
n +1 ) E d n +1
April 1, 2015
Carlos Agelet de Saracibar
49
1D Plasticity Algorithms > Rate Dependent Plasticity Models
Consistent elastoplastic tangent modulus
The consistent elastoplastic tangent modulus takes the form,
d n=
d
+1
trial
n +1
trial
trial
sgn
E
q
( E + K + H + t ) df ntrial
( n+1 n+1 )
+1
1
d n +1 = E d n +1
trial
q
( E + K + H + t ) sgn 2 ( ntrial
n +1 ) E E d n +1
+1
1
d n +1 = E 1 E ( E + K + H + t )
) d
n +1
d n +1 = E ep d n +1 , E ep := E 1 E ( E + K + H + t )
April 1, 2015
Carlos Agelet de Saracibar
)
50
1D Plasticity Algorithms > Rate Dependent Plasticity Models
1D Plasticity algorithm
1D Plasticity algorithm
Step 1. Given the strain at time n+1 (strain driven problem), and
the stored plastic internal variables at time n (plastic internal
variables database)
Step 2. Compute the trial state at time n+1
p
:
=
Enp+,trial
E
1
n
trial
p ,trial
:
Een,+=
E
E
1
n +1
n +1
e ,trial
p ,trial
p
:
Strial
=
C
E
=
C
E
E
=
C
E
E
(
)
(
n +1
n +1
n +1
n +1
n +1
n )
trial
trial
trial
f ntrial
q
:
=
+
+1
n +1
n +1
Y
n +1
April 1, 2015
Carlos Agelet de Saracibar
51
1D Plasticity Algorithms > Rate Dependent Plasticity Models
1D Plasticity algorithm
Step 3. Check the trial yield function at time n+1
if f
trial
n +1
0 then set
( )n+1 =
( )n+1 ,
E ep =
E and exit
trial
Step 4. Compute discrete plastic multipliet at time n+1
n +1 = ( E + K + H + t ) f ntrial
+1
1
April 1, 2015
Carlos Agelet de Saracibar
52
1D Plasticity Algorithms > Rate Dependent Plasticity Models
1D Plasticity algorithm
Step 5. Return mapping algorithm (closest-point-projection)
trial
Sn +1 =Strial
C
S
f
(
n +1
n +1
n +1 )
Strial
n+1
1
trial
trial
trial
n +1= ntrial
+
+
+
E
K
H
t
f
E
sgn
q
(
) n+1
(
+1
n +1
n +1 )
trial
trial
=
+
+
+
q
:
q
E
K
H
t
f
(
)
n +1
n +1
n +1 K
1
trial
trial
trial
trial
=
+
+
+
+
q
:
q
E
K
H
t
f
H
sgn
q
(
) n+1
( n+1 n+1 )
n +1
n +1
April 1, 2015
Carlos Agelet de Saracibar
53
1D Plasticity Algorithms > Rate Dependent Plasticity Models
1D Plasticity algorithm
Step 6. Update plastic internal variables database at time n+1
Enp+1 =Enp + n +1 Strial f (Strial
n +1 )
n+1
trial
trial
np+1 = np + ( E + K + H + t )1 f ntrial
sgn
q
(
n +1
n +1 )
+1
trial
E
K
H
t
f
=
+
+
+
+
(
) n+1
n +1
n
1
trial
trial
trial
E
K
H
t
f
sgn
q
=
+
+
+
(
)
( n+1 n+1 )
n +1
n
n +1
Step 7. Compute the consistent elastoplastic tangent modulus
E := E 1 E ( E + K + H + t )
ep
April 1, 2015
Carlos Agelet de Saracibar
54
1D Plasticity Algorithms > Rate Dependent Plasticity Models
Nonlinear isotropic hardening
Nonlinear isotropic hardening
q := = ( ) = ( )
Exponential saturation law + linear hardening
( ) = ( Y ) (1 exp ( ) ) + K
q := := ( Y ) (1 exp ( ) ) K
April 1, 2015
Carlos Agelet de Saracibar
55
1D Plasticity Algorithms > Rate Dependent Plasticity Models
Nonlinear isotropic hardening
Time discrete nonlinear isotropic hardening
qn +1 := ( n +1 ) = ( n + n +1t )
trial
:
qntrial
( n+1 ) = (n ) = qn
+1
trial
:
q=
q
n +1
n +1 ( n + n +1t ) + ( n )
April 1, 2015
Carlos Agelet de Saracibar
56
1D Plasticity Algorithms > Rate Dependent Plasticity Models
Nonlinear isotropic hardening
Plastic loading: Yield function at time n+1
n ) ) n +1 > 0
=
f n +1 f ntrial
+1 n +1t ( E + H ) ( ( n + n +1t ) ( =
Nonlinear residual scalar equation on the plastic multiplier at
time n+1
g n +1 :=g ( n +1 ) =f n +1 n +1 =0
g=
f
n +1 :
trial
n +1
April 1, 2015
n +1t E + H + ( ( n + n +1t ) (=
0
n ))
t
Carlos Agelet de Saracibar
57
1D Plasticity Algorithms > Rate Dependent Plasticity Models
Nonlinear isotropic hardening
Newton-Raphson iterative solution algorithm
Step 1. Initialize iteration counter and plastic multiplier
=
k 0,=
nk+1 0
Step 2. Compute the residual g at time n+1, iteration k
Step 3. While the absolute value of the current residual at time
n+1, iteration k, is greater than a tolerance
Step 4. Solve the linarized equation
g nk+1 + Dg nk+1 nk+1 =
0
Step 5. Update the plastic multiplier at time n+1, iteration k+1
+1
k
k
nk=
+1
n +1
n +1
Step 6. Compute the residual g at time n+1, iteration k+1
Step 7. Increment iteration counter k=k+1 and go to Step 3
April 1, 2015
Carlos Agelet de Saracibar
58
1D Plasticity Algorithms > Rate Dependent Plasticity Models
Nonlinear isotropic hardening
Newton-Raphson iterative solution algorithm
g nk+1 + Dg nk+1 nk+1 =
0
t E + H + ( n + nk+1t ) ( n )
g=: f
t
k
k
E + H + n +1t ( n + nk+1t ) nk+1t
Dg n +1 n +1 :=
t
k
:= E + ( n + n +1t ) + H + n +1t
t
nk+1 := nk++11 nk+1
k
n +1
April 1, 2015
trial
n +1
k
n +1
Carlos Agelet de Saracibar
59
1D Plasticity Algorithms > Rate Dependent Plasticity Models
Nonlinear isotropic hardening
Consistent elastoplastic tangent modulus
The consistent elastoplastic tangent modulus is computed taking
the variation of the stress at time n+1, yielding,
trial
trial
sgn
n +1 = ntrial
t
E
q
( n+1 n+1 )
n +1
+1
trial
trial
sgn
d n +1 = d ntrial
d
t
E
( n+1 n+1 )
+1
n +1
where the variations of the trial stress tensor, and plastic
multiplier, at time n+1, have to be computed.
April 1, 2015
Carlos Agelet de Saracibar
60
1D Plasticity Algorithms > Rate Dependent Plasticity Models
Nonlinear isotropic hardening
The variation of the plastic multiplier at time n+1 is computed
setting the variation of the residual at time n+1 equal to zero,
g=
f
n +1 :
trial
n +1
=
dg
n +1 : df
n +1t E + H + ( n + n +1t ) (=
0
n)
t
trial
n +1
=: df ntrial
+1
d n +1t E + H + ( n + n +1t ) d n +1t
t
d n +1t E + ( n + n +1t ) + H + =
0
t
d n +1=
t E + ( n + n +1t ) + H + df ntrial
+1
April 1, 2015
Carlos Agelet de Saracibar
61
1D Plasticity Algorithms > Rate Dependent Plasticity Models
Nonlinear isotropic hardening
Substituting the variations shown before, the following discrete
tangent constitutive equation is obtained,
d n +1 = Enep+1d n +1
where the consistent elastoplastic tangent modulus at time n+1
is given by
1
ep
E=
E 1 E E + ( n + n +1t ) + H +
n +1
April 1, 2015
Carlos Agelet de Saracibar
62
1D Plasticity Algorithms > Contents
Contents
Contents
1. Introduction
2. 1D Rate independent plasticity models
1.
2.
3.
4.
Return mapping algorithm
Consistent elastoplastic tangent modulus
Step by step algorithm
Nonlinear isotropic hardening
3. 1D Rate dependent plasticity models
1.
2.
3.
4.
Return mapping algorithm
Consistent elastoplastic tangent modulus
Step by step algorithm
Nonlinear isotropic hardening
4. 1D Computational plasticity assignment
April 1, 2015
Carlos Agelet de Saracibar
63
1D Plasticity Algorithms > 1D Computational Plasticity Assignment
1D Computational plasticity assignment
Implement in MATLAB the BE time-stepping algorithm for 1D
rate-independent/rate-dependent hardening plasticity
models, including linear and nonlinear isotropic hardening,
and linear kinematic hardening
Perform the numerical simulation of uniaxial cyclic plastic
loading/elastic unloading examples for the following cases:
o Rate-independent/rate-dependent perfect plasticity
o Rate-independent/rate-dependent linear isotropic hardening plasticity
o Rate-independent/rate-dependent nonlinear isotropic hardening
plasticity, considering an exponential saturation law
o Rate-independent/rate-dependent linear kinematic hardening
plasticity
o Rate-independent/rate-dependent nonlinear isotropic and linear
kinematic hardening plasticity
April 1, 2015
Carlos Agelet de Saracibar
64
1D Plasticity Algorithms > 1D Computational Plasticity Assignment
1D Computational plasticity assignment
For the perfect plasticity models, plot the stress-strain curves
For the linear isotropic/linear kinematic hardening models,
plot the stress-strain curves showing the influence of the
isotropic/kinematic hardening parameters
For the nonlinear isotropic hardening model, plot the stressstrain curves showing the influence of the exponential
coefficient of the exponential saturation law on the stressstrain curves
For the rate-dependent plasticity models, plot the stressstrain, and the stress-time curves showing the influence of the
viscosity parameter and the loading rate.
Show that the rate-independent response can be recovered
from the rate-dependent results using very small values for
the viscosity or the loading
rate
April 1, 2015
Carlos Agelet de Saracibar
65
1D Plasticity Algorithms > 1D Computational Plasticity Assignment
1D Computational plasticity assignment
Write a comprehensive deliverable report (10 pages)
providing the data of the cyclic loading and material
properties considered, the stress-strain curves, and the
stress-time curves for the rate-dependent plasticity examples.
Add suitable comments on the results, comparing the
influence of the different material parameters and loading
conditions.
Add a printed copy of the subroutines as an Appendix
April 1, 2015
Carlos Agelet de Saracibar
66
