5.4 Three-dimensional Multisurface Viscoplasticity: the Cailletaud Single Crystal Model

Back

Single crystals are advanced metallic materials consisting of just one crystal. Substantial

progress made in the last two decades in casting technology enables the manufacturers to

control crystal growth in the liquid metal by carefully monitoring the cooling conditions.

Thus, a highly anisotropic material ensues, in contrast with the usual metallic materials

(polycrystals), in which the different orientations of the many crystals assure the isotropic

properties. In this section, the focus will be on nickel-base alloys, exhibiting a face cube

centered (FCC) crystal structure. A good reference on crystalline plasticity is (Havner

1992).

5.4.1 Theoretical considerations

In single crystals, viscoplasticity is mainly due to a crystallographic dislocation slip. Other

mechanisms will not be considered here. The crystallographic slip planes and directions

for nickel-base superalloys at high temperature are known (Mґeric et al. 1991) and can

be divided into octahedral slip systems (12 systems consisting of 4 {111} planes with 3

< 011 > directions per plane, Figure 5.7) and cubic slip systems (6 systems consisting of

3 {001} planes with 2 < 011 > directions per plane, Figure 5.8).

Accordingly, the slip planes and directions are explicitly known and are generally

denoted by their normal nβ and unit vector lβ respectively. Here, β stands for any of the

18 slip systems. For each slip system, an orientation tensor mβ is defined by

mβ := 12

(nβ lβ + lβ nβ). (5.163)

Whether dislocations move along a slip system basically depends on the shear stress component

τ β in the slip direction:

τ β = mβ : σ = (nβ lβ) : σ = nβT σ lβ. (5.164)

This is really a one-dimensional system and we can fall back on yield-surface formulations

such as in Equation (5.41):

f β(σ , q) = |τ β + q

β

2

| + q

β

1

0. (5.165)

INFINITESIMAL STRAIN PLASTICITY 245

Figure 5.7 Octahedral slip systems

Figure 5.8 Cubic slip systems

246 INFINITESIMAL STRAIN PLASTICITY

However, now we are dealing with 18 yield surfaces at the same time, which may intersect

each other at so-called corner points. This is an example of multisurface viscoplasticity.

The underlying theory has been developed in the 1950s and 1960s (Koiter 1960), and is a

straightforward extension of the one-dimensional derivation in Section 5.2. The governing

equations for m slip systems are

1. Elastic stress–strain relations

σ = _

_e (5.166)

2. Internal variable relationships

q = h(α) (5.167)

3. Yield surfaces

f β(σ , q) = 0 (5.168)

4. Evolution equations

˙_

p =

m

_

β=1

γ˙ β f β(σ, q)

σ

(5.169)

˙α

=

m

_

β=1

γ˙ β f β(σ, q)

q

(5.170)

5. Kuhn–Tucker conditions

γ˙ β 0, fβ(σ, q) 0, γ˙ βf β(σ, q) = 0 (5.171)

6. Consistency conditions

γ˙ β ˙ f β (σ, q) = 0. (5.172)

q is a vector of internal variables with a size that is usually a multiple of the number

of slip systems. In Equations (5.169) and (5.170), γ˙ β is only nonzero for the active slip

systems. In what follows, we will concentrate on a particular single crystal viscoplastic

model developed by Georges Cailletaud and coworkers, (Mґeric et al. 1991), (Mґeric and

Cailletaud 1991). For other single crystal models, see (Fedelich 2002) and (Meissonnier

et al. 2001). The Cailletaud model does not completely fit into the theory described by

Equations (5.166) to (5.171) because of the following two aspects:

1. The model is not associative, which means that the evolution equations are derived

from a function that differs from the yield surface.

2. The evolution equation, Equation (5.170), is modified by a term depending on the

total accumulated plasticity.

INFINITESIMAL STRAIN PLASTICITY 247

Denoting the yield surfaces by hβ (not to confuse with h in Equation (5.167)), the extra

term for ˙α by w and the creep term (viscous term) by gβ(γ˙ β), one obtains

1. Elastic stress–strain relations

σ = _

_e (5.173)

2. Internal variable relationships

q = h(α) (5.174)

3. Yield surfaces

hβ (σ, q) = gβ(γ˙ β) (5.175)

4. Evolution equations

˙_

p =

m

_

