5.3 Three-dimensional Single Surface Viscoplasticity

Back

5.3.1 Theoretical derivation

The governing equations for three-dimensional applications are the same as for onedimensional

applications and their derivation is similar (cf Equations (5.29)–(5.35)):

1. Elastic stress–strain relations

σ = _

_e (5.76)

2. Internal variable relationships

q = h(α) (5.77)

3. Yield surface

f (σ, q) = 0 (5.78)

4. Evolution equations

˙_

p = ˙γ

f (σ, q)

σ

(5.79)

˙α

= ˙γ

f (σ, q)

q

(5.80)

236 INFINITESIMAL STRAIN PLASTICITY

5. Kuhn–Tucker equations

γ˙ 0, f (σ, q) 0, γ˙f (σ, q) = 0 (5.81)

6. Consistency condition

γ˙ ˙ f (σ , q) = 0. (5.82)

If viscous effects are to be taken into account, Equation (5.78) is replaced by

f (σ, q) = g(˙_p) (5.83)

and f > 0 is feasible. Although the form of the equations is similar, due to the tensorial

character of the quantities involved, one arrives at a much larger set of nonlinear equations

than in the one-dimensional case. Indeed, Equations (5.78) to (5.80) lead to a set of nonlinear

equations the size of which is the sum of the number of independent stress components

and the number of internal variables plus one. Here, an example covering the usual isotropic

metal plasticity will be given.

To fully determine the plasticity model defined by Equations (5.76) to (5.83), one has

to define the yield surface, choose the internal variables and possibly define a creep law.

In three dimensions, the yield surface is not so trivial as in the one-dimensional case.

Indeed, one has to find one scalar equation connecting the tensorial quantities σ and q.

Furthermore, for isotropic materials, the yield surface should contain invariants only, that

is, I1σ , I2σ and I3σ . Practical observations have shown that the hydrostatic pressure p does

not significantly lead to plasticity in metals (does not apply to soils!). Since

p := 13

I1σ (5.84)

the first invariant does not enter the yield condition explicitly. Therefore, a new stress

tensor, the deviatoric stress s is defined by

s := dev σ := σ + pI . (5.85)

On the basis of the deviatoric stress, a new invariant is defined called the von Mises stress

σvm by

σvm := _32

_s_2 = _32

_dev σ _2. (5.86)

Since

_s_2 = _σ_2 13

(I1σ )2

= 23

(I1σ )2 2I2σ (5.87)

one finds

σvm = _(I1σ )2 3I2σ . (5.88)

INFINITESIMAL STRAIN PLASTICITY 237

The von Mises stress is a measure of the shear energy. The factor

2/3 is introduced such

that the von Mises stress in a tensile test coincides with the applied tensile stress. Including

isotropic and kinematic hardening (internal variables q1 and q2 respectively), the following

yield surface is proposed (Huber–von Mises yield surface):

_dev (σ) + q2_ + _23

q1 = 0. (5.89)

In deviatoric space (axes s1, s2 and s3), Equation (5.89) is a sphere with radius

2/3 q1 and

center q2 (Figure 5.6). The inside of the sphere is the elastic range. During plasticity, the

sphere can both expand (isotropic hardening) and move (kinematic hardening). The internal

variable q1 is a scalar, whereas q2 is a tensor. Since only the deviatoric part of the stress

is relevant in Equation (5.89), the hydrostatic part of q2 remains arbitrary throughout the

analysis.

Defining

ξ := dev (σ) + q2 (5.90)

and since

q2 = dev (q2) (5.91)

ξ

_ξ_ = ξ

_ξ_ (5.92)

σ

(dev (σ)) = I 13

I I (5.93)

where

(I)

ij

kl := 12

(δi

kδ

j

l

+ δi

lδ

j

k) (5.94)

is the fourth-order identity tensor and I is the second-order identity tensor, straightforward

application of Equations (5.79) and (5.80) yields

˙_

p = ˙γ

ξ

_ξ_ (5.95)

s1

s2

s3

_23

q1

q2

Figure 5.6 Yield surface in deviatoric stress space

238 INFINITESIMAL STRAIN PLASTICITY

˙α

1 = ˙γ_23

(5.96)

˙α

2 = ˙γ

ξ

_ξ_. (5.97)

