On the numerical treatment and analysis of finite deformation ductile single crystal plasticity

Computer methods in applied me~hanics and engineerlnD EI.~FVIER

('omput. Methods Appl. Mech. Engrg. 129 (1996) 235-254

On the numerical treatment and analysis of finite deformation ductile single crystal plasticity Paul Steinmann, Erwin Stein' h~stitt~t ft~r Baumet'h,znik m~d Numerische Mechanik. UJliversitiit ttamtover. Appc/strafle 9A, 30167 Hannover, Germany

Received 18 April 1994:reviscd 3 July 1995

Abstract This contribution is concerned with the numerical analysis and simulation of instability phenomena such as shearband localization in plastic flow pr,~cesscs within the framework of ductile single (fee) crystal viscoplasticity. To this end. a family of implicit cc~nstitttlivc :dg~)rithms and their implemenlaticm are considered in detail. Issues of accuracy a~d conservation of plastic incompressibility arc examined. Moreover. the implementation within an overall Newton-Raphson finite element solution strategy relics on the correct evaluation of the algorithmic tangent moduli which result from linearizing the slress update integrator. Likewise. these algorithmic tangent moduii are employed within the numerical analysis of the spatial localization tensor. Examples demonstrating the performance of the proposed strategy at the local constitutive level as well as al the global :;tructural level are given.

1. Introduction The underlying motivation of this contribution is to describe structural softening behavior of metals, e.g. necking or localizat;an, by the evolution of nonuniform deformations of ductile single crystals. Within this work single crystals are modeled following the traditional lines of the finite deformation constitutive framework for crystal slip plasticity, (see e.g. [1-4] among others). Thereby, we consider a classical power law type rate-dependent viscoplastic relation based on the Schmid law together with at. associated Taylor-type hardening law. h is widely agreed that material rate-dependence regularizes (to a certain extent) the impending local loss of ellipticity which results in strong mesh dependence of localization computations. The single crystal model is embedded within the structure ot finite strain multipi,~ative elasto-viseoplasticity. The multiplicative decomposition of the deformation gradient into an elastic and a plastic part is usually credited to Lee [5]. As a consequence, the cons'!tutive relations are derived with reference to the plastic intermediate conliguration, for an account of these ideas see e.g. [6] and references therein. For the model considered hardening evolution and the plastic flow rule follow from a dissipation potential and :.he hypothesis of general associativity. Emphasis is directed towards the implicit numerical implementation of the constitutive model within the framework of finite element analysis of initial inelastic boundary value problems. This is in contrast to the explicit numerical treatment, e.g. in the frame of a rate tangent method pursued by Peirce et al. [7]. It is well known, that the set of constitutive equations to be integrated is considerably stiff and "Corresponding author. 0045-7825/96/$15.00 ~) 1996Elsevier ScienceS.A. All rights reserved SSDI 0045-7825(95)00913-2

P. Stebtmann.E. Ste#t/Compat. MethodsAppl. Mech. Engrg..129(1996)235-254

236

poses severe difficulties for the numerical integration algorithm. For planar double slip, Rashid and Nemat-Nasser [SJ tackled this problem by employing a weighted residual method for the constitutive integration algorithm within an explicit finite element environment. Mohan et al. [9] pursued an explicit integration of the plastic flow in combination with an implicit approximation of the slip system shear rate within a forward gradient type treatment of the global equilibrium equations. By the way of contrast, we focus on the fully implicit incremental algorithmic counterpart of the inelastic constitutive model and the exact linearization of the resulting integration algorithm. Thereby, the constitutive algorithm relies on a family of integration methods which approximate the tensor exponent of the non-symmetric algorithmic flow direction. Satisfaction of the plastic incompressibility constraint is highlighted as a major source for the accuracy of the integration algorithm. The resulting scheme is incorporated into an implicit Newton-Raphson solution strategy on the global structural level. Localization of plastic flow in one or several shear bands is often observed in ductile single crystals, e.g. (fcc) AI-Cu alloys, when a homogeneous state of smoothly varying deformations has' reached a certain critical strain or stress level. For a rate-independent solid the condition for the onset of localization, i.e. a bifulcation into a localized band mode, is obtained from a material instability analysis (see e.g. [10-12]). This bifurcation coincides with the local loss of ellipticity of the governing quasi static equations. The acoustic, or more precisely the localization tensor inherits all informations concerning the onset of localization and the bifurcation mode. It is often argued that for rate-dependent viscoplastic solids bifurcation into a localized band mode is effectively excluded, thus the governing equations remain elliptic (see e.g. [13]). Nevertheless, in the limit of low rate sensitivity the tendency towards localization is practically identical to that observed for the corresponding rate-independent model. Another objective of this work is therefore to examine the localization behavior of the a~gorithmic single crystal model sketched above. To this end, the elliptieity condition is investigated for the spatial algorithmic tangent moduli which result from the linearization of the stress integration algorithm. These moduli serve as the kernel for the spatial localization tensor. This is in contrast to the traditional localization analysis based on the continuum description of single crystals pursued, e.g. by Asaro and Rice [2], Asaro [4] and others. As a model problem we examine the performance of the proposed integration algorithm and the resulting localization behavior for the example of simple shear. Finally, a finite element analysis of the inelastic initial boundary value problem of a single crystal in planar double slip extension involving strain localization is presented and discussed in detail.

2. Continuum mechanics framework

To set the stage we briefly review the kinematical background of finite multiplicative elastoplasticity. Assume the non-linear deformation map x = ~(X, t) : B0 x 2 ~ ~3 mapping particles labeled with X in the reference configuration to their actual position x in the deformed configuration B within the time interval 2. Then, F = Vx~ with J = detF > 0 denotes the deformation gradient with local multiplicative decomposition F = F~ - Fp into an elastic Fe and a plastic part Fp (see Lee [5]). The spatial velocity gradient is denoted by i =/~" • F - i and decomposes additively

! =]:~.F~. l +F~.F~ ,Fp I .F~, I = !~ +F~ .L o.F~ 1.

(1)

The elastoplastic decomposition introduces the notion of the plastic intermediate configuration/3o assumed as macro stress free. Then, we can define a set of strain measures associated with the multiplicative decomposition. Among them the most important one for the subsequent derivations is the elastic right Cauchy-Green tensor Cc defined in the intermediate configuration Br

co = F'~. F~.

(2)

Geometrically, this fundaalental covariant strain measure is coneected via elastic or plastic pull-back and push-forward with the spatial metric g and the right Cauchy-Green tensor, respectively

C~=Ft~.g. F c = F ~ ' . C . F ~ I.

(3)

I? Stei,mmmz, E. Stein/Comput. Ah.d~ods Appl. Mech. Enqrg. 129 (1096)235-254

237

Guided by this consideration, the reduced form of the Clausius-Duhem inequality is expressed either in terms of the Kirchhoff stress r in B and rate of deformation tensor I -'m = ~£v(g), the Second PiolaKirchhoff stress I; 0 in Bp and plastic Lie derivative of the elastic right Cauchy-Green tensor £g(Ce) or the Second Piola-Kirchhoff stress 2"_~in B~ and rate of the right Cauchy-Green tensor {7 1 ~Zp • £,;(c~)

1

~r" £v(g) - + =

l v,

- ~' = ~ . .

• C - ~, > 0.

(4)

Due to the restrictions of the principle of material frame indifference the macro part of the Helmholtz free energy ~-m~r,, is formulated most generally in terms of the elastic right Cauchy-Green tensor C¢. For an account on this argument refer, e.g. to Miehe and Stein [6]. An additional scalar internal variable ~, is responsible for isotropic hardening and is the argument of the micro part of the Helmholtzfree energy t]tmicm lit = t]trnacr°(cc) + lltmicr°(y).