β=1

γ˙ β f β(σ, q)

σ

(5.176)

˙α

=

m

_

β=1

γ˙ β f β(σ, q)

q

+ wβ __ t

0

γ˙ β dt

__ (5.177)

5. Kuhn–Tucker conditions

γ˙ β 0, hβ (σ, q) 0, γβhβ(σ , q) = 0 (5.178)

6. Consistency conditions

γ˙ β˙hβ(σ, q) = 0. (5.179)

Specifically, in the Cailletaud model, there are 2m internal dynamic variables, which

will be denoted by q

β

1 , q

β

2 , β = 1, . . . ,m. The internal variable relationships take the form

q

β

1

= bβQβα

β

1 (5.180)

q

β

2

= cβα

β

2 (5.181)

where bβ , Qβ and cβ, β = 1, . . . ,m are constants. Consequently, Equations (5.180) and

(5.181) are linear. The yield surfaces are defined by

hβ (σ, q) :=_

_

_

τ β + q

β

2

___

r

β

0

+

m

_

α=1

Hβαqα

1 . (5.182)

The parameters r

β

0 are the initial yield values and Hβα are the interaction coefficients

between the slip systems. Equation (5.182) is a slightly more complicated form than

Equation (5.165). The potential function for the evolution equations reads

f β (σ, q) :=_

_

_

τ β + q

β

2

___

+ q

β

1

+ dβ

2cβ

(q

β

2 )2 + 1

2Qβ

(q

β

1 )2. (5.183)

248 INFINITESIMAL STRAIN PLASTICITY

This is the yield function in Equation (5.165), augmented by quadratic terms in q

β

1 and q

β

2 .

The parameters dβ, β = 1, . . . ,m are constants, cβ and Qβ already appeared in the internal

variable relationships. The only functions left are gβ and wβ. They will be specified in the

next section.

5.4.2 Numerical aspects

The numerical procedure to solve Equations (5.173) to (5.179) is totally similar to the

methods treated in the previous sections. The two basic considerations in the analysis are

1. The procedure is strain-driven, that is, we start from a given increment __n+1 =

_n+1 _n and look for the corresponding stress and tangent-stiffness matrix.

2. A trial-and-error procedure is used, starting from the assumption that the step is

purely elastic. A verification of the yield condition tells us whether this assumption

is right.

Consequently, we assume in step n + 1

αn+1 = αn (5.184)

_

p

n+1

= _p (5.185)

γ

β

n+1

= γ β

n, β= 1, . . . ,m. (5.186)

The stress is obtained from Equation (5.173). Now, the yield surfaces are verified. If

hβ (σ, q) 0, β {1, . . . ,m} (5.187)

then the step is elastic and the solution is found. Equation (5.173) yields the stress, the

tangent-stiffness matrix is the elasticity tensor. If

B(0)

act := {β|hβ (σ, q) > 0} _= (5.188)

plastic flow takes place. B(0)

act is the initial set of active slip planes. Equations (5.173) and

(5.174) lead to

˙σ

= C : (˙_ ˙_p) (5.189)

˙q

= h

α

: ˙α =: D : ˙α. (5.190)

For the Cailletaud model, D is a constant matrix. Using a backward Euler scheme, Equations

(5.189) and (5.190) and Equations (5.175) to (5.177) can be rewritten as

_σ n+1 = Cn+1 : (__n+1 __

p

n+1) (5.191)

_qn+1 = Dn+1 : _αn+1 (5.192)

hβ (σn+1, qn+1) = gβ(_γ

β

n+1), β B(0)

act (5.193)

__

p

n+1

= _

βB

(0)

act

_γ

β

n+1σ f β(σn+1, qn+1) (5.194)

INFINITESIMAL STRAIN PLASTICITY 249

_αn+1 = _

βB

(0)

act

_γ

β

n+1qf β(σn+1, qn+1). (5.195)

The abbreviations

σ f β := f β

σ

(5.196)

and likewise for the derivative with respect to q were used, and the functions wβ were

dropped for now. Notice that only the active slip planes are considered in Equations (5.193)

to (5.195)! Substituting Equations (5.191) and (5.192) into Equations (5.193) to (5.195) one

finally obtains

hβ (σn+1, qn+1) = gβ(_γ

β

n+1), β B(0)

act (5.197)

__n+1 C

