6.6 Stress Update Algorithm

Back

6.6.1 Derivation

The equations describing viscoplasticity were summarized in the previous section. Ultimately,

we would like to transform these equations into a numerical algorithm yielding the

solution at time-step n + 1 if the solution at t = tn is known. To this end, the backward

288 FINITE STRAIN ELASTOPLASTICITY

Euler rule will be applied to the time derivatives. It is an implicit unconditionally stable

scheme expressing the time derivative at t = tn+1 in terms of the function values at t = tn

and t = tn+1:

( ˙ f )n+1 ≈ fn+1 − fn

_t

. (6.143)

Applying this to the flow rule, Equation (6.136), yields

μJ

−2/3

n+1 DEVn+1(Cp−1

n+1

− Cp−1

n ) = −2(γn+1 − γn)μn+1Nn+1 (6.144)

where

DEVn+1A := A − 13

(A : Cn+1)C

−1

n+1. (6.145)

Henceforth, the following definition will be used:

_γn+1 := γn+1 − γn. (6.146)

The hyperelasticity equation, Equation (6.132), yields

DEVn+1Sn+1 = μJ

−2/3

n+1 DEVn+1(Cp−1

n+1). (6.147)

Consequently, the flow rule, Equation (6.144), can be written as

DEVn+1Sn+1 = μJ

−2/3

n+1 DEVn+1Cp−1

n

− 2_γn+1μn+1Nn+1. (6.148)

Equations (6.138) and (6.137) yield

μn+1 = 13

μJ

−2/3

n+1 TRn+1Cp−1

n+1

+ 13

J

−2/3

n+1 TRn+1Q2,n+1. (6.149)

The auxiliary equations, Equations (6.140) and (6.142), lead to

Cp−1

n+1 : Cn+1 = Cp−1

n : Cn+1 (6.150)

Q2,n+1 : Cn+1 = Q2,n : Cn+1. (6.151)

Hence, Equation (6.149) can be transformed into

μn+1 = 13

μJ

−2/3

n+1 TRn+1Cp−1

n

+ 13

J

−2/3

n+1 TRn+1Q2,n. (6.152)

In a similar way, Equation (6.141) is transformed into

J

−2/3

n+1 DEVn+1(Q2,n+1) = J

−2/3

n+1 DEVn+1(Q2,n) −

2h

eq’

2

3μ

_γn+1μn+1Nn+1. (6.153)

Defining

T trial

n+1 := μJ

−2/3

n+1 DEVn+1Cp−1

n (6.154)

Atrial

n+1 := −J

−2/3

n+1 DEVn+1(Q2,n) (6.155)

_trial

n+1 := T trial

n+1

− Atrial

n+1 (6.156)

T n+1 := DEVn+1Sn+1 (6.157)

An+1 := −J

−2/3

n+1 DEVn+1(Q2,n+1) (6.158)

_n+1 := T n+1 − An+1 (6.159)

FINITE STRAIN ELASTOPLASTICITY 289

the evolution equations, Equations (6.148) and (6.153), can be rewritten as

T n+1 = T trial

n+1

− 2_γn+1μn+1Nn+1 (6.160)

An+1 = Atrial

n+1

+

2h

eq’

2

3μ

_γn+1μn+1Nn+1. (6.161)

Notice that as soon as the displacement field at t = tn+1 is known, the trial functions can

be calculated. They represent the stress state at t = tn+1 in the assumption that step n + 1

is purely elastic. Next, the yield-surface condition is checked (Equation (6.133)). Using

Equations (6.157), (6.158) and (6.159), the yield condition at t = tn+1 can be expressed as

_n+1 ≤ _23

h1(_

peq

n+1). (6.162)

Assuming at first that there is no plastic flow, or

_n+1 = _trial

n+1 (6.163)

_

peq

n+1

= _

peq

n (6.164)

Equation (6.162) reduces to

_trial

n+1

 ≤ _23

h1(_

peq

n ). (6.165)

If this equation is satisfied, the state is purely elastic and the trial functions are the solution.

Furthermore, the plastic internal variables do not change. If, on the other hand,

Equation (6.165) is not satisfied, plastic deformation occurs and Equation (6.134) applies,

where f = 0 for nonviscous deformation. At t = tn+1, this equation reads

_n+1 = _23

h1(_

peq

n+1) + _23

f (_

peq

n+1). (6.166)

Now _n+1 satisfies

_n+1 = T n+1 − An+1 (6.167)

= T trial

n+1

− Atrial

n+1

− 2μn+1

_

1 + h

eq’

2

3μ

_

_γn+1Nn+1 (6.168)

= _trial

n+1

− 2μn+1

_

1 + h

eq’

2

3μ

_

_γn+1

_n+1

_n+1. (6.169)

Equation (6.169) reveals that _n+1 and _trial

n+1 are parallel, and accordingly,

Nn+1 =

_trial

n+1

_trial

n+1

 

(6.170)

which means that Nn+1 can be calculated using the trial state. Substituting Equation (6.170)

into Equation (6.169) and taking the norm, one finds that

