5.2 The General Framework of Plasticity

Back

5.2.1 Theoretical derivation

Throughout the present chapter, it will be assumed that the strains and rotations are

small, such that the infinitesimal quantities ˜eKL and ˜EKL can be used. One recalls

The Finite Element Method for Three-dimensional Thermomechanical Applications Guido Dhondt

2004 John Wiley & Sons, Ltd ISBN: 0-470-85752-8

226 INFINITESIMAL STRAIN PLASTICITY

(Section (1.14.4)),

˜E

KL ˜eKLδk

Kδl

L (5.1)

SKL σ klδ K

k δ L

l . (5.2)

Furthermore, it will be assumed that a rectangular spatial coordinate system and a rectangular

material coordinate system are used and that both coincide. To emphasize this, the

infinitesimal strain tensor will be represented by the new symbol _kl .

How can we characterize a viscoplastic material? A simple stress–strain test on steel

yields Figure 5.1. One notices that the material is elastic for small stresses. For growing

stress (σ > σA), the curve deviates from the elastic straight line by an amount that will be

called the plastic strain _p. It is this amount that is not recovered after unloading (point

C). When loading again, the material remains elastic up to σ = σB before accumulating

further plastic strain. Accordingly, the elastic range depends on the previous plastic flow.

From these considerations, the following assumptions (some of which are valid for the

infinitesimal range only) seem plausible:

1. The total strain can be decomposed in an additive way into an elastic part and a

plastic part:

_ = _e + _p. (5.3)

2. Plastic flow starts as soon as the stress (axial stress in a tensile specimen, an appropriate

stress invariant in the two- and three-dimensional case) reaches a specific

value, which can change as a function of the previous amount of plastic deformation.

Accordingly, this value is a function and will be denoted q. The amount of plastic

flow is represented by α. Hence,

q = q(α). (5.4)

A

B

C

σ

_

_e _p

Figure 5.1 Elastoplastic stress–strain curve

INFINITESIMAL STRAIN PLASTICITY 227

q is an internal dynamic variable (comparable to a stress), α is the kinematic counterpart

(comparable to a strain). The material remains elastic if

σ < q σ +q < 0. (5.5)

The function

f (σ, q) = σ + q = 0 (5.6)

represents the boundary of the elastic range and is called the yield surface.

In theory, plastic deformation takes place if the yield surface is exceeded. However,

because of the plastic flow, the yield surface is expanded such that the stress stays on the

yield surface (in the absence of viscosity). This can be made obvious by simply performing

an unloading–loading experiment, as shown in Figure 5.1. Consequently, the increasing

stress drags the yield surface along, and Equation (5.6) keeps its validity during the plastic

flow. The reverse statement is not true: Equation (5.6) does not necessarily imply plastic

deformation. Imagine we load until σ = q and freeze the loading at that point: this is a

purely elastic process. Summarizing, during plastic deformation, one has

σ = E(_ _p) (5.7)

q = q(α) (5.8)

σ + q = 0. (5.9)

Assume that we know the total strain and want to know the stress. Equations (5.7) to

(5.9) yield three equations in the four unknowns σ , _p, q and α: one equation is lacking.

Physically, we do not know at what rate the plastic flow is accumulated: an evolution

equation for _p (and actually, also for α) is lacking. This equation will be obtained by

maximizing the plastic dissipation.

Recall from Chapter 1 that the entropy rate was obtained by substituting an expression

for the free energy, Equation (1.390), into the entropy inequality, yielding Equation (1.392)

and subsequently Equation (1.396). For our infinitesimal (one-dimensional) considerations,

Equation (1.411) reduces to

ψ = ψ(_, θ,θ,X). (5.10)

In the previous analysis, it was shown that the total strain of the elastic theory should

be replaced by _e. Furthermore, a dependence of ψ on the internal kinematic variable is

assumed. Hence,

ψ = ψ(_ _p, α, θ,θ,X). (5.11)

Note that the time derivative of the internal kinematic variable α is not necessarily continuous:

at unloading, the time rate of α discontinuously drops to zero. This also applies to

_p. Taking the time derivative yields

˙ψ

= ψ

(_ _p)

˙_

ψ

(_ _p)

˙_