1

n+1 : _σ n+1 = _

βB

(0)

act

_γ

β

n+1σ f β (σn+1, qn+1) (5.198)

D

1

n+1 : _qn+1 = _

βB

(0)

act

_γ

β

n+1qf β (σn+1, qn+1). (5.199)

If mact is the number of active slip planes, Equation (5.197) represents mact equations,

Equation (5.198) represents 6 equations and Equation (5.199) represents 2 × mact equations

in the unknowns σn+1 (6), qn+1 (2 × mact) and _γ

β

n+1 (mact). Hence, we obtain 3 × mact +

6 equations in 3 × mact + 6 unknowns. For the inactive slip planes, Equations (5.184) to

(5.186) apply. Equations (5.197) to (5.199) are the basis for our further consideration.

5.4.3 Stress update algorithm

The stress can be determined by solving Equations (5.197) to (5.199) for _σn+1. Since

these equations are nonlinear, a Newton–Raphson iterative technique is used for their solution

(cf Section 3.1). Assume that we have an intermediate solution denoted by a superscript

(k). To obtain a better approximation, the Equations (5.197) to (5.199) are linearized at the

solution (k) and solved. Denoting

h

β

n+1 := hβ (σn+1, qn+1) (5.200)

and similarly for f and g, linearization yields

_h

β(k)

n+1

g(k)

n+1

_ + σ h

β(k)

n+1 : _σ (k)

n+1

+ qh

β(k)

n+1 : _q(k)

n+1

_γ g

β(k)

n+1__γ

β(k)

n+1

= 0 (5.201)



_

p

n+1

+ _

p

n + _

βB

(k)

act

_γ

β

n+1σ f

β

n+1



(k)

+ C

1(k)

n+1 : _σ (k)

n+1

+ _

βB

(k)

act

__γ

β(k)

n+1 σ f

β(k)

n+1

+ _

βB

(k)

act

_γ

β(k)

n+1

_2

σσ f

β

n+1 : _σ n+1 + 2

σqf

β

n+1 : _qn+1_(k)

= 0 (5.202)

250 INFINITESIMAL STRAIN PLASTICITY



αn+1 + αn + _

βB

(k)

act

_γ

β

n+1qf

β

n+1



(k)

+ D

1(k)

n+1 : _q(k)

n+1

+ _

βB

(k)

act

__γ

β(k)

n+1 qf

β(k)

n+1

+ _

βB

(k)

act

_γ

β(k)

n+1

_2

qσ f

β

n+1 : _σ n+1 + 2

qqf

β

n+1 : _qn+1_(k)

= 0 (5.203)

where

_σ (k)

n+1 := σ (k+1)

n+1

σ (k)

n+1 (5.204)

_q(k)

n+1 := q(k+1)

n+1

q(k)

n+1 (5.205)

__γ

β(k)

n+1 := _γ

β(k+1)

n+1

_γ

β(k)

n+1 . (5.206)

Notice that in the above derivation, __n+1 as well as all quantities with the superscript (k)

are assumed to be given (the process is strain-driven). The first terms in square brackets

in each equation are the function values of Equations (5.193) to (5.195), equivalent to

f (x0) F in Equation (3.2). If the equations are satisfied, these function values should be

zero. Therefore, they are also called the residual. The other terms in the equations are the

gradients. Equations (5.201) to (5.203) are linear in __γ

β(k)

n+1 , _σ (k)

n+1 and _q(k)

n+1. Defining

_R(k)

n+1

_ := __

p

n+1

+ _

p

n

αn+1 + αn

(k)

+ _

βB

(k)

act

_γ

β(k)

n+1

_σ f

β

n+1

qf

β

n+1

_(k)

(5.207)

which is the residual of Equations (5.202) to (5.203),

_A

(k)

n+1

_

1

:=

C

1

n+1

+_βB

(k)

act

_γ

β

n+12

σσ f

β

n+1 _βB

(k)

act

_γ

β

n+12

σqf

β

n+1

_βB

(k)

act

_γ

β

n+12

qσ f

β

n+1 D

1

n+1

+_βB

(k)

act

_γ

β

n+12

qqf

β

n+1

(k)

(5.208)

and finally

_F

β(k)

n+1

_ :=



σ f

β

n+1

qf

β

n+1



(k)

(5.209)