_n+1 = _trial

n+1

 − 2μn+1

_

1 + h

eq’

2

3μ

_

_γn+1 (6.171)

290 FINITE STRAIN ELASTOPLASTICITY

and the yield condition, Equation (6.166) leads to

_trial

n+1

 − 2μn+1

_

1 +

h

eq’

2 (_

peq

n+1)

3μ

_

_γn+1 = _23

h1(_

peq

n+1) + _23

f (˙_

peq

n+1). (6.172)

In this equation, _trial

n+1

 and μn+1 are known, _

peq

n+1 and _γn+1 are unknowns, related by

Equation (6.121) , or, equivalently,

_

peq

n+1

= _

peq

n + _23

_γn+1. (6.173)

Accordingly, Equation (6.172) yields

_trial

n+1

 − 2μn+1



1 +

h

eq’

2

_

_

peq

n + _23

_γn+1

_

3μ



_γn+1

= _23

h1

_

_

peq

n + _23

_γn+1

_ + _23

f

__23

_γn+1

_

. (6.174)

This is a nonlinear equation in _γn+1, which can be solved by the Newton–Raphson

technique. Once _γn+1 is known, all other quantities can be calculated using the equations

in this section. Indeed, the definition of DEV leads to

DEVn+1Cp−1

n+1

= Cp−1

n+1

− 13

(Cp−1

n+1 : Cn+1)C

−1

n+1 (6.175)

DEVn+1Cp−1

n

= Cp−1

n

− 13

(Cp−1

n : Cn+1)C

−1

n+1 (6.176)

and consequently, since

TRn+1Cp−1

n+1

= TRn+1Cp−1

n , (6.177)

we find

DEVn+1Cp−1

n+1

− DEVn+1Cp−1

n

= Cp−1

n+1

− Cp−1

n . (6.178)

Equation (6.144) can be transformed into

DEVn+1Cp−1

n+1

− DEVn+1Cp−1

n

= −2_γn+1

μn+1

μ

J

−2/3

n+1 Nn+1. (6.179)

Hence,

Cp−1

n+1

= Cp−1

n

− 2_γn+1

μn+1

μ

J

−2/3

n+1 Nn+1. (6.180)

Similarly, Equation (6.153) yields

Q2,n+1 = Q2,n

− 2heq’

3μ

_γn+1μn+1J

−2/3

n+1 Nn+1. (6.181)

Finally, the second Piola–Kirchhoff stress follows from Equations (6.148) and (6.132):

Sn+1 = T trial

n+1

− 2_γn+1μn+1Nn+1 + K

2

(J 2

n+1

− 1)C

−1

n+1. (6.182)

FINITE STRAIN ELASTOPLASTICITY 291

6.6.2 Summary

Given: Cp−1

n , Q2,n, γn, _

peq

n , Un+1.

1. Step 1: Geometric update

Cn+1, Jn+1, C

−1

n+1.

2. Step 2: Elastic prediction

TRn+1Cp−1

n

= Cp−1

n : Cn+1 (6.183)

DEVn+1Cp−1

n

= Cp−1

n

− 13

(TRn+1Cp−1

n )C

−1

n+1 (6.184)

TRn+1(Q2,n) = Q2,n : Cn+1 (6.185)

DEVn+1(Q2,n) = Q2,n

− 13

TRn+1(Q2,n)C

−1

n+1 (6.186)

T trial

n+1

= μJ

−2/3

n+1 DEVn+1Cp−1

n (6.187)

Atrial

n+1

= −J

−2/3

n+1 DEVn+1(Q2,n) (6.188)

_trial

n+1

= T trial

n+1

− Atrial

n+1. (6.189)

3. Step 3: Check for yielding

If

_trial

n+1

 − _23

h1(_

peq

n ) ≤ 0 (6.190)

(·)n+1 = (·)trial

n+1 and EXIT.

4. Step 4: Radial return scheme

μn+1 = 13

J

−2/3

n+1

_μTRn+1(Cp−1

n ) + TRn+1(Q2,n)_ (6.191)

Nn+1 =

_trial

n+1

_trial

n+1

 

(6.192)

_trial

n+1

 − 2μn+1



1 +

h

eq’

2

_

_

peq

n + _23

_γn+1

_

3μ



_γn+1

= _23

h1

_

_

peq

n + _23

_γn+1

_ + _23

f

__23

_γn+1

_ (6.193)

from which _γn+1 can be determined.

292 FINITE STRAIN ELASTOPLASTICITY

5. Update of the plastic state variables

_

peq

n+1

= _

peq

n + _23

_γn+1 (6.194)

Cp−1

n+1

= Cp−1

n

− 2_γn+1

μn+1

μ

J

2/3

n+1Nn+1 (6.195)

Q2,n+1 = Qn

−

2h

eq’

2

3μ

_γn+1μn+1J

2/3

n+1Nn+1. (6.196)

6. Update of the stress

Sn+1 = T trial

n+1

− 2_γn+1μn+1Nn+1 + K

2

(J 2

n+1

− 1)C

−1

n+1. (6.197)