Similar to the von Mises stress, an equivalent plastic strain is defined by

_peq := _23

__p_. (5.98)

The factor

2/3 is introduced such that in a tensile test, the equivalent plastic strain is

equal to the plastic strain in the tensile direction. Indeed, since a hydrostatic pressure does

not lead to plasticity, the plastic deformation is volume preserving, and hence,

_

p

11

+ _

p

22

+ _

p

33

= 0. (5.99)

Consequently, for uniaxial plastic strain, one has

_

p

22

= _

p

33

= 12

_

p

11 (5.100)

and

_peq = _23

_(_

p

11)2 + (_

p

22)2 + (_

p

33)2 = _

p

11. (5.101)

From Equations (5.95) and (5.98), it is apparent that the physical meaning of

2/3 γ˙ is

the equivalent plastic strain rate.

Finally, the relationships q1(α1) and q2(α2) are left to be defined. The variable q1

means the von Mises stress with respect to the reference stress q2 at yield (cf Equation (5.89)),

α1 is the accumulated plastic strain. This relationship must be obtained from experiments

and will be written as

q1 = h1(α1). (5.102)

The second set of internal variables q2 is tensorial and the relationship

q2 = h2(α2) (5.103)

is more difficult to obtain. Time differentiation of Equation (5.103) yields

˙q2 = h2

α2

: ˙α2 (5.104)

which implies that the tensor ˙q2 is not necessarily parallel to ˙α2. This complicates the subsequent

analysis. Since the material is isotropic, it seems convenient to assume that the kinematic

hardening is also isotropic, that is, we write the relationships for Equations (5.103)

and (5.104) for the equivalent properties:

q

eq

2

= h

eq

2 (α

eq

2 ) (5.105)

leading to

˙q

eq

2

= h

eq

2

α

eq

2

˙α

eq

2 (5.106)

INFINITESIMAL STRAIN PLASTICITY 239

or

_32

_˙q2_ = h

eq

2

α

eq

2

_23

_ ˙α2_. (5.107)

Equation (5.107) suggests the following isotropic tensorial relationship:

˙q

2 = 23

h

eq

2

α

eq

2

˙α

2. (5.108)

Comparison with Equation (5.104) leads to

h2

α2

= 23

h

eq

2

α

eq

2

II . (5.109)

One finally obtains

˙q

1 = −˙h1 (5.110)

˙q

2 = 23

h

eq

2

α

eq

2

γ˙

ξ

_ξ_

= _23

˙h

eq

2

ξ

_ξ_ (5.111)

since

˙α

eq

2

= _23

γ˙ (5.112)

and Equation (5.96). From Equations (5.95) to (5.97), it is obvious that

_peq = α1 = α

eq

2 . (5.113)

5.3.2 Numerical procedure

Just as in the one-dimensional case, the total strain is assumed to be given, all quantities are

known at t = tn, they are to be determined at t = tn+1. Again the trial-and-error procedure

is used. At first, it is assumed that no plasticity occurs:

_

p

n+1

= _

p

n (5.114)

q1,n+1 = q1,n (5.115)

q2,n+1 = q2,n (5.116)

γn+1 = γn (5.117)

σ n+1 = _

_e

____

n+1

. (5.118)

If

_dev (σn+1) + q2,n+1_ + _23

q1,n+1 0 (5.119)

240 INFINITESIMAL STRAIN PLASTICITY

the solution is found. Else, the following equations have to be solved at t = tn+1:

σ = _

_e (6 equations) (5.120)

_dev + q2_ + _23

q1 = 23

g(˙_peq) (1 equation) (5.121)

˙_

p = ˙γ

ξ

_ξ_ (5 equations) (5.122)

˙q

1 = −˙h1 (1 equation) (5.123)

˙q

2 = _23

˙h

eq

2

ξ

_ξ_ (5 equations). (5.124)

These are 18 equations in 18 unknowns: σ (6), _p (5), q1 (1), q2 (5), γ˙ (1) (recall that _p

and q2 are deviatoric in nature). If we assume that the material is isotropic in the elastic

regime, the equations can be further simplified. Indeed, the elastic stress–strain relationship

for a linear elastic isotropic material satisfies

σ = λtr(_e)I + 2μ_e. (5.125)

Hence,

s = dev σ = 2μdev _e. (5.126)

