5.5 Anisotropic Elasticity with a von Mises–type Yield Surface

Back

In the previous section, we introduced the Cailletaud model for single crystals. In order to

use the model, 21 parameters must be determined (3 elastic constants and 9 viscoplastic

constants for each slip system) for the relevant temperature range. This is a huge and

expensive task. Therefore, one frequently resorts to the following approximation: the elastic

range is properly described by the anisotropic elasticity tensor. The yield surface, however,

is assumed to be isotropic of the von Mises form. In this respect, the equations are similar to

the ones in Section 5.3.1. However, because of the anisotropic elastic behavior, the solution

method is more complex and closely linked to the solution procedure in the Cailletaud

model.

5.5.1 Basic equations

The governing equations are merely a concretization of Equations (5.76) to (5.83):

1. Elastic stress–strain relation

σ = C : (_ _p). (5.280)

2. Internal variable relationships

Two internal variables are used: an isotropic scalar variable q1 denoting the radius

of the elastic domain in deviatoric stress space and a kinematic tensor variable q2

denoting its center. The relationship between the internal variables in stress space

{q1, q2} and the corresponding ones in strain space {α1, α2} is assumed to be linear.

A generalization to other functional relationships is no problem.

q1 = d1α1 (5.281)

q2 = 23

d2α2. (5.282)

The factor 23

is introduced such that the equivalent quantities satisfy (cf Equation

(5.10))

q

eq

2

= d2α

eq

2 . (5.283)

3. Yield surface (Equation (5.89))

_dev (σ) + q2_ + _23

(q1 r0) = 0. (5.284)

The parameter r0 is the yield stress at zero-equivalent plastic strain.

INFINITESIMAL STRAIN PLASTICITY 263

4. Evolution equations

In the associative theory, they are derived from the yield surface in the form of

Equations (5.79) and (5.80) and, for a von Mises type surface, Equations (5.95) to

(5.97):

˙_

p = ˙γ n (5.285)

˙α

1 = ˙γ_23

(5.286)

˙α

2 = ˙γ n (5.287)

where

n := ξ

_ξ_ (5.288)

and

ξ := dev (σ ) + q2. (5.289)

5. Kuhn–Tucker equations

γ˙ 0, f (σ, q1, q2) 0, γ˙f (σ, q1, q2) = 0. (5.290)

6. Consistency condition

γ˙ ˙ f (σ, q1, q2) = 0. (5.291)

Viscous effects are taken into account by a Norton-type law

˙_

peq = A(σvm)n (5.292)

or

σvm = g(˙_ peq) = _˙_peq

A

_(1/n)

(5.293)

and Equation (5.284) is replaced by

_dev (σ ) + q2_ + _23

_q1 r0 g(˙_ peq)   = 0. (5.294)

Finally, recall that (Equations (5.112) and (5.113))

α1 = α

eq

2

= _peq = _23

γ. (5.295)

5.5.2 Numerical procedure

Starting from known quantities at time t = tn, the solution at t = tn+1 is what

is being looked for. Using the trial-and-error procedure explained in previous

264 INFINITESIMAL STRAIN PLASTICITY

sections, we first assume that no plasticity takes place in [tn, tn+1]. Consequently

(Equations (5.114)–(5.118)),

_

p

n+1

= _

p

n (5.296)

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

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

γn+1 = γn (5.299)

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

p

n+1). (5.300)

If

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

(q1,n+1 r0) 0 (5.301)

the assumption was right and the solution at t = tn+1 is found. If Equation (5.301) is not

satisfied, the following set of 24 equations in 24 unknowns, obtained by backward Euler

discretization of Equations (5.279) to (5.282), (5.285) to (5.287) and (5.294) has to be

solved:

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

p

n+1) (6 equations) (5.302)

_q1,n+1 = d1_α1,n+1 (1 equation) (5.303)

_q2,n+1 = 23