p + ψ

α

˙α

+ ψ

θ

˙ θ + ψ

θ

∇˙θ. (5.12)

228 INFINITESIMAL STRAIN PLASTICITY

Substituting this expression into Equation (1.389) yields (ρ0 ρ)

1

θ

_ρ

ψ

_e

+ σ

_ ˙_ + ρ

θ

ψ

_e

˙_

p ρ

θ

ψ

α

˙α

ρ

θ

_ψ

θ

+ η

_ ˙ θ ρ

θ

ψ

θ

∇˙θ 1

θ2 qθθ 0

(5.13)

where qθ stands for the heat conduction (the superscript θ is introduced in this chapter to

avoid confusion between the heat conduction qθ and the internal dynamic variable q). This

inequality is satisfied if

σ = ρ

ψ

_e (5.14)

η = ψ

θ

(5.15)

ψ

θ

= 0 (5.16)

and

σ˙_p ρ

ψ

α

˙α

1

θ

qθθ 0. (5.17)

The first two terms do not vanish because of the time history of plastic deformation. As in

Equation (5.14), q is now defined by

q = ρ

ψ

α

(5.18)

(or, the other way around, Equation (5.18) can be viewed as the definition of ψ(α) through

Equation (5.4)). Inequality (5.17) is satisfied if

σ˙_p + q˙α 0 (5.19)

and

1

θ

qθθ 0. (5.20)

Now, the evolution equations for ˙_ p and ˙α are obtained by postulating that for a given ˙_p

and ˙α, the state (σ, q) will prevail, which maximizes (Simo and Hughes 1997) (Halphen

and Nguyen Quoc Son 1975)

dp = σ˙_ p + q˙α (5.21)

subject to

f (σ, q) 0. (5.22)

This postulate is accepted to hold for metals, but does not necessarily hold for all kinds of

material. This also amounts to minimizing

lp := σ˙_p q˙α + ˙γf(σ, q) (5.23)

INFINITESIMAL STRAIN PLASTICITY 229

with respect to σ and q where

γ˙ 0 (5.24)

and subject to

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

(Luenberger 1989). By taking the derivative of lp with respect to σ and q, one obtains the

evolution equations

˙_

p = ˙γ

f (σ, q)

σ

(5.26)

˙α

= ˙γ

f (σ, q)

q

. (5.27)

Equations (5.22), (5.24) and (5.25) are called the Kuhn–Tucker conditions. Although there

are two new equations, Equations (5.26) and (5.27), there is also a new unknown γ˙ , called

the consistency parameter (sometimes called the plastic rate parameter). By Equation

(5.24), γ˙ cannot be negative, and by Equation (5.25), it can be strictly positive only if

yielding takes place. The notation γ˙ was chosen because it is easiest to consider it as a

rate of accumulated plastic flow. To emphasize that, for plastic deformation to occur, the

stress state has to persist on the yield surface the following equation is added

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

This is also called the consistency condition. Summarizing, the following equations apply:

1. Elastic stress–strain relations

σ = E(_ _p) (5.29)

2. Internal variable relationship

q = h(α) (5.30)

3. Yield surface

f (σ, q) = 0 (5.31)

4. Evolution equations

˙_

p = ˙γ

f (σ, q)

σ

(5.32)

˙α

= ˙γ

f (σ, q)

q

(5.33)

5. Kuhn–Tucker equations

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

230 INFINITESIMAL STRAIN PLASTICITY

6. Consistency condition

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

The evolution equations are also called the flow rule (Equation (5.32)) and the hardening

law (Equation (5.33)). Equations (5.29) to (5.35) are quite general and are easy to extend

to higher dimensions. Notice particularly that the maximum-dissipation principle implies

that the evolution equations can be derived from the yield surface. This kind of model is

called associative.

Returning to our basic yield surface, Equation (5.9), extending it to negative stresses

in the form

|σ| + q 0 (5.36)

and substituting this equation into (5.32) and (5.33) yields

˙_

p = ˙γ sgn(σ ) (5.37)

˙α

= ˙γ . (5.38)

Equations (5.37) and (5.38) reveal that, for our simple example, ˙ γ and ˙α are the magnitude

of the plastic strain rate:

˙ γ =_

_

˙_

p_

_ . (5.39)

Hence, since σ = h(α) and α = _p for monotonic loading (˙α = ˙_p and let α = 0 for _p =

0), one can derive h(α) from Figure 5.1 by subtracting the elastic strain at a given stress

level (Figure 5.2).

The type of hardening considered so far is called isotropic hardening, since it equally

applies to positive and negative stresses (Figure 5.3). In practice, the Bauschinger effect

is frequently observed: after plastic deformation in the tensile range, plastic deformation

in the compressive range takes place at higher stress levels than expected. Pure kinematic

hardening implies that the size of the elastic range has not changed, just its origin: BD =

2 OA in Figure 5.3. Introducing an internal variable q2 for the center of the elastic range,

this amounts to a yield surface of the form

|σ + q2| + c 0 (5.40)

h(α)

_p, α

Figure 5.2 Isotropic hardening curve

INFINITESIMAL STRAIN PLASTICITY 231

O

A

B

C

E

σ

_

D Kinematic hardening

F Isotropic hardening

G Combined isotropic/Kinematic hardening

Figure 5.3 Types of hardening

where c is a constant. If BF = 2 BE, pure isotropic hardening applies. Point G symbolizes

a state in between, that is, combined hardening. For this general case, the yield surface

looks like

|σ + q2| + q1 0 (5.41)

leading to the flow rule

˙_

p = ˙γ sgn(σ + q2) (5.42)

and the evolution equations

˙α

1 = ˙γ (5.43)

˙α

2 = ˙γ sgn(σ + q2). (5.44)

From Equation (5.42), one again notices that γ˙ is the rate of the accumulated plastic

strain (in absolute value). The only functions left to be determined are q1(α1) and q2(α2),

both of which are material characteristics to be obtained by experiments. Notice that the

inclusion of kinematic hardening does not substantially change the governing relations

(Equations (5.29)–(5.35)): instead of q, one now deals with q1 and q2 or, equivalently, a

two component vector q.

The theory derived so far can also be extended to include viscous effects. For many

materials, one observes that during loading the stress state exceeds the yield surface and

slowly creeps back with time until the yield surface is reached from above. In Figure 5.4,

232 INFINITESIMAL STRAIN PLASTICITY

σ

_

t = t1

t = t2

t = 0

Figure 5.4 Viscoplasticity concept

the external loading (forces, temperature etc.) is frozen at t = t1. Because of the rate

at which the loading was applied, the material did not have time to accumulate enough

plastic deformation to attain equilibrium conditions. Physically, this means that since plastic

flow is largely equivalent to dislocation motion, the dislocations did not have enough

time to redistribute in accordance with the applied load. As time goes by, an equilibrium

configuration is reached. This phenomenon is also called creep and is usually modeled by

a creep law of the form

˙_

p = g

1(σ 0) (5.45)

where σ 0 is the overload. This is the stress amount by which the yield surface is exceeded

and which constitutes the driving force for the creep process. Consequently, if plastic

deformation occurs, Equation (5.31) has to be replaced by

f (σ, q) = g(˙_p). (5.46)

An example of a creep law is the Norton law, which is of the form

˙_

p = A(σ0)n. (5.47)

The parameter A is usually a small number (depending on the unit of σ 0), and n > 1, for

example, n = 5. Both A and n are material parameters.

5.2.2 Numerical implementation

Ultimately, Equations (5.29) to (5.35) have to be solved; for viscous problems, Equation

(5.31) is to be replaced by Equation (5.46). Notice that Equations (5.34) and (5.35)

have a kind of regulating character: they characterize the discrete split into elastic deformation

and plastic deformation.

The usual finite-element procedure consists of the creation of a tangent-stiffness matrix

at an instantaneous deformation field and solving the resulting linear equations to obtain

INFINITESIMAL STRAIN PLASTICITY 233

a new deformation field, after which the procedure can be repeated until convergence.

Thus, two major tasks ensue: the calculation of a consistent tangent-stiffness matrix and

the calculation of the stress corresponding to a given deformation field in order to check

convergence (satisfaction of the equilibrium conditions). In both cases, the deformation field

is given. Consequently, the total strain is known. This is an extremely important fact and