Sometimes, Equation (6.193) is written in a different way. Defining

f trial

n+1 := _trial

n+1

 − _23

h1(_

peq

n ) (6.198)

one gets

f trial

n+1

= _23

           

h1

_

_

peq

n + _23

_γn+1

_ − h1(_

peq

n )

 

+ _23

f

__23

_γn+1

_ + 2μn+1



1 +

h

eq’

2

_

_

peq

n + _23

_γn+1

_

3μ

_

γn+1. (

6.199)

For linear hardening and creep laws of the form

h1

_

_

peq

n + _23

_γn+1

_ = h1(_

peq

n ) + h

 

1_γn+1_23

(6.200)

h

eq

2

_

_

peq

n + _23

_γn+1

_ = h

eq

2 (_

peq

n ) + h

eq

2 _γn+1_23

(6.201)

f

__23

_γn+1

_ = η_23

_γn+1 (6.202)

where h

 

1, h

eq

2 and η are constants. Equation (6.199) further reduces to

2μ_γn+1 =

f trial

n+1

1 + h

eq’

2

3μ

+ h

 

1

3μ

+ η

3μ

. (6.203)

For nonlinear laws, Equation (6.199) is first written as

g(_γn+1) := f trial

n+1

− _23

           

h1

_

_

peq

n + _23

_γn+1

_ − h1(_

peq

n )

 − _23

f

__23

_γn+1

_

− μn+1

_2_γn+1 + _23

           

h

eq

2

_

_

peq

n + _23

_γn+1

_ − h

eq

2 (_

peq

n )

 1

μ

_ (6.204)

FINITE STRAIN ELASTOPLASTICITY 293

since (backward Euler)

h

eq

2

_

_

peq

n + _23

_γn+1

_ ≈

h

eq

2

_

_

peq

n + _23

_γn+1

_ − h

eq

2 __

peq

n _

_23

_γn+1

. (6.205)

For a Newton–Raphson type solution of Equation (6.204), the derivative of g is also

needed:

dg

d(_γn+1)

= −23

h

 

1

_

_

peq

n + _23

_γn+1

_ − 23

f

 __23

_γn+1

_

− 2μn+1

            1 + _ 1

3μ

_

h

eq

2

_

_

peq

n + _23

_γn+1

_ (6.206)

where h

 

1, h

eq

2 and f

 denote derivatives with respect to their arguments. Since h1, h

eq

2

and f are user-defined functions, the derivatives can be determined too (analytically or

numerically). The Newton–Raphson scheme can be started with _γ(0)

n+1

= 0. Subsequent

iterations yield

_γ(k+1)

n+1

= _γ(k)

n+1

−

g(_γ (k)

n+1)

g(_γ (k)

n+1)

(6.207)

until _γ(final)

n+1 is small enough. Occasionally, depending on the form of the creep and

hardening functions, the Newton–Raphson procedure does not converge (cf Chapter 3).

Then, other techniques such as bisection (Press et al. 1990) (LЁuhrs et al. 1997) can be used.

6.6.3 Expansion of a thick-walled cylinder

Consider the expansion of a long, thick-walled cylinder with inner radius ri of 10 mm

and an outer radius ro of 20 mm, subject to internal pressure p. The material constants

are E=11050 MPa, ν = 0.454 and σvm = 0.5 MPa at zero equivalent plastic strain. In the

plastic range, the material does not harden (perfect plastic behavior). Consequently, the

von Mises stress at the zero equivalent plastic strain applies to the complete plastic range.

A quarter of the cylinder is modeled with three 20-node brick elements in the radial

direction and 5 in the circumferential direction. In the axial direction, only one element

layer is modeled, with its upper and lower layers of nodes fixed in the axial direction (plane

strain assumption). Reduced integration is used throughout. Instead of applying an internal

pressure, the nodes at the inner radius are moved in the radial direction in a uniform way.

The reason for this is shown in Figure 6.4: as soon as the cylinder is fully in the plastic

regime, the internal pressure steadily decreases. Therefore, it cannot be used as the loading

parameter. Also shown is the thickness of the cylinder.

Figure 6.5 shows the change in volume. Notice that during plastic deformation, the

volume decreases slightly. Accordingly, the plastic flow is not completely isochoric. This

is discussed in more detail in Section 6.8.

Comparison with the results published by Simo (Simo 1988b) shows good agreement.

Accordingly, 20-node brick elements with reduced integration can be used for large strain

plasticity. The use of fully integrated 20-node brick elements leads to divergence.

294 FINITE STRAIN ELASTOPLASTICITY

0

0

0.5

1

10 20 30 40 50 60 70 80

0.1

0.2

0.3

0.4

0.6

0.7

0.8

0.9

Cylinder thickness

Internal pressure

ri(mm)

d/10(mm), p/pmax (−)

(Simo 1988b)

Figure 6.4 Variation of the internal pressure and thickness

−0.5

0

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

10 20 30 40 50 60 70 80

ri(mm)

(V0 − V )/V0 (‰)

Figure 6.5 Variation of the volume

Навигация