_H

β(k)

n+1

_ :=



σ h

β

n+1

qh

β

n+1



(k)

(5.210)

INFINITESIMAL STRAIN PLASTICITY 251

then Equations (5.202) and (5.203) can be written as

_R(k)

n+1

_ + _A(k)

n+1

_

1

:



_σ (k)

n+1

_q(k)

n+1



+ _

βB

(k)

act

__γ

β(k)

n+1

_F

β(k)

n+1

_ = 0 (5.211)

which is equivalent to

_A(k)

n+1

_ : _R(k)

n+1

_ +



_σ (k)

n+1

_q(k)

n+1



+ _

βB

(k)

act

__γ

β(k)

n+1

_A(k)

n+1

_ : _F

β(k)

n+1

_ = 0. (5.212)

From Equation (5.201), one gets

_h

β(k)

n+1

g

β(k)

n+1

_ + _H

β(k)

n+1

_T

:



_σ (k)

n+1

_q(k)

n+1



_γ g

β(k)

n+1__γ

β(k)

n+1

= 0. (5.213)

Premultiplying Equation (5.212) by _H

α(k)

n+1

_T

and inserting Equation (5.213) leads to

_H

α(k)

n+1

_T

: _A

(k)

n+1

_ : _R

(k)

n+1

_ + _γ gα(k)

n+1__γα(k)

n+1

_hα(k)

n+1

gα(k)

n+1

_

+ _

βB

(k)

act

__γ

β(k)

n+1

_H

α(k)

n+1

_T

: _A

(k)

n+1

_ : _F

β(k)

n+1

_ = 0, α B(k)

act . (5.214)

Defining

(Gαβ)(k)

n+1 := _H

α(k)

n+1

_T

: _A

(k)

n+1

_ : _F

β(k)

n+1

_ + _γ g

β(k)

n+1 δαβ (5.215)

Equation (5.214) can be rewritten as

_

βB

(k)

act

(Gαβ )(k)

n+1__γ

β(k)

n+1

= _hα(k)

n+1

gα(k)

n+1

_ _H

α(k)

n+1

_T

: _A

(k)

n+1

_ : _R

(k)

n+1

_ , α B(k)

act .

(5.216)

These are m(k)

act linear equations in m(k)

act unknowns. Their solution yields __γ

β(k)

n+1 , β B(k)

act .

Substituting into Equation (5.212) yields an expression for _σ (k)

n+1 and _q(k)

n+1 and iteration

(k) seems to be finished. However, there is one further consideration to be taken into

account. Equation (5.188), which defines the active slip systems, is not completely correct

in the sense that it constitutes a necessary condition to be an active system but not a

sufficient one. Indeed, because of the presence of corner points in the yield surface, the

consistency parameter after iteration k + 1 (Equation (5.206)),

_γ (k+1)

n+1 := _γ(k)

n+1

+ __γ(k)

n+1 (5.217)

is not necessarily positive. For details the reader is referred to (Simo and Hughes 1997). All

active planes for which _γ (k+1)

n+1

0 have to be removed from B(k)

act and Equation (5.216)

has to be solved again until for all active slip systems _γ (k+1)

n+1 > 0. Accordingly, the

252 INFINITESIMAL STRAIN PLASTICITY

number of active slip systems can decrease from iteration to iteration, which is symbolized

by the superscript (k) on B(k)

act .

What form do the above equations take in the Cailletaud model? Recall that the potential

for the evolution equations is defined by

f β(σ, q) := |σ : mβ + q

β

2

| + q

β

1

+ dβ

2cβ

(q

β

2 )2 + 1

2Qβ

(q

β

1 )2. (5.218)

Consequently,

σ f

β

n+1

= mβsgn(τ

β

n+1

+ q

β

2,n+1) (5.219)

qf

β

n+1

=





0

...

0

1 + q

β

1

Qβ

sgn(τ β + q

β

2 ) + dβ

cβ q

β

2

0

...

0





row(2β 1)

row(2β)

row(2m)

. (5.220)

In Equation (5.177), qf

β

n+1 was modified by a function wβ taking the total plasticity into

account. In the Cailletaud model, wβ is defined by

wβ :=





0

...

0

(ϕβ 1)sgn(τ β + q

β

2 )

0

...

0





row(2β) (5.221)

where