Equation (5.122) shows that ˙_ p is deviatoric, consequently,

˙_

p = dev (˙_p) = dev (˙_ ˙_e) = dev (˙_) dev (˙_e) (5.127)

and

2μdev (˙_ ) s˙ = 2μγ˙

ξ

_ξ_. (5.128)

Replacing the time derivatives by backward Euler differences and defining

__n+1 := _n+1 _n (5.129)

and similarly for the other expressions, one obtains

2μdev (__n+1) sn+1 + sn = 2μ_γn+1

ξ n+1

_ξ n+1_. (5.130)

Now, strial

n+1 is obtained from sn by assuming that __ is purely elastic, that is,

strial

n+1

= sn + 2μdev (__n+1) (5.131)

which leads to

strial

n+1

sn+1 = 2μ_γn+1

ξ n+1

_ξ n+1_ (5.132)

INFINITESIMAL STRAIN PLASTICITY 241

for Equation (5.130). Backward Euler for Equation (5.124) yields

q2,n+1 q2,n

= _23

_h

eq

2 (_

peq

n+1) h

eq

2 (_

peq

n )       

ξ n+1

_ξ n+1_. (5.133)

Subtracting Equation (5.132) from Equation (5.133), one gets

ξ n+1 ξ trial

n+1

= 2μ_γn+1 _23

_h

eq

2,n+1

_ ξ n+1

_ξ n+1_ (5.134)

where

ξ trial

n+1 := strial

n+1

+ q2,n. (5.135)

Equation (5.134) shows that the vectors ξ n+1 and ξ trial

n+1 are parallel (therefore, the algorithm

is sometimes called the radial return method). This result is crucial in the present derivation.

If the kinematic hardening had not been isotropic, this simplification would not apply! Since

all terms in Equation (5.134) are parallel, the equation applies to their size equally well:

_ξ n+1_ _ξ trial

n+1

_ = 2μ_γn+1 _23

_h

eq

2,n+1. (5.136)

There is only one equation left to be satisfied: the yield condition, which reads

_ξ n+1_ + _23

q1,n+1 = _23

g(__

peq

n+1). (5.137)

One finds

h1,n+1 = h1(_

peq

n + __

peq

n+1) = h1

_

_

peq

n + _23

_γn+1

_

. (5.138)

Finally, one gets for the yield condition

_ξ trial

n+1

_ _2μ_γn+1 + _23

 

h

eq

2

_

_

peq

n + _23

_γn+1

_ h

eq

2 (_

peq

n )

_

_23

h1

_

_

peq

n + _23

_γn+1

_ = _23

g

__23

_γn+1

_

. (5.139)

Consequently, we finally arrive at one nonlinear equation in _γn+1. This equation can be

solved using a Newton–Raphson technique. Denoting the initial value for the unknown

_γn+1 by _γ(0)

n+1 and writing

_γ(k+1)

n+1

= _γ(k)

n+1

+ __γ(k)

n+1 (5.140)

linearization of Equation (5.139) about _γ(k)

n+1 yields

_ξ trial

n+1

_ _2μ_γ(k)

n+1

+ _23

 

h

eq

2

_

_

peq

n + _23

_γ(k)

n+1

_ h

eq

2 (_

peq

n )

_

_23

h1

_

_

peq

n + _23

_γ(k)

n+1

_ _23

g

__23

_γ(k)

n+1

_

2μ + 23

(h1 + h

eq

2 )

_peq

_____

_

peq

n +_23

_γ

(k)

n+1

+ 23

g

__peq

____

_23

_γ

(k)

n+1

__γ

(k)

n+1

= 0. (5.141)

242 INFINITESIMAL STRAIN PLASTICITY

Once _γn+1 is known, one finds for the other variables

_

peq

n+1

= _

peq

n + _23

_γn+1 (5.142)

q1,n+1 = h1(_

peq

n+1) (5.143)

q2,n+1 = q2,n

_23

_h

eq

2 (_

peq

n+1) h

eq

2 (_

peq

n )       

ξ trial

n+1

_ξ trial

n+1

_

(5.144)

_

p

n+1

= _

p

n + _γn+1

ξ trial

n+1

_ξ trial

n+1

_

(5.145)

σ n+1 = _

_e

____

n+1

. (5.146)