the starting point of most numerical algorithms in present finite-element implementations.

Accordingly, _n+1 is known, where t = tn+1 is the new time step and all quantities at t = tn

are known. Further strategy consists of trial and error. One first assumes that no plasticity

occurs in the present step, that is,

_

p,trial

n+1

= _

p

n (5.48)

qtrial

n+1

= qn (5.49)

γ trial

n+1

= γn. (5.50)

Hence (Equation (5.29)),

σ trial

n+1

= E(_n+1 _

p

n). (5.51)

If the new stress state lies within the elastic range, that is, if

f (σtrial

n+1, qtrial

n+1) 0 (5.52)

then the solution is found

_

p

n+1

= _

p,trial

n+1 (5.53)

qn+1 = qtrial

n+1 (5.54)

γn+1 = γ trial

n+1 (5.55)

σn+1 = σ trial

n+1 (5.56)

and the tangent modulus is the elastic one. If this is not the case, plasticity takes place and

the following equations must be satisfied at t = tn+1:

f (σ, q) = g(˙_ E

1 ˙σ) (5.57)

˙_

E

1˙σ = ˙γ

f (σ, q)

σ

(5.58)

˙

h1(q) = ˙γ

f (σ, q)

q

. (5.59)

These are (for the one-dimensional case) three equations in the three unknowns σ , q and

γ˙ . Now, a backward Euler scheme is applied to turn the differential equations (5.57) to

(5.59) into difference equations. To this end, the equations are evaluated at t = tn+1 and

the first derivative ˙_ is replaced by

˙_

= _n+1 _n

_t

+ O(_t) (5.60)

234 INFINITESIMAL STRAIN PLASTICITY

and similarly for the other first derivatives. One can prove that the backward Euler scheme

is first-order accurate and unconditionally stable. Applying this to Equations (5.57) to (5.59)

leads to

f (σn+1, qn+1) = g[__n+1 E

1(σn+1 σn)] (5.61)

__n+1 E

1(σn+1 σn) = _γn+1

f (σn+1, qn+1)

σ

(5.62)

h

1(qn+1) h

1(qn) = _γn+1

f (σn+1, qn+1)

q

(5.63)

where

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

_γn+1 := γn+1 γn. (5.65)

Equations (5.61) to (5.63) are three nonlinear equations in three unknowns σn+1, qn+1

and _γn+1. They can be solved using the customary mathematical techniques to solve

sets of nonlinear equations. To concretize the further derivation, the yield surface of

Equation (5.36) is taken, no creep is assumed and a linear isotropic hardening law is

chosen of the form

q = σ0 Kα (5.66)

where σ0 and K are constants. Consequently,

|σn+1| + qn+1 = 0 (5.67)

__n+1 E

1(σn+1 σn) = _γn+1sgn(σn+1) (5.68)

(σ0 qn+1) (σ0 qn) = K_γn+1. (5.69)

The sign of σn+1 is the same as that of σ trial

n+1 (Figure 5.5).

sgn(σn+1) = sgn(σ trial

n+1). (5.70)

Accordingly (Equations (5.68) and (5.69)),

σn+1 = E__n+1 E_γn+1sgn(σ trial

n+1) + σn (5.71)

qn+1 = qn K_γn+1 (5.72)

yielding (Equation (5.67))

_γn+1 = 1

E + K

_E__n+1sgn(σ trial

n+1) + σnsgn(σ trial

n+1) + qn_ . (5.73)

Substitution of _γn+1 into Equations (5.71) and (5.72) yields σn+1 and qn+1. In particular

(Equation (5.71)),

σn+1 = E(_n+1 _n) E2(_n+1 _n)

E + K

Eσn

E + K

Eqnsgn(σ trial

n+1)

E + K

+ σn (5.74)

INFINITESIMAL STRAIN PLASTICITY 235

σ

_

σ trial

n+1

σn+1

_n

σn

_n+1

__n+1

_γn+1

Figure 5.5 Trial-and-error method

from which the plastic tangent is obtained:

dσn+1

d_n+1

= EK

E + K

. (5.75)

Because of the particular choice of q and f , the resulting equation, Equation (5.73), was

linear and easy to solve.