ϕβ := φβ + (1 φβ)eδβ & t

0

γ˙ β dt . (5.222)

The parameters φβ and δβ are material constants. In the numerical procedure, the accumulated

plasticity in step n + 1 can be approximated by

_ t

0

γ˙ β dt

n

_

i=1

_γ

β

i

+ _γ

β(k)

n+1 . (5.223)

INFINITESIMAL STRAIN PLASTICITY 253

In our derivation, the effect of wβ will be incorporated into a modified qf

β

n+1:

qf

β

n+1

=





0

...

0

1 + q

β

1

Qβ

ϕβsgn(τ β + q

β

2 ) + dβ

cβ q

β

2

0

...

0





row(2β 1)

row(2β)

row(2m)

. (5.224)

This modified value will be used instead of the original value in all previously derived

formulas. Notice that theoretically we now have

qf

β

n+1(σn+1, qn+1,_γ

β

n+1) (5.225)

that is, qf

β

n+1 is not only a function of σn+1 and qn+1, but also of _γ

β

n+1. Accordingly,

the linearization of Equation (5.199) is not completely correct any more. This effect will

be neglected. It will at most slow down convergence. The second derivatives yield

2

σσ f

β

n+1

= 0 (6 ×6 matrix) (5.226)

2

σqf

β

n+1

= 0 (6 × 2m matrix) (5.227)

2

qσ f

β

n+1

= 0 (2m ×6 matrix) (5.228)

2

qqf

β

n+1

=



0 0 0

0

1

Qβ 0

0 dβ

cβ

0

0 0 0



(2m × 2m matrix) (5.229)

where the submatrix in Equation (5.229) occupies rows and columns (2β 1) and 2β.

From Equations (5.180), (5.181) and (5.190) one finds

D = Diag(b1Q1, c1, b2Q2, . . . , bmQm, cm). (5.230)

254 INFINITESIMAL STRAIN PLASTICITY

Consequently (Equation (5.173)),

_F

β(k)

n+1

_ =





mβsgn(τ

β

n+1

+ q

β

2,n+1)

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

0

...

0

1 + q

β

1,n+1

Qβ

ϕβ sgn(τ

β

n+1

+ q

β

2,n+1) + dβ

cβ q

β

2,n+1

0

...

0





(k)

row(2β + 5)

row(2β + 6)

(5.231)

and (Equation (5.209))

_A

(k)

n+1

_

1

= Diag_

C

1

n+1

...

1

b1Q1

+ _γ1

n+1

1

Q1 ,

1

c1

+ _γ1

n+1

d1

c1 ,

1

b2Q2

+ _γ2

n+1

1

Q2 , . . .

_

(k)

.

(5.232)

The sum over B(k)

act in Equation (5.208) was replaced by a sum over all slip systems since

_γn+1 = 0 for an inactive slip system. Combining Equations (5.231) and (5.232), one

obtains

_A

(k)

n+1

_ : _F

β(k)

n+1

_ =





Cn+1 : mβ sgn(τ

β

n+1

+ q

β

2,n+1)

. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

0

...

0

Qβ+q

β

1

1

bβ

+_γ

β

n+1

ϕβ cβ sgn(τ

β

n+1

+q

β

2,n+1)+dβ q

β

2,n+1

1+_γ

β

n+1dβ

0

...

0





(k)

row(2β + 5)

row(2β + 6)

. (5.233)

Now recall the expression for the yield surface in the Cailletaud model:

hβ (σ, q) := |σ : mβ + q

β

2

| r

β

0

+

m

_

α=1

Hβαqα

______________1 . (5.234)

INFINITESIMAL STRAIN PLASTICITY 255

Differentiation yields

σ h

β

n+1

= mβsgn(τ

β

n+1

+ q

β

2,n+1) (5.235)

qδ

1

h

β

n+1

=

m

_

α=1

Hβαδαδ = Hβδ (5.236)

qδ

2

h

β

n+1

= sgn(τ

β

n+1

+ q

β

2,n+1) δβδ. (5.237)

In matrix form,

_H

β(k)

n+1

_ =





mβsgn(τ

β

n+1

+ q

β

2,n+1)

. . . . . . . . . . . . . . . . . . . .

Hβ1

0

Hβ2

0

...

Hββ

sgn(τ

β

n+1

+ q

β

2,n+1)