5.3.3 Determination of the consistent elastoplastic tangent matrix

The consistent elastoplastic tangent matrix is the instantaneous slope of the stress–total

strain relationship. The term “consistent” points to the fact that the slope has to be determined

for the actual numerical scheme used, that is, it depends on the numerical procedure

(using another scheme, e.g. the midpoint rule instead of backward Euler, will lead to

another slope). It is well known (Simo and Hughes 1997) that the slope derived for the

present numerical scheme deviates from the continuum tangent. Consistency of the slope

is a prerequisite for the quadratic convergence of the Newton–Raphson scheme.

For materials that are linear in the elastic range, Equation (5.146) reduces to

σ n+1 = C : _e

n+1

= C : (_n+1 _

p

n+1) (5.147)

where

C := 2_

_e_e . (5.148)

Substituting Equation (5.145), one now arrives at the following stress–strain relationship:

σ n+1 = C : (_n+1 _

p

n _γn+1nn+1) (5.149)

where

nn+1 := ξ n+1

_ξ n+1_

=

ξ trial

n+1

_ξ trial

n+1

_

. (5.150)

Thus, the tangent relation at t = tn+1 takes the form

dσ n+1 = C : (d_n+1 d_γn+1nn+1 _γn+1 dnn+1). (5.151)

For an isotropic elastic material, C amounts to

C = 2μI + λI I (5.152)

INFINITESIMAL STRAIN PLASTICITY 243

where λ, μ are Lamґe’s constants. Since n is deviatoric I : n = 0, one can write

C : n = 2μn (5.153)

and Equation (5.151) reduces to

dσn+1 = C : d_n+1 2μ(d_γn+1nn+1 + _γn+1 dnn+1) (5.154)

or

dσn+1 = C 2μ

_nn+1 _γn+1

_n+1

+ _γn+1

nn+1

_n+1

__ : d_n+1. (5.155)

The consistent elastoplastic tangent is the expression in square braces. To determine

_γn+1/_n+1, Equation (5.139) is differentiated with respect to _n+1:

_ξ trial

n+1

_

_n+1

_2μ + 23

_

peq

n+1

h

eq

2

+ 23

_

peq

n+1

h1 + 23

__

peq

n+1

g_ _γn+1

_n+1

= 0. (5.156)

Since

_ξ trial

n+1

_

_n+1

=

_ξ trial

n+1

_

ξ trial

n+1

:

ξ trial

n+1

_n+1

(5.157)

and (combining Equation (5.135) with Equation (5.131))

ξ trial

n+1

= sn + 2μdev (_n+1 _n) + q2,n (5.158)

one obtains (Equations (5.92) and (5.93))

_ξ trial

n+1

_

_n+1

= nn+1 : 2μ_I 13

I I_ = 2μnn+1. (5.159)

Consequently,

_γn+1

_n+1

= _2μ + 23

_

peq

n+1

h

eq

2

+ 23

_

peq

n+1

h1 + 23

__

peq

n+1

g_

1

2μnn+1. (5.160)

The derivative in the last term of Equation (5.155) yields

nn+1

_n+1

= nn+1

ξ trial

n+1

:

ξ trial

n+1

_n+1

=

_ 1

_ξ trial

n+1

_

I 1

_ξ trial

n+1

_2

ξ trial

n+1

ξ trial

n+1

_ξ trial

n+1

_

_

: 2μ_I 13

I I_

= 2μ

_ξ trial

n+1

_

(I nn+1 nn+1) : _I 13

I I_

= 2μ

_ξ trial

n+1

_

_I 13

I I nn+1 nn+1_ . (5.161)

244 INFINITESIMAL STRAIN PLASTICITY

Substituting Equations (5.160) and (5.161) into Equation (5.155) finally yields

Cep = C (2μ)2_γn+1

_ξ trial

n+1

_

_I 13

I I nn+1 nn+1_

(2μ)2nn+1 nn+1 _2μ + 23

_

peq

n+1

h

eq

2

+ 23

_

peq

n+1

h1 + 23

__

peq

n+1

g_

1

. (5.162)

It can be shown that for _γn+1 = 0, the continuum elastoplastic tangent is obtained (Simo

and Hughes 1997). Accordingly, because of the finite size of the increments, the consistent

numerical tangent deviates from the continuum tangent.