d2_α2,n+1 (5 equations) (5.304)

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

 

q1,n+1 r0 g

__23

_γn+1

__ = 0 (1 equation)

(5.305)

__

p

n+1

= _γn+1nn+1 (5 equations) (5.306)

_α1,n+1 = _γn+1_23

(1 equation) (5.307)

_α2,n+1 = _γn+1nn+1 (5 equations) (5.308)

in the unknowns _σ n+1 (6), _q1,n+1 (1), _q2,n+1 (5), _γn+1 (1), __

p

n+1 (5), _α1,n+1

(1) and _α2,n+1 (5). Because of the anisotropic character of Equation (5.302), the solution

method of Section 5.3 cannot be used. However, notice the similarity of Equations (5.302)

to (5.308) to Equations (5.191) to (5.195). Indeed, since

nn+1 = σf (σn+1, q1,n+1, q2,n+1) (5.309)

the present set of equations can be considered as a special case of Equations (5.191) to

(5.195) for h = f , q := {q1, q2}, g replaced by

2/3 g and just 1 slip system. Focusing

on the solution method starting at Equation (5.207), one obtains the following residual:

_R

(k)

n+1

_ =



_

p

n+1

+ _

p

n

α1,n+1 + α1,n

α2,n+1 + α2,n



(k)

+ _γ(k)

n+1





nn+1

_23

nn+1





(k)

. (5.310)

INFINITESIMAL STRAIN PLASTICITY 265

For the determination of _A(k)

n+1

_

1

, the second derivatives of f with respect to σ and q

are needed. Using Equations (5.90) to (5.93) one arrives at

2

σσ f =

σ

ξ

_ξ_

= 1

_ξ_

σ

(ξ ) 1

_ξ_2 ξ

σ

_ξ_

= 1

_ξ_

_I 13

I I_ 1

_ξ_2 ξ ξ

_ξ_ (5.311)

= 1

_ξ_

_I 13

I I n n_ (5.312)

and in a similar way,

2

q2q2

f = 2

σq2

f = 2

q2σ f = 2

σσ f = χ

_ξ_ (5.313)

where

χ := I 13

I I n n. (5.314)

Let B be a fourth-order tensor of the form

B := aI + bI I + cn n (5.315)

with a, b, c R,

_n_ = 1 (5.316)

and n is deviatoric:

n : I = 0. (5.317)

Since fourth-order tensors in three-dimensional space can be viewed as 9 × 9 matrices,

we know that tensor contraction (A : B) is associative and that there is a neutral element

I. However, for tensors of the form in Equation (5.315), contraction is also commutative.

Indeed,

(a1I + b1I I + c1n n) : (a2I + b2I I + c2n n)

= a1a2I + (a1b2 + b1a2 + 3b1b2)I I + (a1c2 + c1a2 + c1c2)n n (5.318)

which is symmetric in the indices 1 and 2. Straightforward calculation shows that

B

1 = 1

a

_I b

a + 3b

I I c

a + c

n n_ (5.319)

and

A : χ = χ : A = aχ. (5.320)

266 INFINITESIMAL STRAIN PLASTICITY

Now, Equation (5.208) reduces to the following form:

_A(k)

n+1

_

1

=



C

1

n+1

+ _γn+1

_ξ n+1_χn+1 0

_γn+1

_ξ n+1_χn+1

0

1

d1

0

_γn+1

_ξ n+1_χn+1 0

3

2d2

I + _γn+1

_ξ n+1_χn+1



(k)

. (5.321)

Defining

a := 3

2d2

(5.322)

b :=

_γ(k)

n+1

_ξ n+1_ (5.323)

(notice that b is a function of n and k, although not explicitly indicated!) and dropping the

indices n + 1 and (k) for simplicity, one arrives at

_A      

1 =



C

1 + bχ 0 bχ

0 d

1

1 0

bχ 0 aI + bχ



. (5.324)

In the further derivation, _A     will be needed and it is clearly numerically advantageous if