Hβ(β+1)

...

Hβm

0





(k)

. (5.238)

Hence (Equation (5.215)),

(Gββ)(k)

n+1

= mβ : Cn+1 : mβ + Hββ

_Qβ + q

β(k)

1,n+1

_

_ 1

bβ

+ _γ

β(k)

n+1

_

+

_ϕβcβ + dβq

β(k)

2,n+1sgn(τ

β(k)

n+1

+ q

β(k)

2,n+1)_

_1 + _γ

β(k)

n+1 dβ_

+ _γ g

β(k)

n+1 (5.239)

(Gαβ)(k)

n+1

= mα : Cn+1 : mβsgn(τ

β(k)

n+1

+ q

β(k)

2,n+1)sgn(τ α(k)

n+1

+ qα(k)

2,n+1)

+ Hαβ

_Qβ + q

β(k)

1,n+1

_

_ 1

bβ

+ _γ

β(k)

n+1

_

, α_= β. (5.240)

256 INFINITESIMAL STRAIN PLASTICITY

For the right-hand side of Equation (5.216), the following term is needed:

'T α

n+1( := __H

α(k)

n+1

_T

: _A

(k)

n+1

__T

=





Cn+1 : mαsgn(τ α

n+1

+ qα

2,n+1)

. . . . . . . . . . . . . . . . . . . . . . . . . . .

Hα1Q1

(1/b1)+_γ1

n+1

0

Hα2Q2

(1/b2)+_γ2

n+1

0

...

HααQα

(1/bα)+_γα

n+1

cαsgn(τ α

n+1

+qα

2,n+1)

1+_γα

n+1dα

...

HαmQm

(1/bm)+_γm

n+1

0





(k)

. (5.241)

The only function left to be specified is the creep function g. In the Cailletaud model,

the relationship between the viscous shear strain rate and the shear stress in a slip plane is

defined by

___

˙_

pβ

τ _

__

= ) τ β

Kβ

*nβ

(5.242)

where < x >= x for x 0 and < x >= 0 for x < 0. Kβ and nβ are material constants.

The total viscous strain rate is related to the slip plane shear strain rates by (Koiter 1960)

˙_

p =

m

_

β=1

˙_

pβ

τ mβ . (5.243)

Now, combining Equations (5.169) and (5.219) yields

˙_

p =

m

_

β=1

γ˙ βmβsgn(τ

β

n+1

+ q

β

2,n+1). (5.244)

Comparison of Equations (5.243) and (5.244) yields

γ˙ β =__

_

˙_

pβ

τ __

_

(5.245)

since the applied shear stress and resulting shear strain rate have the same sign. Hence,

Equation (5.242) can be written as

γ˙ β = ) τ β

Kβ

*nβ

. (5.246)

INFINITESIMAL STRAIN PLASTICITY 257

In the viscoplastic theory, the stress gβ by which the yield surface is exceeded is to be

relaxed by creep, that is,

< τβ >= gβ. (5.247)

Consequently,

gβ(γ˙ β) = Kβ (γ˙ β)(1/nβ ) (5.248)

= Kβ __γβ

_t

_(1/nβ )

(5.249)

and

_γ g

β(k)

n+1

= Kβ

nβ_t

__γ

β(k)

n+1

_t

_

1

nβ

1

. (5.250)

The last expression is used in Equation (5.239).

Summarizing, one arrives at the following algorithm to obtain σ n+1 from σn:

1. Compute the elastic predictor and the value of the yield surfaces

_

p,trial

n+1

= _

p

n (5.251)

qtrial

n+1

= qn (5.252)

_γtrial

n+1

= 0 (5.253)

σtrial

n+1

= Cn : (_n+1 _

p,trial

n+1 ) (5.254)

h

β,trial

n+1

=_

__

mβ : σ trial

n+1

+ qtrial

2,n+1

___

r

β

0

+

m

_

α=1

Hβαqα,trial

1,n+1. (5.255)

Notice that gβ(_γ trial

n+1) = 0.

2. Check for plasticity.

If h

β,trial

n+1

0, β: step n + 1 is elastic, that is,

σn+1 = σtrial

n+1 (5.256)

the values at t = tn+1 are the trial values: the solution is found.

else

B(0)

act = _β {1, . . . ,m}|h

β,trial