(5)

Taking into account the definition of the plastic Lie derivative of the elastic right Cauchy-Green tensor £g(C~) = F ~ ' . C . F ; ' = C~ + 21C~. tpl ~:m

(6)

the elastic constitutive law, the isotropic hardening law and the remaining dissipation inequality follow as Og,........ Og,"i .... Xp=2 ~ , g = : ~ and D = [ C ~ . 2 ' p ] : L p - g - j , ) 0 . (7) Next, the framework of an associated structure for crystal plasticity is briefly outlined. To this end, the principle of maximum dissipation in its penalty formulation with ~/0 a reference shear strain rate, i.e. the penalty parameter, and $"(.) monotonic increasing penalty functions is written as 5'o~ ~ (r",g) --+ stat

- D ( ~ ' e • -~p, g) + ~1 £

(8)

6r

with the Sehmid resolved shear stress r" on slip system a which is characterized by the slip plane normal N " and slip direction M " in the plastic intermediate configuration Bp, i.e. r" = [C~. ~;p]' Z"

with

Z" = M '~ ® N " .

(9)

Please note that plastic flow is assumed as plastic volume preserving since we have simple shear with M " _L N " for each crystallographic slip system a. Typically, for fcc crystals we have a = 12 slip systems which may be activated in a three-dimensional deformation. Then the principle of maximum dissipation renders the associated flow rule for the plastic velocity gradient on Bp Lp = ~ ~

y0&,q J Z

t~

,_.,

(10)

~t

Additionally, the associated evolution law for the hardening variable ,/is obtained as 1

5' = - ~ ~

(11)

~,,,0,,'~~.

In the sequel the penalty functions q~' are chosen in the sense of a Norton creep law as ,¢

4,,, = ~-2-~ 2g g

1

+

(12)

to render the evolution equations 5'" = 5'0 sgn(r')

g

and

5' = 7 7 ] ~ ,or

g

(13)

P Steinmann, E. Stein lC~mput. Methods Appl. Mech. Engrg. 129 (1996)235-25.1

238

Finally, as K increases the following approximations are commonly made for the hardening variable evolution law K 5'"~" ~/=-Z~-IE g ~'E l~/'l

EXAMPLE 1. For the micro part

K for--~1~¢+1

i~'"l and-----,l.g

(14)

of the Helmholtzfree energy lltmicr°, we may choose, e.g. for AI-Cu

alloys ~micn,(y) = z07 +

h0

In cosh ~ \rx, - t o / /

(15)

with ~1 and ~',o the initial shear yield stress and the saturation strength and 17o the initial hardening rate. Then, the hardening law, i.e. the current shear flow stress, and its derivative a~g = h follow as g(T)=~)+[r

_r0]tanh (

h07 '),

\r~o

r0/

-

h(T)=hocosh 2( hot ~. \r~

-

(16)

'r0/

REMARK 1. Often the hardening law is stated more generally in rate form for each slip system a, e.g. Asaro [14] or Hill [15], with h,~0 the hardening moduli as 8'~("/' Y~) = Z h~'~(r)l'J'el with ~, = E I~'"1.

(t7)

Following, e.g. Hutchinson [16] the hardening moduli may be specified to a good approximation with g the latent hardening ratio as h,~(y) = h(-/)[e + [1 - e]8,,,l

with 1 <~e ~< 1.4

(18)

This type of hardening law describes latent hardening for e = 1 as well as overshoot phenomena for e > 1. A specific form of purely self hardening or Koiter hardening may be retrofitted by setting g = 0. For isotropic or Taylor hardening with g = 1 and h,,~(7) = h(y) we can integrate the hardening law in closed form to obtain g = g" -- g"(y).

3. An implicit constitutive algorithm

In this section we focus on the incremental algorithmic counterpart of the inelastic constitutive model. The set of equations to be integrated over a finite time step consists of the flow rule Lp - E ~M'~ ®N" = Z ~t'"Z'~ with

Lp = Fp,

F~. 1

(19)

together with the viscoplastic evolution law for the shear rate on slip system o~ 5'" = % sgn(r")

with r" = [C~. 27j,] • Z".

(20)

The hardening law is stated in closed form which corresponds to isotropic Taylor hardening g(y)=ro+[1.x_ro]tanh(

ho~, \r~o

-

h(y)=hocosh2( ho.____7_7 ~.

(21)

\ r x - r0/

"~,/

Within an incremental time stepping procedure the solution for the plastic deformation gradient is advanced in time via "+lFp = exp(A). "Fp = Pp. "F v with A = ~ AT" Z".

(22)

or

Here, we introduced the algorithmic flow direction A and the incremental algorithmic plastic deformation

P. Steinmarm. E. Stein/Cmnput. Metltods Appl, Mech, Engrg. 129 (1996)235-254

239

gradient/~p. The tensor exponential of the algorithmic flow direction defines a map from the set of deviatoric tensors onlo the group of isochoric tensors exp(.) • sl(3) ---, SL(3) and has the important property of conserving the plastic volume Jp : det Fp = 1. The tensor exponential has recently been extensively applied within isotropic multiplicative elastoplasticity based on the spectral representation of symmetric tensors (e.g. [6. 17. 18] among many others). However, crystal plasticity is anisotropic and the algorithmic flow direction is non-symmetric. This leads to severe difficulties for the evaluation of the tensor exponent and especially its derivative with respect to its argument. In order to bypass these difficulties in the following, the tensor exponent is approximated by the one parameter family of integrators which result from applying the generalized midpoint rule to the generating equation Lp = A exp(A) ,~ [I - 0A l t. [I + [l - 0]A 1 - / v v .

(23)

R E M A R K 2. This set of integrators contains a number of well-known methods. Firstly, by setting 0 = 1 the flow rule i,.~integrated with the fully implicit lirst-order accurate Euler backward method to give "+lFp = "Fp ~.-A . ::'"~'p --: [I - A] t. "Fp.

(24)

Secondly, 0 = renders the (1. 1) Pad6 approximation, see Weber i19]. to the tensor exponent which corresponds to the second-order accurate midpoint rule 'r+IFp= "Fp+-~A.[n+IFI~+

p]=

1-

A

- !+

A

. "Fp.

(25)

Finally, by setting 0 = 0 the conditionally stable Euler forward method is retrieved "+lFp = "Fp + A . "Fp = !1 + A]" "Fp.

(26)

This approximation has fiequently been used, e.g. in Mohan et al. [91. R E M A R K 3. Questions of A-s~ability or the stronger requirement of B-stability of integration algorithms have up to now been investigaf;ed within the context of the geometrically linear theory of elastoplastieity. Results obtained under this asriumption point to the use of 0/> ~ in order to achieve non-linear stability. MOTIVATION 1. As a moti'vation for the integration algorithms outlined above consider the case of shearing along a single slip system resulting in simple shear with Lp = y M ® N

,.~

Fp = I + y M ® N ,

F~ i =_I - r M ® N.

(27)

The evolution of the plastic part of the deformation gradient Fp = ~ M ® N is traced by straightforward integration to yield n, I 1

,,+l~ ~ --: "F v .~. j,l

~,,dt ~4 ® N = "Fp + A 3, M ~ N.

(28)

Remarkably, for this simple case the tensor expooent of the algorithmic flow direction A = Ay M ® N and its approximations render the identical result exp(A) =

,!

=[I-OA]

I.[I+[t-OIA]=I+ATM®N.

(29)

IL~{)

Clearly, this coincidence relies on M _L N and does not hold as soon as more than one slip system is taken in!o account. I N T E R L U D E I. Since plastic flow is assumed to conserve the plastic volume, i.e. Lp : I = 0, the integration algorithm has to preserve plastic incompressibility, which is exactly satisfied for detFp = ]p = 1. After some straightforward manipulations taking into account Lt = A : 1 = 0, L2 = ½A : A and

L3 = det A = .~A2 • A we arrive at ]p = 1 - [1 - O]2L:, + [1 - 0]3L3 _ 1. 1 - 02L~ - 03L3

(30)

P. Steinmann, E. Stein/Comput. Methods AppL Mech. Engrg. 129 (1996)235-254

240

For planar double slip A defines a rank two tensor with det A = 0. Therefore, setting 0 = ½ renders the only member of the above family of integration algorithms which exactly conserves the plastic volume. In a general case of multiple slip with det A # 0 a plastic volume preserving algorithm is designed by solving Eq. (30) for 0 1 1Lz 4- , / [ ! L z ] 2 1 O=2-gL--~.~ V[3 L3J-1"~

(31)

The solution of this equation with ½- ~

~< 0 ~< ½+ ~

renders an optimal adaptation of the integration

algorithm to the constraint Jp = 1. With these preliminaries at hand the elastic right Cauchy-Green deformation tensor Ce follows straightforward in terms of the incremental plastic deformation gradient/~'p and the trial elastic right Cauchy-Green tensor "C = "Fp t . "+It . "Fp I as "+lC~ • = "+tFl~' - "+lC- "+IF~1 = p~' • "co

.~'~'.

(32)

Additionally, the remaining evolution equations for the shear rates on the slip systems a are integrated fully implicitly to obtain ,,+ty. = ,,y,,+Ay.

l

n + l Tot K

with Ay" = Ay0sgn("+~¢'') ~

.

(33)

Under the assumption of an elastically stiff response of metals the Second Piola-Kirchhoff stress in Bp is given by the hyperelastic isotropic St. Venant relation

n*l.Yp = [~A["+lC~ : l - 3] - l.t] l + lzn+iCe = "+'Stl + l.t "÷tCe.

(34)

with abbreviation ,,÷l~t = ½A["÷tCe : I - 3] - g. Then. the algorithmic Schmid resolved shear stress on slip system tr is obtained as ,,+1¢, = [ ' " C ~ . "+t.Yp] : Z" = [ "+l at "+l C~ + tt"+1C2¢] : Z ".

(35)

Finally, the sum of the accumulated slip strains on all slip systems, i.e. the hardening variable, is computed dS "+lT-'ny+A"

7

with A y - - - ~ I A ~ , %

(36)

tt

Then, the system of equations R" = 0 has to be solved ~iteratively for the unknown AT" with n+l ,/.a '¢

R" = -A'y" + AT0 sgn( ,+1 ~.,,) ~

.

(37)

First, however, it is appropriate to discuss the influence of the rate sensitivity parameter ~c on the difficulties encountered in solving the equations R" = 0. Phenomenoiogically, besides a wide range of rate-dependent material behavior the rate-independent case is approached in the limit as K ~ ~ . To model rate-independence as the limiting case of a rate-dependent behavior circumvents the ambiguity encountered in truly rate independent formulations with regard to which slip systems are active and which are not. Unfortunately, as K is increased the set of constitutive equations to be integrated becomes considerably stiff. For example, a very nearly rate-independent response is recovered by setting K = 200. As a consequence, the slightest overshoot in I"+~¢'l/'+lg will prevent the residuum R" from being evaluated numerically, i.e. the computer will indicate a floating point overflow. Mc_eover, as numerical experience demonstrates, the domain of attraction for the local Newton method is drastically reduced. This observations pose severe restrictions on the applied time step size. To overcome these calamities and to assure convergence of the Newton method even for large time steps we reformulate the residual R -q a s

P.. Stebmmnn, E. Stein/Comput. Methods Appl. Merit Engrg. 129 (1996) 235-254 R" = -'V/ Ay0 sgn( a~,,,,,+1r") + I""z"l ,,+----5~-i

241 (38)

,

This system is effectively solved with the Newton method in combination with a simple linesearch that prevents the violation of sgn(Ay") = sgn( ,,+1r"). Nevertheless, as soon as A'y" ~ 0 this system is again not well conditioned since then the diagonal entry in its Jacobi matrix tends to infinity, i.e. o%:R" .--. oo. This case occurs for example in problems involving localization when large parts of the solution domain unload elastically or when the slip planes rotate very drastically. Therefore, for slip planes a with AT" ---, 0 the algorithm switches to the solution of the original equation R" = 0. The integration algorithm is summarized in Appendix A. For very constraint situations this strategy fails to converge and we employ a solution technique summarized in Appendix B which gradually increases the rate sensitivity parameter ~r from an initially low value to its final value in a number of steps. This strateg3 corresponds to the stepwise satisfaction of a constraint within an augmented Lagrange optimization raethod and converges under all conditions. Nevertheless, the total number of iterations, i.e. the sum of ti~e iterations in each step, increases drastically. To finalize this section, the derivatives Oar,, ,,~t r" and 0~: ,,+lg needed for the computation of the Jacobi rna~rix J-,t~ = -cga: R" follow as

O,+lt,, O,+t ~-,, O,,+lCe OAT.° - OCt. : OAT~

and

O,,*lg ~=h(y)sgn(AT

~)

(39)

where the remaining derivatives are given with P" = [Z"] ~y"~ by

OC~

- zX["+'C~" " P"]l + "+t s t e " +2t~[""C~ "P"]~'"' -

0"+1Ce 0A3,t~ - 2

O[~vt 0A3,t3

(40)

I 'I S'¢nl

OFp' [ /~pt. "C,,. 0A3,t~j

"

(41)

[l + [ 1 - o]a]-' . ZtJ . [Ol + ll - OlPi,' ] , •

(42)

4. Linearization of the constitutive algorithm

In this section we focus lJ the consistent linearizati-m of the constitutive algorithm which constitutes a necessity for optimal convergence within the global Newton-Raphson solution strategy. As a point of departure consider the weak form of the balance of linear momentum weighted with 8x-,~ I~ = VxSx expressed in the plastic intermediate configuration with c'°' =

2rp-IF'

dV =

E"rc.

" ol dV

(43)

Here, ~p eFel . ~r • eFet denotes the pull-back of the current Kirchhoff stress to the fixed plastic intermediate configuration "Bp at time "t. Next, note that the driving force for the evolution of ~p, and equivalently for X r, is the elastic trial deformatiot~ gradient "F~. respectively the elastic trial right Cauchy-Green tensor "Co. Then. consider linearizations defined by A(.) = d, (-(x + eAx))],=0 =

A"F~ . la- eF~ .

. AI~ = -!~ . la,

A"C~

2"F~t..:(rwm. "Fc

(44)

Finally, the linearization of the weak form of the balance of linear momentum follows in symbolic notation as AGint= f ~ I ~ : [ " F e ® " F c ] ' ~ : [20~p "Fc®"Fe]:l;i+~"[Ita.t~]dV. Therefore, the algorithmic material operator in the spatial configuration is given by

(45)

P. Steinmann. E. Stein/Comput. Methods AppL Mech. Eagrg. 129 (1096)235-254

242

AG int - f B l ~ : g~ "1~ + o" "[Ira. I ~l dv

(46)

with the spatial algorithmic material operator 2o,~r J£~ = [~Fe ® "Fe] " ~--~,-~ '[ "F~ ® "Fe].

(47)

In other words, to obtain the spatial algorithmic material operator one has to linearize ~p (the puff-back of the current Kirchhoff stress to the fixed plastic intermediate configuration "Bp at time "t) with respect to the elastic trial right Cauchy-Green tensor "Ce and subsequently perform a push-Jbrward with the elastic trial deformation gradient ~Fe to the carrent spatial configuration. In the sequel it is useful to note the relation between the spatial material operators £'f and E~ I , ' El" !~ = la" g~ : ia + o'" lI t . 141

"~

14 dr.

AGin' = fB la"

(48)

With the plastic right Cauchy-Green tensor tTpl -- "F; l . be. "F;" = ~'~t. Fp' referred to the fixed plastic intermediate configuration -"B v at time "t we have for the specific model under consideration the following implicit dependence of ..vp on eCe "~'2"p( "C~., ay°("C~)) = "*'~51~7p'+/zCr' • "C~. ~,~t.

(49)

Therefore, the algorithmic material operator in the spatial configuration follows straightforward as [

J£~ = )~be ~ b~ + 2#Zb,. + 2 Z t,

"F~.

0Ayt~

.

]

0ATtj

,],vo

Here, we made use of the chain rule to obtain ,3a,/~

= 0-9~-~ I..c~ Cp ~ + ~l O~.----rg +

(50)

..c~

"

"

,' , c

together with

o[c~ : tl

_ t

oct' -

[ ot,~'

"1 s v m

]

"

The derivatives rg~Ayzt follow easily from the push-forward of/~" = 0 with R" = R"( eCo Ay0("C¢)). i.e.

j

0ATt~ -~l OR" --t ,,t~3-7~~ = Fp . ~ . Fp

with l,,t~ -

OR" i 0'-XVt~ .c,.'

(53)

All remaining derivatives are known from the stress update algorithm. R E M A R K 4. Note that the algorithmic material operator coincides with the elastic material operator if the shear rates on the slip systems ~x are fixed, i.e. ["F~ ,~" "F~] : 20~p 0 "C~ ]a~,, " ["Fc ,':,, "F~.] = Ab~ ,~ b~ + 2~,2"t,,.

(54)

$. Localized failure at the constitutive level

A failure phenomenon which is frequently observed in laboratory experiments as well as in nature is the localization of deformations within narrow bands. Thereby, the underlying physical mechanisms vary with the particular material under consideration. Large inelastic strains accumulate inside the localization zone and form a precursor to fracture.

R Stemmamt. E. Steht/Comput. Methods AppL Mech. Engrg. 129 (19961235-254

243

On the theoretical side, the micro structure of lhe material is traditionally homogenized into a continuum, allowing point measures of stress and strain to describe the response of a given configuration under prescribed load or deformation histories. Classical continuum theory describes localization at the constitutive level as a bifurcation problem of the velocity gradient I, i.e. the spatial velocity gradient may develop a discontinuity across a singularity surface with orientation n in in the spatial configuration /3 (see [10-12]). Discontinuous bifurcation, which is synonymous witl-, the loss of local ellipfiei~y in the quasi static case. is reflected by a singularity of the spatial localization tensor q. To this end, it is assumed that the homogeneously deformed solid is subjected to quasi static increments of deformation allowing for a resulting homogeneous overall deformation state as a primary mode. Bifurcation of the spatial velocity gradient while the velocity v = #,¢(X.t) c ¢ ~ and the deformation gradient F = V x ¢ ( X , t ) maintain conlir~uity corresponds to the formation of a singular surface of order two with spatial surface normal n as a :;ccondary deformatit~n mode. With quantities (,) at the plus and minus side of the singular surface indicated by superindices (.)" and (-) and I'l denoting the jump (.)+ - (.)- of (.), this situation is characterized by ~v~=0

and

IF~=0

:rod

~l]=~'m:.~'n~0.

(55)

Thereby, the structure of the j:Jmp in the velocity gradient is determined by the Maxwell compatibility conditions where ~" denotes the loca!:,ation amplitude and the vector m describes the mode of the discontinuity. As an example, n ± m designates a shear failure mode and m [I n indicates a tensile failure mode. By requiring equilibrium for the nominal traction rate across the discontinuity and by considering a linear comparison solid in ~bc sense of Hill [20] wc conclude that localization initiates as soon as the spatial localization tensor becomes singular q.m=0

,-~

detq=(I

with

q--n.£1.n.

(56)

Thus. the direction of the vector m is determined by the underlying eigenproblem for q. Here. the fourthorder tangential material operator ~1 relates the nominal stress rate o'. i.e. the stress rate governing the spatial incremental balance equation of linear momentum, to the spatial velocity gradient o'= O - i- tr + div v or = El : 1

(57)

and includes the geometric stiffness contributions l.x "L'l "I~ = I.~ • E2 "I,~ + tr" Ilk. i~]. Within the context of quasi static boundary value problems the first occurrence of the singularity of the localization tensor det q = 0 is referred to as the loss of ellipticity, while det q = 0 indicates !he loss of hyperbolic:,ty within dynamic initial boundary value problems. Note that for the rate sensitive case a possible jump in the stress rates is connected to the jump in the strain rates via the elastic material operator thus prohibiting effectively the singularity of the localization operator (see e.g. [13, 21]). Nevertheless, the elastoplastic limit i~; gradually regained if the rate sensitivity tends to zero.

6. Numerical localization analysis Within this section we briefly discuss details pertaining to the numerical localization analysis of the algorithmic counterpart of the single crystal model employed. For the numerical analysis the localization condition for the plane strain case ~s evaluated in ter, ns of the angle 0, which is enclosed by the normal to the discontinuity surface n I = [cos O, sin O]t and the fixed background e ~ direction. Thus the coefficients qi/of the spatial localization tensor q = q~pei ®e f with respect to the fixed background coordinate system e; are evaluated as qll = [COSO]2Eit ~ [sin O]ZE33 +cos Osin O[Ei3 + E3ll ql2 = [cos 012Ei4 + [sin 012E32 + cos ,gsin O[Eiz + E34]

q21 = [cos O]2E4t + [sin tg]2E23 + cos Osin OIEzl + E43] q2." = [cos 012E44 + [sin 0t2E22 + cos Osin 0[E24 + E42].

(58)

244

P. Steinmann, E. Stein/Comput. Methods Appl. Mech. Engrg. 129 (199t51235. 254

Here, the fourth-order algorithmic spatial material operators g~' and g~ have been recast into a matrix representation with £~' ---, E r~,/]l i

|s~2,,

o¢/122 £;---....

~112

,El 121]

L~ I,,

s~ ,22 s~ ,,2 s],2tj

~2,2 e~22l/ s,-~,,- s,-,~, /

[0~ +

o

,,,2

o

cr22

0

cr~l

1

!

0 0_,2 0 " trl 2 0 o'l Ij

(59)

Then, the zeros of det q are determined by the roots of the following euartic polynomials 4 4 detq = [cos 0] 4 E ai[tan O]i = [sin O]4 E ai[tan O] [ 4+ii - 0. i=O i=0

(60)

The extrema of det q render the critical orientations for n and are given by the roots of the following quartic polynomials 4 4 0 -- [cos 0] 4 E Ai[tan O]i = [sin 0] 4 E A~[tan O] [~4÷iI. i=0 i =0 ,, Thereby, the abbreviations det E~ = EikEit - EikEi! render the constants ai and Ai as

(61)

a0 = +det Elt4a a~ = + det E'2,4+ det E ~ + det El 4. + det E~ 12 3,* .~4 - det E23 14 - d e t El4 2s a2 = +det EL,+det E34 +det E~.42 + det EL, a3 = - d e t Ei~.- det E~:~- det E ~ - det E~4

(62)

a~ = +det E ~ with A0 = al, Ai = 2a2 -4a0, A2 = 3a3 - 3al, A3 = 4 a 4 - 2a2, and A4 = -a3. This structure is invoked in our numerical analysis of localized failure of ductile single crystals at the constitutive level to determine either the zeros or the minima of det q.

7. Example: simple shear To assess the performance and the accuracy of the proposed integration algorithm and to investigate numerically the localization properties at the constitutive level the case of simple shear of a crystal with planar double slip kinematic is examined in detail. To this end, the deformation history is applied in a number of timc steps n with n = 10, n = 100 and n = 1000 to achieve a maximum amount of shear ~max where the deformation gradient has the following representation in terms of the fixed background coordinate system with base vectors ei

F=l+tel®e2

with 0 ~ < ~ < ~ .....

(63)

The material parameters are chosen to model an AI-Cu alloy with Lam6 constants A = 35104.88 N / m m 2 and tx = 23427.25 N/ram 2, initial shear yield stress To = 60.84 N / m m 2, saturation strength r~ = 109.51 N / m m 2, initial hardening rate h~ = 541.48, rate sensitivity parameter • = 200 and reference shear strain rate Y.o = 10 -~. These material constants do not incorporate a softening behavior into the constitutive description. The deformation history is applied with a shear velocity & -- 1 antil a maximum amount of shear ¢~,,x = 5 is obtained (Fig. 1). For the kinematics of the crystal we as ,ume the planar double slip model where the initial slip planes are symmetrically oriented about the ,': direction with e~ • M" = ±½, i.e. the angle enclosed by the e2 axis and the slip directions M " is ±30 °. First, we focus on the accuracy of the integration algorithms. The components ~'~t, ~'22 and r~2 of the Kirchhoff stress computed for n = 10,100,1000 timesteps with integration algorithr~s 0 = 0, ½,1 are

Sa,inmaml, E. Sa'in/('omput. Metltods Appl. Mech. Engrg. 129 (1996)235-254

245

e2 ._

r

el

Fig. 1. Simple shear: tile cifcct of ~m:,~ :: 5.

presented in Figs. 2-4, Figs. 10-12 and Figs. 6-8. It becomes obvious that for this application both, the explicit Euler forward and the implicit Euler backward method, render extremely erroneous results for large time steps. The sign and the amount of r~l and r2~ are computed completely wrong. Even the n = 1000 timestep solution does not reproduce the excellent results of the midpoint rule which in contrast produces nearly accurate results even with n - 10 timesteps. Sign and amount of ~'lI and ~'22are computed more or less exactly with n = 10 timesteps, only the detail s of the initial phase cannot be resolved. The explanation for this apparently weird behavior is the conservation property for the plastic volume of the algorithms, It was stated above that the only member of the family of integration algorithms which exactly preserves plastic incompressibility is the midpoint rule. The evolution of the determinant of the inverse plastic deformation gradient Jp I is plotted in Fig. 5 and Fig. 9 for the explicit Euler forward and the implicit Euler backward method, respectively. Clearly, exact plastic incompressibility is not even maintained for n = 1000 by both algorithms. In contrast the midpoint rule renders exactly jpi = 1 throughout all calculations. Therefore, the algorithm to be preferred for planar double slip is definitely the 0 = ½ algorithm. Next the evolution of the plastic process is briefly discussed. To this end, the resulting incremental shear rates for slip systems 1 and 2 computed with n = 1000 and 0 = ~I are depicted in Fig. 13. After an initial phase of symmetric double slip with

(.{

(l

7500.

4590.

1500.

750,

I

450.

~

"rll

150.

L_

-150[

-150

-450

-,150{

-750

-750(

0.

I.

2.

3. Amount of Shear

4.

5.

=='-~'--"T~ ' ' w ' - t 0. I. 2.

,

[ 3.

Amount of Shear

Fig. 2. S i m p l e shear; Kirchhoff .stress for 0 : 0 and t, = IO timcsteps. Fig. 3. Simple she;tr: Kirchhoff stress for ,~ -= O and it : IOO l i m e s t e p s .

4. e

5.

P Swimnamr, E. Stein/Comput, Methoth AppL Mech, Engrg, 129 (1996)235-254

246 1.15

250,

150.

1.fi9l 51),

tt = tO n = 100 n = 1000

1.03

0,97

-50,

0.91

- 150

0.85 (

,

0.

,

i

I.

,

2. 3. Amount of Shear t

~

-2~°. I ~

4.

,

o

5.

,

~.

,~

~.'

.i.

Amount of Shear

Fig. 4. S i m p l e s h e a r : K i r c h h o f f s t r e s s for 0 = 0 and n - IO1}O t i m e s t c p s . Fig. 5. S i m p l e s h e a r : d e t e r m i n a n t of i n v e r s e plastic d c f o r m a t k m g r a d i e n t for o = O.

7500•

750.

1500

450.

151}0

150. ..................... ~................................

- 150t

-1501

~

~

- -150( -750(

'-

!

i

0.

i

1.

r

,

'

_

i

2. 3. Amount of Shear t

4,

I

O.

1.

I

-- '

]

l

I

4.

2. 3,. Amount of Shear e

'1 5,

Fig. 6. S i m p l e s h e a r : K i r c h b o f f s t r e s s for 0 = I and n = I1) t i m e s t c p s . Fig. 7. S i m p l e s h e a r : Kirc,t off s t r e s s f o r 0 = I a n d .o ..-. 100 l i m c s l c p s . [.15

250.

150.

I.!1,'1

511.

37 t

n = 1000 n = 100 n=tO

-- 51).

-[50

0.91

-250.

0.85 '

0

i

I

r

I

'

I

'

2, 3, Amount o[ Shear ~.

I

4.

'

i

5.

O.

Fig. 8, S i m p l e s h e a r : K i r c h h o f f s t r e s s for 0 = i a n d ,

I.

i

i

i

2. 3, 4 Amountof Shearc

= IOOO t i m e s t c p s .

|:ig, 9. S i m p l e s h e a r : d e t e r m i n a n t of i n v e r s e plastic d c [ o r m a t i o n g r a d i e n t for 0 = I.

P. St:,imnann, E..~'win/('mnput. Method.~ AppL Mech. I?~grg 129 (1096)235-254

247

250.

250.

A -50' I

150.

~

-50,

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

-~50.

-150

~25[ 0.

I.

2.

3.

,1

o--"

5

~.

i

~. .

•

4.

Amo,mt of Shear

AlllOllllt O~'51war,

Fig, lb. Simple shear: Kirchlmlf ,~lrc~s tor 0

~ and t~

10 t i m e s l c p s .

[:ig. 11. Simple shorn: Kiachht~ff s~rcss for it

~ and n

Illll timcsWps.

2S0.

0.010

Ifi0.

0.006 rll

50.

0.002 A7 ~

0 0021 -50' 1 *0 PO

-lfiO

-0.011]

-2,50.

I

O,

l.

2. 3. Amount of Shear e

.I.

5.

O,

1.

.-

'

2.

T

3,

'

I

4.

Amount of Shear t

Fig. 12. Simple shear: KffehholT strcs,~ Lor 0 ~ 'and n 11100timcstcps. Fig. 13. Simple shear: incrcmcmal sh, ~r talcs ~11 slip s,,',~lenls I and 2.

which corresponds to the linearized format of simple shear strain, the rotation of the slip systems forces the incremental shear rates on slip system 1 to vanish gradually until they turn to zero for an amount of shear of e ~ ~. For this deformation slip system 2 is exactly aligned with the e, axis. Upon further deformation, slip-system 1 becomes active again with opposite sign of the shear rates while slip system 2 slowly tends to zero. In the tinal stage the whole deformation is carried by slip system 1 which tends to be aligned with the e~ axis. Clearly, for this situation we have ~/i ~ Ayl/At = i = 1 and Ir~l .~ ~/,V~r,¢ = 113.36 N/ram -~. Finally, the tendency towards localization at the constitutive level for the example of simple shear is examined numerically with the strategy sketched in the previous section. To this end, the n = 1000 timestep solution with 0 = ~ is considered as a numerically exact reference solution. Fig. 14 displays tile two local minimum values det qmm and the minor local maximum value of the determinant ot" the spatial localization tensor q normalized with the according global maximum value detqm x. Recall that a zero root det qmin = 0 indicates the possible onset of localization. In the first stage of the deformation op to e -,~ ~ the minimum value of detqmin is even negative, thus indicating a whole fan of possible discontinuit-es. In the following deformation regime localization is prohibited with detqn, i. > 0 until the

P. Steinmann, F. Stein/Comput. Methods Appl. Mech. Engrg. 129 (1996)235-254

248

240, 220.

O,+l 200.1~

I

's°l'WN{ ~

0.2 d,.t~

jl

-0.2

16°t

I

I40'-~ I

-0.6

IOO..~.~ -I.0 O.

," 1 , ,i'-~ f , n I.

2.

r

3.

Amount of

Sheare

4.

5.

0.

i.

2.

3.

Amountof Sheart

4.

,5.

Fig, 14. Simple shear: localization condition for 0 = ~ and n = 10OO timestcps, Fig. 15, Simple shear; discontinuity surface normal for O = ~ and n --: IIHH~ timcstcps.

satisfaction of the localization condition is gradually regained as e ---, 5 for both local minima of detqmi o. The corresponding critical angles ~ of the discontinuity surface normal n are depicted in Fig. 15. For the deformation regime with det qmin < 0 the angles O for the roots of det qm~, -- 0 are plolted in addition. Clearly, these critical angles serve as upper and lower bounds for det qmi, -- 0. A zoom of this scenario is provided in Fig. 16. Starting from a discontinuity with a surface normal characterized by ,9 ~ 135 ° t In the following deformation regime a for • = 0 the critical surface normal tends to O ~ 85 ~ for • ..~ ~. critical direction of 0 = 90 ° is obtained as e ~ 5. In summary, we conclude that especially the first part of the deformation history is extremely susceptible with respect to localization. Nevertheless the spatial orientation of the slip planes and the possible discontinuity surface do not coincide in general, which is in agreement with experimental observations [22, 23].

1.0

0.6 Hill.

0,2 det~} -0.2

tt

_..

-0.6~

-1.0 t O.

0.5

0.

i

b 2.

I.

Amount of

n

L

3.

4.

Shear

Fig. }6. Si,:nplc shear: z(mm o f the critical disconlinuity orient:ilion, Fig. 17. Simple shear: localization condition for fl = ~ and n = I(Y~) timcslcps, K = 2.

5,

P. Stt'immmn. E. Stein I Compta. Methods Appl. Met'h. Engrg. 129 (1996)235-~254

249

In order to assess the influence of the rate sensitivity parameter • we recalculate the same problem with K set to • = 2. Fig. 17 presents the results for the two local minimum values detqrni n normalized with the according global maximum value det qm,,.~"Obviously. localization is excluded throughout the whole deformation history since det q,,, > 0. Therefore, we conclude that high rate-dependency rules out localization. Nevertheless. one should bear in mind that for tinite deformations the spatial material operator £1 which constitutes the kernel of the spatial localization tensor q contains as well the geometric stiffness contributions. Thus, the localization behavior is not exclusively governed by the rate form of the constitutive response.

8. Planar double slip tension localization problem A homogeneous plane strain tension problem of a single crystal was studied earlier by Peirce et al. [24] in the context of rate-independent localization computations. Localization is typically a material instability phenomenon within solid mechanics. Different failure modes participate within this problem. In the first part of the load history the boundary conditions allow the specimen to remain in an essentially homogeneous deformation pattern. In the post peak regime a localization mode develops in the form of shear band formation which is superposed by diffuse necking at high strain levels. Thereby, the shear bands are oriented with ~ 40 ° with respect to the axis of tension. The plastic deformations accumulate within these bands and lead to the final failure of the specimen. Here, we adopt the height to width ratio of the specimen used by Peirce et al. [24] who studied a similar problem with EST triangle elements within a quasi explicit rate tangent method approach and the material constants of the previous example. Contrarily, we investigate this problem with the implicit 0 = ~.. stress update pursued in this paper within the overall Newton-Raphson solution strategy. It was emphasized above that 0 = ~ is the only algorithm which results in the exact conservation of the plastic volume for planar double slip crystal kinematics. For the spatial discretization we employ displacement gradient enhanced QIE4 elements developed by Simo and Armero [251 and further improved by Simo et al. [26]. The QIE4 element expansion proves to be extremely flexible in bending dominated problems which is essential in capturing localized failure modes (see [27]). Clearly, in order to circumvent the deterioration of the spatial representation of the problem by significantly deforming finite elements an ALE solution method would be more appropriate. Nevertheless, this question currently constitutes a major research activity and is until now not completely settled for path depended materials. One quarter of the specimen with width of 20 mm and height of 60 mm is discretized into 20 × 40 finite elements and is subjected to prescribed displacements with maximum elongation of 7.5 mm at the

5000, |.uads'tep

.t,,,,t,.jl ~ .

,50 ..............

o a , ¢

tst,.p 75

2111111.

HI{10

0.

• i i i 2, 4, 6, 8. Displacement of Top Surface

Fig. 18. Planar double slip tension localization problem.

250

P St¢inm.n., E. Stein lComput, Methods Appl. Mech. Engrg. 129 (1996) 235-254

HIIIIIIIIII|I|IIlIIIIIIIIIIIII$111111 iI I l l l l l l l l l l l l l l l l l l l itlnOlllUllllllll IiIIIIIIllllllllllUhllUllllltlllnllll I IIIIIII IIIIIIIIIII Ull Illll|lilll

dmmihm,m!!hm!mmm,l | I IIIIIIIIIIIIHIIII II II II|III|I l,hlb..,,,.,..,.1.1hd,,.,,.,] iiIIIIiIiilllllllUllllllllllUllnlllll]

dI III|IIIlIIIIIBIIIUlllIIIIIIIIIIIIIII lllllnllllllllllnl II nll iii iiii ,I....................... lh. ......... B I illllllllllllllllnllllllnllllnlllll ihlllllllllllllllllllllllllllllll|llll~l IIIIIIIII IIIIllllllllllllllllllnllllll lllll}lllllllllllUllllnlllllllllllllllJ iiiiiiiiiiiiiiiii iiiiiiiii iiiiiiiiIiii lilllilnlilll|llhllllllllhllnlllll$11 I illlllllllllllllllllllll i i i Ii i i i l t l ,Il O.................... ,,,,I,,,I,,I,,,,,II lli$llllll|lllllllllllll lUlIIIIIIIII I}IIIIIIIIIllllilIIIIlulIhlIIlU|UlIII

llllllllllllfillllglllllllllllllllllllll iiinnllllnllllllllUlllllllllllllllll] |lllnlll$1nllllllllllllllllllll|lllll lllnlllllllllllUllllllllUllllllllllU l|llllllllllllnilllllllllllllllllllllll] IIIIIII+ l l l l l l l l l l l l l | l l | i n | l l l l l l l l l l IIIIIIIIIIIIIIIItllUlIIIIIII}IIIIIlUll UlUlIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIUll IIIUlIIIIIIIIIIIIIIIII I n l l l i l l l l l II+ ililllllllllnllUllll~lnllllllllllln iiiiiiitllllllllllllllllll/llllll IIII I iiiiiiiiilUlllllllllllli!llllnlldllllll/ iiiiIiiiiiiinllllllll ii nil n i i IlilillllllllllnlllllllllllllllUllllil iIIIIIIIlUlnillilnlllllUlnllllllln illilllilnllnllllllllllnUilllllliin illlllliUlllllnlllllllllllllllltllll lilliliunlniillllllllllinlUiilltUi lilllllllllllllllllllllllltllillllllilll lillllllllllllllllllllUllllllll}llllll IlnlllllllllllllllllllllllllllUllilll,

..,,.m,,,..,.d.hmh,llttdl

lllllllllllllllllllllllllllllllllllllll:

>

IIIIIlllliUllllllll IIInllnlnllllll lullnlnllUlllllllllllllllliilllUllllll llnllllllllllllllllllllllllnllllllllll IIIIIIIiiiiiiiiiiiiiii¢iiiiiiiiiiiiiiii llllllllllillllnlllllllllllllllnlllll llllllllllllUlllllllllllllllllllllllll

Illlllllllllllllllllllllllllllllllllllll ilIIlUlIIlUlilll IIIIIlUIIILUlUlUl

innllnllllli illlllllllllnUlll iIIIilllnilli~llltllllllllnlllll lilllHiiiiiiiiiiiiiIiinlllllnllllllll IIIIIfllllllUlllllllllnlllllllllllllll iiiil¢illlllillllllllllltlllllllinlllll i}llll|nllltlllllllllllllllllllllllllll ilIIIIIIIIIITIIIlIIIIIUlIIITIInlIIIlU IIIlilllllllllllllll$1nlllllllllllllll Illllllllllllllll IIIIIIIIIUl llllllll IIIllllIIIIIIIIIlUlIIIIIIIIIIIIUlIIil lilllilll¢lltllUlHIIIIIIInlllllllllJI IllllUllllilllHIIIIIlllllnllllllnlll IllllllllllnllillUllUlnnlllUllllll II | I I I l l l l l n n l n l l l l l l l l l l l l l l l l ~ l n l I l l l l l II Illlll ll n I I I 111 } liiin lilUilllillllllllUnlllUllllll llllllllllnlnlllllllll$llllllllll|lll IIIIII IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

1:0.1509 2:0,18215 3:0.2142 4:0,24,59 5:0.2776 6:0.3092

2:0 4588 3:0.7246 4:0.990.'1 5:1.2560 6:1.5218 7:1.787,5

IIllllhilllllllllllillllllllllllllllll

IlIIIlIIIIIIUlUUlIIIIIIIIIIIIIInlII llllllnlllllUlllllllllllllllllllllllll llllllnlllllUlUnllUllllllllllllllll

Fig. It). Lt}adstcp 50: deformed conliguration and distribution of hardening variable. Fig. 20. Loadstcp 75: ¢Átfft+rm~:dconliguralion and distribution of hardening variabtc.

top surface. A total process times of T = 15 s subdivided into 150 equal time increments with At = 0.1 s is considered, thus resulting in a process velocity of the top surface of v = 0.5 mm/s. The material is modeled as nuarly rate insensitive with the constants of the previous example. Note that there is no constitutive softening included in the model. Localization is triggered by assigning a reduced initial yield strength 0.%'o and saturation strength 0.9~'x, to the lower leo element. The result|rig load displacement curve is reported in Fig. i8. After reaching the peak load, which together with the elongation at this instant essentially coincides with the results reported in [24], the load displacement curve drops dov,n steeply. Clearly, in the post peak regime we obtain a slightly different load carrying behavior in comparison to the classical results in {24] due to the difference between the truly rate-independent computation in [24] and the viscoplastic single crystal approach pursued in this contribution. Moreover, the new implicit algorithm allows to tra:e the loading history into considerably more advanced regimes of elongation of the specimen. The deformed configurations and the distribution of the internal variable y are depicted for loadstep 50, 75, 100 and 150 in Figs. 19-22. Obviously, the specimen remains in an essentially homogeneous deformation mode until loadstep 50 which corresponds approximately to peak load, see Fig. 19. Then, at loadstep 75 a shear band localization mode has clearly developed as demonstrated in Fig. 20. With ongoing deformation this shear band mode and the accumulation of plastic deformations become more pronounced at loadstep 100 which is displayed in Fig. 2I. The final configuration is depicted in Fig. 22 highlighting the localized failure mode which seems to be superposed by diffuse necking. The remaining parts of the specimen behave more or less rigidly. The sequence of Figs. 19-22 neatly monitor the evolution of the distribution of the internal variable 3'

P. Steinmann, E. Stehl/(_'omp.t. M,.thods Appl. Mech. Et2grtf. I29 (1906)235-254

251

i Fig. 21. I..oadslep llll): deformed conliguration and distribution ~ff hardening variable. Fig. 22. Loadstep 15(}: deformed conliguralit~n and dislribution of hardening variable.

which drastically accumulates in a narrow band. Thus, the tendency to develop a localized failure mode is clearly reflected.

9. Conclusions An implicit strategy for the numerical simulation of deformation processes of ductile single crystals has been proposed, The integration of the very stiff constitutive equations is performed by a simple and reliable algorithm. One of the main features of the constitutive algorithm is its versatility to adapt to the constraint of plastic incompressibility. For planar double slip this condition is exactly satisfied if the tensor exponent of the algorithmic flow direction is approximated by the second-order midpoint rule. The dramatic influence on the accuracy of the stress predictions within a given deformation history has been demonstrated for the example of simple sJaear. Another important aspect for the implementation within the global Newton-Raphson finite element solution strategy is the exact or rather consistent linearization of the stress update. It has been argued that for a constitutive formulation with reference to the plastic intermediate configuration the elastic trial deformation gradient eFe is the driving force for the constitutive algorithm. As a consequence, the algorithmic material operator follows by linearizing the

252

P Steinmamt, E. Stein/Comput. Methods Appl. Mech. Engrg. 120 (1996)235-254

constitutive response with respect to the elastic trial right Cauchy-Green tensor eC~. Within the finite element computations these algorithmic tangent moduli proved to ensure a quadra~:ic rate of convergence for the equilibrium iteration. Moreover, the algorithmic material operator constitutes the kernel of the spatial localization tensor. The numerical localization analysis at the constitutive level was performed for the case of simple shear. It turns out that especially the early stages of simple shearing are extremely susceptible with respect to localization as long as both slip systems are active with equal sign of the slip rates. Remarkably, no softening behavior was incorporated into the crystal model, The intuitive result of a discontinuity surface parallel to the direction of shearing is only approached for very high amounts of shear. Finally, the example of a planar double siip localization problem in extension demonstrates the applicability of the proposed strategy within a general non-linear finite element environment. After the onset of localization the post failure regime has been traced and the deformed configurations and the plastic zone at different ioadsteps have been reported.

Appendix A. Implicit integration algorithm for crystal plasticity

I. Trial step with F = "++F, set i = 0, il+itialize A3 ,~ = 0 C=R'.R

"F;'.C-

"C, =

~F;'

2. Set iterati.l~ c o u n t e r I = ~ + I a n d line s~,arch p a r a m e t ~ r s =

3,

]

I ) e t e r m h , . F~ "j a n d u p d a t e F ~ ' ~nd C, w~th A = ~_~,, ' A ' f ' Z ~

F'7' ----[! + [I - O]A]-'. I! - OAJ F,.' = "/';'. F; z C, = F~-t- "C.. Y~" 4.

(~ompute

SCriM|I) ~tr~',~s oi| each

~,

l

= [;:,Ic,:

2

slip ~y.~tt'ln ¢~

1 .- al - ~,1

~" = I~,c, + ,C~l: z"

5. Initialization for I = I

~';'0sgn(r~)

~oto :~

'A-~ ° = ' A ~ ° - s ~ A ? ~'

goto 3

'A'r" = &AT~ =

6, Linesearch ifsgn(r°) ~ sgn(A~.~) l s = ~.~ 2

7. U p d a t e hardening variable ~ = "'~ + ~ .

I~oi

a n d current shear flow stress g(~,)

8. Compute residuum • if

A 7 ~ _< tol

and

then

[r°[ < y /~ =-A~"

+ A % s g n ( . ¢ a} ~ V

• else

9, Check convergence if

IR[ <_ to|

then update ~" and exit

10. Compute JACOBI matrix and so|re [or new iterate L~ a J~

= -0a~o/t °

A A ~ ,~ = J ~ ' J R °

'+iA~,~ = ' A T a -}- A A ' y ~

goto 2

R Stehunann, E, St('it~IComput, Metllod.~" Appl, Mech, Engrg. 129 (19901 235-254

253

Appendix B. Implicit integration algorithm for crystal plasticity

1 Trial ~tep wld~ F = "+~b",.set ~ = 0, im'.ializ~,.-X'F'= 11and ^" = h C = F' F '•!. St,t

'C,

"F~,' G "F~'

iter~ttiorLcounter ~ ~+ I

3 I)~.termit,,/~,~tand ~tp,l:~t,'F~,I and C, with A : ~.~, A:."Z" .v:'-[/~Is 't.

01At ~ [ / - o ^ l

('IHII],IIII' 5('II:¢HD

~,

St[,'>~

=

,~11 ~,,~th ..,]iil

F, ~

"F,~' F,. ~

~>q,'ill

,,

I.! yc, 1 :q - ~,] ..

•j

C.:F;' 'O,.F;'

~" : i~c, + ~c~!: z"

*'J. t.llJdat~ ha[delmLq vmt.d,le "; = "'r + ~.~, J_X-, '[ alld , 111reta shear ll,.~w stress

q [

t

t,1[7) l 1

6. Compute resi(luv., h"' : - ~3.'" 4 3,-r,,:,gn(,") ~' i" 7. Check convergence • if

[Ri < tol

and

~: =

~'~:

update r and exit • if

IRJ
and

~<~'~*~

update ~ = min(2~,gTM}and goto fi 8. Compute .l^coBl matrix and solve for t~ewiterate A'F' d,,~ = -Oa.~'.~'

~'r'" = JS~ l'l'~

'+) A3" = '~'~" + A~7 ~' goto 2

References Ill R. Hill and J.R. Rice. ('onstitulive analysis ()f elaslo-plastic crystals at arbitrary strain. J. Mech. Phys. Solids 20 (19721401-413. [2] R.J. Asaro and J,R, Rice, Strain h~alizalion in ductile singl,: crystals. J. Mech. I'h~,,s.Solids 25 (19771 309-338. [3] R. Hill and K.S. Havner. I:'crspecfivcs in the mechanics of clasmplastic crystals, J. Mech, Phys. Solids 30 (19821 5-22. [4] R.J. Asaro, Crystal plasticity. J. Appl. Mech. 5(1 (19831 921-93.,t. [51 E.H. Lee, Elastic-plastic deformations at tinitc strains. J. Appl. Mech,. ASME 36 (19691 1-6. [6] C. Miehe and E. Stein, A canonical model of muhiolicativc clasto-plasticity: Formulation and aspects of the numerical implementation. Eur. J, Mech. A/Solids 11 (1992} 25-43. [7] D. Pierce. R.J. Asaro and A. Ncedlcman. Material rate-dependence and Iocafizcd deformation in cryslalline solids. Acla Melall. 31 (19831 1951-1976. [8] M.M. Rashid and S. Ncm~:!-N'v~ser. A consfilUliVe algorithm for r::le-d,'nen,lcn~ crwtal plasticity. Comput. Methods Appl. Mech, Engrg. 94 {1992j 2(11-228. i91 R. Mohan. M. Ortiz and ('.E Shih. An analysis of cracks in ductile single crystals. J. Mech. Phys, Solids 40 (19921291-337. [1()1 R. Hill. Acceleration waves in solids. J. Mech. Phys, Solids I(J (19621 1-16. Ill] J. Mandel. Conditions de slabilitd el poslulal dc I)ruckcr, in: L Kravtchenko and P.M. Sirieys, cds.. Rheology and Soil Mechanics {Springer-Verlag, Berlin, 1%61. [12] J.R. Rice. The localization of plastic deformation, in: W.'I: Koiter. cd.. Theoretical and Applied Mechanics (North-Holland, Amsterdam. 19771. [13] A. Needleman. Material talc-dependence and mesh sensitivity in localization problems, Comput. Methods Appl, Mech, Engrg. 67 (19881 69-86. [14] R.J. Asaro. Geometrical effects in the inhomogcneous deformation of ductile single crystals. Acta Metall. 27 (19791445-453. 1151 R. Hill, Generalized constitutive relations fl,r inc,'emental deformation of metal crystals by muttislip, J. Mech. Phys. Solids 14 (1966195-102. [i61 J.W. Hutchinson. Elaslic-plasqc behavior of polycryslallinc metals and composites, Phil. Prt~. R. Soc, London A318 (19701 247-272.

254

P. Steinmann, E, Stehl/Compl,t. Methods AppL Mech. Engrg, 120 (1996)235-254

[17] J.('. Simo, Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the inlinitesimal theory, Comput, Methods Appl. Mech. Engrg. 99 (1992)61-112. [18] P.. Steinmann, C. Miehe and E. Stein. Comparison of different linite deformation inelastic damage models within multiplicative e]astoplasticity lot ductile materials, Comput. Mech. 13 (1994) 458-474, [19] G. Weber. Computational procedures for some new rate constitutive equations for elasl,~-plasticity Ph.D. l)issertati~n, MIT USA (1988). [20] R. Hill, A general theory of uniqueness and stability in elastic-plastic solids, J. Mech. Phys. Solids 6 {I 958) 236-249. [21] B. Loret and J.H. Prevost, Dynamic strain localization in clasto-(visco)plastic solids. Part I: General formulation and onedimensional examples, Comput. Methods Appl, Mech. Engrg. 83 (1990)247-273. I22] Y.W. ('hang and R.J. Asaro, Lattice formationsand localized shearing in single crystals. Arch. Mech. 32 (1981)) 369-388. t231 Y.W. (?hang and R.J. Asaro, An experimental study of shear localization in aluminium-copper single crystals. Acta Metall. 29 (1981) 241-257, [24] D. Pierce, R,J. Asaro and A. Ncedleman, An analysis of nonuniform and localized defi~rmation in ductile single cryst;tls, Acta MetalI. 30 (1982) 1087-1119. [251 J.C. Simo and E Armero. Geometrically non-linear enhanced strain mixed method and the method of incompatible m~des. Int. J. Numer. Methods Engrg. 33 (1992) 1413-1449. [26] J,C. Simo. E Armero and R.L. Taylor. Improved versions of assumed enhanced strain tri-linear elements for 31) linite deformation problems, Comput. Methods Appl. Mech. Engrg. I I0 (1993) 359-386, [271 P. Steinmann and K. Willam. Performance of enhanced finite clement tk)rmulations in localized failure computati~ms, Comput. Methods Appl. Mech. Engrg. 90 (1991) 845-867.