this inversion can be performed in a largely analytical way. Denoting

_A       :=

P 0 R

0 d1 0

Q 0 S

, (5.325)

the submatrices P, Q, R and S satisfy

C

1 + bχ bχ

bχ aI + bχ

_ : P R

Q S

_ = I 0

0 I

_

. (5.326)

To solve this system, the block in the first row and second column of the left matrix will

be reduced by premultiplying the first block equation by (aI + bχ), the second by bχ and

subtracting the second from the first. This results in

(aI + bχ) : (C

1 + bχ) (bχ) : (bχ) 0

bχ aI + bχ

_ : P R

Q S

_ = aI + bχ

0 I

_

.

(5.327)

Notice that the block in the first row and second column actually reads

(aI + bχ) : [(bχ) : Q] (bχ) : [(aI + bχ) : Q] = 0 (5.328)

which is only true by virtue of the associativity and above all the commutativity of the

tensor-contraction operation for this kind of tensors.

INFINITESIMAL STRAIN PLASTICITY 267

P can be obtained from Equation (5.327) by solving the following equation:

_(aI + bχ) : C

1 + abχ_ : P = aI + bχ. (5.329)

In order to find __γ(k)

n+1, the equivalent of Equation (5.216) has to be solved. Since

_F(k)

n+1

_ = _H(k)

n+1

_ =





nn+1

_23

nn+1





(k)

(5.330)

it is clear that P : n, Q : n, R : n and S : n are needed, and not P, Q, R and S. This is

an easier task to accomplish. Since (Equation (5.327))

[aI + bχ] : Q = [bχ : P] (5.331)

[(aI + bχ) : C

1 + abχ] : R = bχ (5.332)

[aI + bχ] : S = [I bχ : R] (5.333)

and

χ : n = χ : I = 0 (5.334)

[aI + bχ]1 = 1

a + b

I + b

3a

I I + b

a

n n_ (5.335)

one arrives at

[(aI + bχ) : C

1 + abχ] : (P : n) = an (5.336)

Q : n = b

a + b

χ : (P : n) (5.337)

R : n = 0 (5.338)

S : n = 1

a

n. (5.339)

Consequently, only one 6 × 6 set of equations must be solved (Equation (5.336), because

of symmetry conditions, the nine equations reduce to six), the other equations are explicit.

The equivalent of Equation (5.216) now reads







n(k)

n+1

_23

n(k)

n+1





T

: _A

(k)

n+1

_ :





n(k)

n+1

_23

n(k)

n+1





+ _23

_γ g(k)

n+1



__γ(k)

n+1

= _

f (k)

n+1

_23

g(k)

n+1

_ _R

(k)

n+1

_T

: _A

(k)

n+1

_T

:





n(k)

n+1

_23

n(k)

n+1





. (5.340)

268 INFINITESIMAL STRAIN PLASTICITY

Since

_A

(k)

n+1

_ :





n(k)

n+1

_23

n(k)

n+1





:=





P : n

_23

d1

Q : n + S : n





(k)

n+1

(5.341)

and

_A

(k)

n+1

_T

= _A

(k)

n+1

_ (5.342)

(_A

(k)

n+1

_

1

is symmetric, Equation (5.324), and the inverse of a symmetric matrix is also

symmetric), one obtains

 

(n : P : n)(k)

n+1

+ 23

d1 + 23

d2 + _23

_γ g(k)

n+1

_

__γ(k)

n+1

= _

f (k)

n+1

_23

g(k)

n+1

_ _R

(k)

n+1

_T

:





P : n

_23

d1

Q : n + S : n





(k)

n+1

(5.343)

where

g

(k)

n+1

=

__23

_γ(k)

n+1

A_t

_

1

n

(5.344)

and

_γ g(k)

n+1

= _23

1

An_t

__23_γ

(k)

n+1

A_t

_

1

n

1

(5.345)