n+1 > 0_ (5.257)

_

p(0)

n+1

= _

p

n (5.258)

α(0)

n+1

= αn (5.259)

_γ

β(0)

n+1

= 0, β= 1, 2, . . . ,m. (5.260)

endif

258 INFINITESIMAL STRAIN PLASTICITY

3. Start of the outer loop

calculate σ(k)

n+1 and q(k)

n+1 from _

p(k)

n+1 and α(k)

n+1 and check if the flow rule and the

evolution equations are satisfied.

σ(k)

n+1

= Cn+1 : (_n+1 _

p(k)

n+1) (5.261)

q(k)

n+1

= Dn+1 : α(k)

n+1 (5.262)

_R(k)

n+1

_ =

__

p

n+1

+ _

p

n

αn+1 + αn

_(k)

+ _

βB

(k)

act

_F

β(k)

n+1

_ (5.263)

h

β(k)

n+1

g

β(k)

n+1

=__

_

mβ : σ(k)

n+1

+ q

β(k)

2,n+1

___

r

β

0

+

m

_

α=1

Hβαqα(k)

1,n+1

Kβ __γ

β(k)

n+1

_t

_

1

nβ

(5.264)

if _

__

h

β(k)

n+1

g

β(k)

n+1

___

< TOL and +

++

R(k)

n+1

+++

< TOL, leave the outer loop.

4. Start of the inner loop

determine __γ

β(k)

n+1 by Equation (5.216) without creep effects.

If _γ

β(k+1)

n+1 := _γ

β(k)

n+1

+ __γ

β(k)

n+1 > 0, β B(k)

act then

recalculate __γ

β(k)

n+1 by Equation (5.216) with creep effects and exit inner

loop.

else

remove the inactive slip planes from B(k)

act and reiterate the inner loop.

endif

End of the inner loop

5. Update the internal variables

Determine

__σ(k)

n+1

_q(k)

n+1

_

from Equation (5.212).

___

p

n+1

_αn+1

_(k)

=

,C

1

n+1 0

0 D

1

n+1

-(k)

:

__σ (k)

n+1

_q(k)

n+1

_

(5.265)

_

p(k+1)

n+1

= _

p(k)

n+1

+ __

p(k)

n+1 (5.266)

α(k+1)

n+1

= α(k)

n+1

+ _α(k)

n+1 (5.267)

_γ

β(k+1)

n+1

= _γ

β(k)

n+1

+ __γ

β(k)

n+1 (5.268)

set k k + 1 and reiterate the outer loop.

INFINITESIMAL STRAIN PLASTICITY 259

6. End of outer loop. Now, the determination of the plastic tangent modulus can start.

Finally, two more remarks:

(a) It is advantageous to substitute Equation (5.265) directly into Equation (5.212)

yielding

___

p

n+1

_αn+1

_(k)

=

,C

1

n+1 0

0 D

1

n+1

-(k)

: _A

(k)

n+1

_ :

: __R(k)

n+1

_ +_βB

(k)

act

_F

β(k)

n+1

___γ

β(k)

n+1

_ (5.269)

where (Equation (5.232)),

C

1 0

0 D

1

_(k)

n+1

: _A

(k)

n+1

_

= Diag

_

I

...

1

1 + b1_γ1

n+1

,

1

1 + d1_γ1

n+1

,

1

1 + b2_γ2

n+1

, . . .

_(k)

. (5.270)

(b) The inner loop is necessary to determine the active slip planes (cf (Simo and Hughes

1997) for more details). In the determination process, the viscous terms are dropped

to make sure that the viscous procedure converges in the limit to the same point on

the yield surface as the inviscid formulation.

5.4.4 Determination of the consistent elastoplastic tangent matrix

The determination of the consistent elastoplastic moduli also starts from Equations (5.197)

to (5.199). We have attained equilibrium for t = tn+1, that is, Equations (5.197) to (5.199)

are identically satisfied and we would like to know how σ changes if _ is perturbed.

Therefore, we differentiate these equations:

σ h

β

n+1 : dσn+1 + qh

β

n+1 : dqn+1 _γ g

β

n+1

d_γ

β

n+1

= 0 (5.271)

d_n+1 + C

1

n+1 : dσn+1 + _

βB

(k)

act

d_γ

β

n+1σ f

β

n+1

+ _

βB

(k)

act

_γ

β

n+1

_2

σσ f

β

n+1 : dσn+1 + 2

σqf

β

n+1 : dqn+1_ = 0 (5.272)

D

1

n+1 : dqn+1 + _

βB

(k)

act

d_γ

β

n+1qf

β

n+1

+ _

βB

(k)

act

_γ

β

n+1

_2

qσ f

β

n+1 : dσn+1 + 2

qqf

β

n+1 : dqn+1_ = 0. (5.273)

260 INFINITESIMAL STRAIN PLASTICITY

These equations are very similar to Equations (5.201) to (5.203). In fact, by replacing

_R

(k)

n+1

_ by _d_n+1

0

, hα(k)

n+1

gα(k)

n+1 by 0, _ by d and dropping the superindex (k), they

are identical. Hence, by comparing with Equation (5.216), one arrives at the following set

of equations:

_

βBact

(Gαβ )n+1 d_γ

β

n+1

= 'T α

n+1(T :

_d_n+1

0

_

, α Bact (5.274)

yielding

d_γ

β

n+1

=

_

αBact

(G

1βα)n+1 'T α

n+1(T

:

_d_n+1

0

_

, β Bact. (5.275)

The equivalent equation of Equation (5.212) reads

_dσn+1

dqn+1

 = _An+1         : __d_n+1

0

 _βBact

_F

β

n+1

_ dγ

β

n+1

 

= _An+1          : _I        _

βBact

_

αBact

(G

1βα)n+1 _F

β

n+1

_ 'T α

n+1(T_ : _d_n+1

0

 

. (5.276)

Only the relationship between dσn+1 and d_n+1 is needed, hence (Equations (5.231),

(5.232) and (5.241)),

'dσn+1( = [Cn+1 _

βBact

_

αBact

Cn+1 : mβsgn(τ

β

n+1

+ q

β

2,n+1)(G

1βα)n+1

mα : Cn+1sgn(τ α

n+1

+ qα

2,n+1)] : 'd_n+1( . (5.277)

Hence, the consistent elastoplastic stiffness matrix C

ep

n+1 satisfies

C

ep

n+1

= Cn+1 _

βBact

_

αBact

(G

1βα)n+1Mβ MαT (5.278)

where

Mα := Cn+1 : mαsgn(τ α

n+1

+ qα

2,n+1). (5.279)

Equation (5.278) shows that each active slip plane modifies the stiffness matrix (without

plastic flow, the consistent elastoplastic stiffness matrix reduces to the elasticity matrix).

Since _G          is not necessarily symmetric (cf Equation (5.240)), Cep is not necessarily symmetric

either. In practice, the matrix is often symmetrized by adding the transpose and

dividing by two.

5.4.5 Tensile test on an anisotropic material

Consider the tensile specimen in Figure 5.9. The axis of the specimen coincides with the

z-axis. The orientation of the anisotropic material is defined by the x

_-y

_-z

_ axis system.

INFINITESIMAL STRAIN PLASTICITY 261

x

x

F

F

h

8 h h

θ

y

z

z

x

_

y, y

_

z

_

Figure 5.9 Geometry of the tensile specimen

0

0

0.5

1

1.5

2

2.5

3

3.5

10 20 30 40 50 60 70 80 90

θ()

γ β ()

Octahedral slip system

Cubic slip system

Figure 5.10 Accumulated plastic slip

The y- and y

_-axes coincide, whereas the z- and z

_-axis include an angle θ. A constant force

F is applied at the end of the specimen in the z-direction. Now we look at what happens

if we vary the angle θ from 0 to 90. In particular, we investigate the accumulated plastic

slip in two different slip systems: the first slip system is octahedral and is characterized

by n = (1, 1, 1) and l = (1, 0,1), the second is a cubic slip system and is defined by

n = (0, 0, 1) and l = (1,1, 0). The slip systems are defined in the local x

_-y

_-z

_ system.

Figure 5.10 shows that the octahedral slip system is activated if the global axes and

the material axes are aligned. Then, the slip direction that is considered includes an angle

262 INFINITESIMAL STRAIN PLASTICITY

of 45 with the loading axis, leading to a large slip. The cubic slip system is activated

especially for angles close to θ = 45. Here again, the angle between the slip direction and

the loading direction is maximized.