represent the viscous effects. Equation (5.343) is a linear equation in __γ(k)

n+1. Once the

correction to the consistency parameter is known, the corrections to the internal variables

can be calculated using the following equivalent of Equation (5.269) (again dropping the

indices n + 1 and (k) for simplicity):

 

__p

_α1

_α2



 =



C

1 0 0

0 d

1

1 0

0 0 aI



:



P 0 R

0 d1 0

Q 0 S



:





R_

Rα1

Rα2





(5.346)

=





C

1 : (P : R_ + R : Rα2 )

Rα1

a(Q : R_ + S : Rα2 )





(5.347)

INFINITESIMAL STRAIN PLASTICITY 269

where



R_

Rα1

Rα2



= 'R( + __γ





n

_23

n





(5.348)

is the update of the residual. Substitution of Equations (5.329) and (5.332) into the first

block equation of Equation (5.347) yields

[(aI + bχ) : C

1 + abχ] : [C : '__p(] = [(aI + bχ) : R_ bχ : Rα2 ]. (5.349)

Notice that the left-hand matrix of Equation (5.349) is the same as in Equation (5.336)

and consequently the LU decomposition (in an upper and lower matrix (Zienkiewicz

and Taylor 1989)) can be reused. Furthermore, C

1, needed to obtain '__p( from C :

'__p(, was already computed to obtain the left-hand side in Equation (5.336). Substituting

Equations (5.331) and (5.333) into the lower block equation leads to

'_α2( = a[aI + bχ]1 : _'Rα2( bχ : 2P : 'R_( + R : 'Rα2(3    . (5.350)

Using the first block equation of Equation (5.347), this is equivalent to

'_α2( = a[aI + bχ]1 : _'Rα2( bχ : 'C : __p( (5.351)

or

'_α2( = a

a + b

I + b

3a

I I + b

a

n n_ : 'Rα2( ab

a + b

χ : 'C : __p( . (5.352)

Accordingly, reintroducing the indices (C is no function of n), one finds for the corrections

of the internal variables

_(aI + bχ(k)

n+1) : C

1 + abχ(k)

n+1

_ : _C : '__p((k)

n+1

_

= _(aI + bχ) : 'R_((k)

n+1

bχ(k)

n+1 : 'Rα2((k)

n+1

_ (5.353)

and

'_α2((k)

n+1

= a

a + b

I + b

3a

I I + b

a

n(k)

n+1

n(k)

n+1

_ : 'Rα2((k)

n+1

ab

a + b

χ(k)

n+1 : _C : '__p((k)

n+1

_ . (5.354)

Recall that b is also a function of k and n. As soon as all the corrections are determined,

the satisfaction of the flow rule

____

f (k+1)

n+1

_23

g(k+1)

n+1

____

< TOL (5.355)

270 INFINITESIMAL STRAIN PLASTICITY

can be checked. If satisfied, convergence is reached and the loop can be left. If not, a

new correction must be determined. Notice that this loop corresponds to the outer loop

in Section 5.4. There is no inner loop since we deal with single surface plasticity. Once

convergence is reached, the consistent elastoplastic tangent matrix can be determined. This

is obtained from the equations equivalent to Equations (5.276) and (5.278):

'dσ n+1( = Pn+1 : _I G

1

n+1[nn+1 Pn+1 : nn+1]_ : 'd_n+1( (5.356)

and

C

ep

n+1

= Pn+1 G

1

n+1(Pn+1 : nn+1) (Pn+1 : nn+1). (5.357)

Here, P : n and P are both needed. For the calculation of P, Equation (5.329) can be used.

Notice that P is needed only after convergence is reached. The quantity Gn+1 is defined

by (Equation (5.343)):

Gn+1 = (n : P : n)n+1 + 23

d1 + 23

d2 + _23

_γ gn+1. (5.358)

Notice that in the absence of kinematic hardening, d2 can be zero. In that case, a is undetermined

(Equation (5.322)). Hence, care must be taken in the implementation to express

the equations in terms of a

1. For instance, Equation (5.336) then reads

_(I + ba

1χ) : C

1 + bχ_ : (P : n) = n. (5.359)

Summarizing, the algorithm runs as follows:

1. Compute the elastic predictor and the value of the yield surface (Equations (5.296)–

(5.300))

2. Check for plasticity (Equation (5.301)). If satisfied, the solution is found. Else, go

to (3).

3. Loop construct

(a) Calculate the residuals of the flow rule, Equation (5.305), and evolution laws,

Equation (5.310). If small enough, exit.

(b) Calculate a correction to the consistency parameter, Equation (5.343).

(c) Calculate a correction to the internal variables, Equations (5.353) and (5.354);

go to 3a.

4. Determine the consistent elastoplastic tangent matrix, Equation (5.357).

5.5.3 Special case: isotropic elasticity

For isotropic materials, the above equations can be substantially simplified. Indeed,

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

INFINITESIMAL STRAIN PLASTICITY 271

where μ and λ are Lamґe’s constants. Using Equation (5.319), one obtains

C

1 = _ 1

9K

1

6μ

_I I + 1

2μ

I (5.361)

where

K := λ + 23

μ. (5.362)

Defining

α := 1

2μ

(5.363)

β := 1

9K

1

6μ

(5.364)

one obtains

C

1 = αI + βI I . (5.365)

From Equation (5.329), the following expression for P results:

P = _(aI + bχ) : C

1 + abχ_

1

: [aI + bχ]. (5.366)

Substitution of Equation (5.365) into (5.366) and taking into account the laws applicable

to tensors of the type at stake (such as Equation (5.319)), one arrives after some algebra at

P = 1

[aα + (a + α)b]

 

(a + b)I + ab 3β(a + b)

3(α + 3β)

I I + ab

α

n n_ (5.367)

and

P : n = 1

α n = 2μn. (5.368)

Accordingly, the coefficient of __γ(k)

n+1 in Equation (5.343) reduces to

Gn+1 = 2μ + 23

d1 + 23

d2 + _23

_γ g(k)

n+1 (5.369)

which is identical to the corresponding coefficient in Equation (5.141) since

_γ g(k)

n+1

= _23

__peqg(k)

n+1. (5.370)

The equivalent consistent elastoplastic tangent matrix takes the form (Equation (5.357)

C

ep

n+1

= Pn+1 4μ2nn+1 nn+1

2μ + 23

d1 + 23

d2 + 23

__pgn+1

. (5.371)

272 INFINITESIMAL STRAIN PLASTICITY

To check whether Equation (5.371) coincides with Equation (5.162), α and β are substituted

into Pn+1 (the index n + 1 is dropped for simplicity):

P = 1

1 + (2μ + a1)b

[2μI + λI I ]

+ 2μb

1 + (2μ + a1)b

_a

1I + K(1 3βa

1)I I + 2μn n_ (5.372)

which is equivalent to

P = (2μI + λI I ) (2μ)2b

1 + (2μ + a1)b

[I 13

I I n n]. (5.373)

Recall that a = 23

d2 and b = _γ/_ξ_. Notice that b contains _ξ_, whereas Equation (5.162)

contains _ξ trial_. The connection between both is given by Equation (5.134), in which the

term _h

eq

2 takes the form

_h

eq

2

= d2__peq = d2_23

_γ (5.374)

in the present context of linear hardening laws. Consequently, Equation (5.134) leads to

_ξ trial_ = _ξ_[1 + b(2μ + a

1)] (5.375)

and Equation (5.373) yields

P = (2μI + λI I ) (2μ)2_γ

_ξ trial_

_I 13

I I n n_ . (5.376)

Equations (5.371) and (5.376) reproduce Equation (5.162). The isotropic case is recovered

as a special case of the anisotropic formulation.