Slip systems and flow patterns in viscoplastic metallic sheets with dislocations

Accepted Manuscript Slip systems and flow patterns in viscoplastic metallic sheets with dislocations Sanda Cleja-Ţigoiu, Raisa Paşcan PII: DOI: Reference:

S0749-6419(14)00067-9 http://dx.doi.org/10.1016/j.ijplas.2014.03.017 INTPLA 1770

To appear in:

International Journal of Plasticity

Received Date: Revised Date:

25 July 2013 18 March 2014

Please cite this article as: Cleja-Ţigoiu, S., Paşcan, R., Slip systems and flow patterns in viscoplastic metallic sheets with dislocations, International Journal of Plasticity (2014), doi: http://dx.doi.org/10.1016/j.ijplas.2014.03.017

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Slip systems and flow patterns in viscoplastic metallic sheets with dislocations Sanda Cleja-T ¸ igoiu and Raisa Pa¸scan University of Bucharest, Romania e-mail: [email protected] and [email protected] Abstract. The work addresses the rate dependent crystal plasticity with hardening dependent on the dislocation density, when the local and non-local (diffusion-like) evolution equations for the dislocation densities are considered. A new procedure is proposed for solving initial and boundary value problems based on the incremental equilibrium equations in terms of variational results, coupled with a rate-type constitutive model. The problems concerning the deformation of a sheet composed of a single fcc-crystal, generated by different slip systems simultaneously activated, are solved numerically for an in-plane stress state. Necessary and sufficient conditions for obtaining the in-plane stress state are proved under the assumptions of multislip flow rule, and cubic elastic symmetry. The variational problem which defines the velocity field in the actual configuration at time t is solved by the finite element method (FEM) and the current state in the sheet is defined via an update algorithm when either local (differential type) or non-local ( diffusionlike) evolution equations describe the dislocation densities. Boundary value problems resulting from the simulation of compressive and tensile tests are solved by considering all eight activated slip systems together with an activation condition. Finally, a comparison between the numerical results obtained by using the non-local (diffusion-like) evolution equations and local (differential type) equations which describe the dislocation densities is carried out. Keywords: Multislip flow rule; Local and diffusion-like evolution of dislocations; Activated slip systems; Variational equality; Metallic sheets; In-plane stress state; FEM and update algorithm.

1

Introduction

The constitutive framework for classical crystal plasticity developed by Taylor (1938), Peirce et al. (1983), and Needleman and Tvergaard (1984) does not involve size dependent effects. The length-scale dependence is associated with the level of non-homogeneities observed in crystal plasticity, see Arsenlis et al. (2004). The physical interpretation of the material scale length is still an open problem, see Kuroda and Tveergard (2008). However, at a small volume level the geometrically necessary dislocations (GND) can be experimentally evidenced, while at a larger volume level the statistically stored dislocations (SSD) are commonly considered.

1

A continuum formulation incorporating the (GND) coupled with the (SSD) has been considered by Arsenlis et al. (2004) in order to describe the plasticity of the crystal at an intermediate length scale, see also Gurtin (2002), Svendsen (2002), Cleja-T ¸ igoiu (2013). Generalizations of the classical crystal plasticity models to account for the scale dependency have been constructed involving higher-order stresses which are work conjugate with the slip gradients (see Gurtin (2002, 2003), Cermelli and Gurtin (2001), Cleja-T ¸ igoiu (2013)), or by providing a work-hardening law dependent on the plastic strain gradient (see Acharya and Bassani (2000), Ohashi (2005)). In the models developed so far, for instance in Arsenlis and Parks (1999), Cermelli and Gurtin (2001), the Burgers vector has been defined by the GND-tensor Fp curlFp in the lattice space. The GND-tensor field has been decomposed in terms of the edge and screw GND densities on the slip system, (α) (α) namely ρG(e) , ρG(s) . In Kuroda and Tvergaard (2008) and Kuroda (2011) (see also Mayeur (α)

(α)

et al. (2011)), non-local evolution laws are used for ρG(e) , ρG(s) , involving the gradient of (α)

the plastic shear rate projected on appropriate directions related to the slip system. ρG(e) (α)

and ρG(s) appear in the expression for the back stress relation through their gradients. This paper deals with the non-local viscoplastic behaviour of metallic sheets of crystalline material in the framework of a model that includes the presence of dislocations. We adopt the constitutive framework of multiplicative decomposition of the deformation gradient F into an elastic and a plastic components, based on the isoclinic configuration, within the constitutive framework of finite elasto-plasticity developed in Cleja-T ¸ igoiu and So´os (1990) F = Fe Fp .

(1)

In finite crystal plasticity, the isoclinic configurations are identified with the lattice structure following Teodosiu and Sidoroff (1976), Mandel (1972). Clayton et al. (2014) includes a third factor that should be incorporated in the multiplicative decomposition ”when accounting for lattice distortion due to the defects within.” The defects existing at the lattice level generate permanent deformations and involve changes of internal structure in crystalline materials during the deformation process. In order to numerically solve the initial and boundary value problems in crystal plasticity, the incremental virtual work principle for an update Lagrangian formulation has been globally formulated by Kuroda and Tvergaard (2008), following McMeeking and Rice (1975), in the form below Z Z ˙St · δLdV = ˆt˙ · δv dA; (2) Ω

∂Ω

see also formula (43) herein provided. Following Triantafyllidis et al. (1982), Needleman and Tvergaard (1984) the incremental form considered by Teodosiu et al. (1993), Inal et al. (2002), Kuroda (2011), (2004) is written in the form Z Z Z Z ˙ ˆt · δvdV − ( ˆt · δvdV). ∆t S˙ t · δLdV = ∆t T · δLdV − (3) Ω

∂Ω

Ω

2

∂Ω

The terms inside the bracket in (3) have been introduced for the correction of the equilibrium equation in numerical computations. Here, the rate of nominal stress is denoted by S˙ t and the rate of prescribed nominal stress on the boundary of the body in ˙ while the virtual velocity gradient δL is expressed through the deformed configuration is ˆt, the virtual velocity field δv. The rate of the nominal stress involves the time derivative of T Kirchhoff stress , where T is the Cauchy stress and ρˆ is the mass density in the actual ρˆ configuration. Various rate (differential) type constitutive relationships are adopted in the literature associated with appropriate spins. Asaro and Needleman (1985), Kuroda (2011), Kuroda and Tvergaard (2008) introduced the elastic spin We = {F˙ e (Fe )−1 }a (i.e. the skew-symmetric part of the elastic rate Le = F˙ e (Fe )−1 ), and considered the fourthorder elasticity tensor in the actual configuration to be constant. Teodosiu et al. (1993) assume that the stiffness elasticity tensor with respect to the isoclinic configuration is constant, and the elastic spin is associated with the elastic rotation tensor, namely with ˙ e (Re )T . R As far as the constitutive algorithmsfor crystal plasticity are concerned, we mention the reformulation of the constitutive framework of (crystal) plasticity with the help of an incremental variational principle has been formulated, by Ortiz and Stainier (1999), Simo and Hughes (2000), Homayonifar and Mosler (2011). In Ling et al. (2005), the performance of the numerical integration schemes in certain constitutive algorithms of rate dependent crystal plasticity has been analyzed. Another numerical algorithm performed by Peirce et al. (1983), Teodosiu et al. (1993), Evers et al.(2002), Cuiti˜ no and α Ortiz (1993), Ling et al. (2005), the slip rates, γ˙ , and the hardening parameters, ζ α , respectively, are obtained using a Newton-Raphson procedure at a given stress state. In Ling et al. (2005), the authors state that ”usually, only one or two slip systems dominate the deformation, while many other system are essentially inactive” and no activation condition is considered, thus there is no need to distinguish between active and inactive slip systems. The algorithms proposed in this paper avoid the iteration procedures based on NewtonRaphson algorithm, and take into account the activation condition, as well as make a distinction between active and inactive slip systems. Moreover, we found that eight systems could be simultaneously active in compression and tensile test, and during the deformation process there exist slip systems that become inactive as can be seen in the example given in Section 7. When the non-local evolution equations (namely partial differential equation) are employed to describe the dislocation densities, the appropriate weak formulations were an¨ alyzed, e.g. by Kuroda (2008, 2011), Arsenlis et al. (2004), Oztop et al. (2013) and the micro boundary conditions are introduced in order to accurately describe the mathematical problem, as well as the physical nature of the plastic deformation at micro-scale. In section 2, the material behaviour is described with respect to the so-called isoclinic configurations, which are identified with the lattice structure following Teodosiu and Sidoroff (1976), Mandel (1972). The differential type constitutive relation for the Kirchhoff 3

stress is prescribed in terms of the elastic rate Le . The evolution in time of the plastic distortion is described by multislip in an appropriate crystallographic system following a Schmid-type law. The evolution equations for dislocation densities are described by nonlocal (diffusion-like) laws, which account for the size effect, see also Cleja-T ¸ igoiu (2013), Cleja-T ¸ igoiu and Pa¸scan (2013). The length scale parameter is introduced into the model through the diffusion-like parameter, the hardening is strongly influenced by dislocations and the back stress is supposed to be dependent on the gradient of dislocation densities. Note that the nonlocal evolution equation incorporates the local (i.e. differential type) evolution equation given by Mecking and Kocks (1981) and Teodosiu et al. (1993), if the diffusion parameter vanishes. The rate form of the constitutive model is provided in Section 3.1. In what follows we present the following: i. We formulate a solution procedure for the initial and boundary value problem in crystal plasticity when a non-local (diffusion-like) evolution equation for the dislocation densities is considered. This is based on the incremental equilibrium equation in terms of variational equality, coupled with the rate-type constitutive model( see Section 3); ii. We provide the necessary and sufficient conditions to realize an in-plane stress state when the evolution of the plastic distortion is described by multislip flow rules in an fcc-crystal assuming that the material possesses elastic cubic symmetry (see Section 4); iii. We propose efficient and flexible algorithms which provide the numerical solution of the initial and boundary value problem by an increment-type procedure: - the velocity field is obtained based on finite element methods (FEM) for solving the variational problem which defines the velocity field in the actual configuration at time t, (see section 5) and - the current state in the sheet is defined via an update algorithm when either local (differential type) or non-local ( diffusion-like) evolution equations describe the dislocation densities (see Section 6). The problems concerning the deformation of the sheet made up from a single fcccrystal, which is generated by different slip systems that could be simultaneously activated, were numerically solved using Matlab in Section 7. For the general FEM procedure we refer to Hughes (1987) and Segerlind (1984). The activation condition of the slip systems is the key point in describing the variations in time of the plastically deformed zone, the plastic shears and the dislocation density. No geometrical imperfections have been considered and only fcc-slip systems are taken into account for the crystal. The velocity field v has been numerically determined at time tn as the solution of the variational equality (see Section 3.2). Further by applying an update procedure we provide the values of the state variables at time tn+1 , i.e. for the fields that occur in the variational equality. In what follows, two updating procedures are provided. The first update algorithm has a local characteristic using the differential type evolution equation for the dislocation densities. The second one involves the non-local (diffusion) evolution law for the dislocation densities and allows the case k zero, which corresponds to the diffusionless models. 4

Both the first and the second (k = 0) algorithms give compatible results, hence showing the flexibility and efficiency of the first one. The qualitative influence of the non-local evolution equations for the dislocation densities on the behaviour of the sheet is numerically investigated under the assumption that the boundary of the sheet is impenetrable to dislocations: (i) When a homogeneous initial distribution of dislocations is considered, compression/tensile tests are solved by considering all the possibly activated slip systems together with the activation condition. (ii) When a non-homogeneity initial distribution of dislocations is considered, then for the sake of simplicity we consider only two possibly activated slip systems. The study performed in single crystal plasticity could also be important for polycrystalline materials, which are considered to be aggregates of grains (e.g. see the procedure described by Ha and Kim (2011), Delaire et al. (2000), Klusemann et al. (2013)). The simulation of large deformability was done in Van Houtte et al. (2005), Homayonifar and Mosler (2011), Keller et al. (2012), Klusemann and Yal¸cinkaya (2013), Man et al. (2010).

2 2.1

Constitutive model with dislocation Elastic constitutive relation

In order to build the constitutive model, we start from the hypothesis that the material behaves like an elastic one with respect to the lattice (the so-called isoclinic configurations), in terms of the symmetric Piola-Kirchhoff stress tensor Π ´ Π 1³ e T e = C[Ee ] with Ee = (F ) F − I . (4) ρ˜ 2 The Piola-Kirchhoff stress tensor, Π, and the Cauchy stress, T, are related by Π = (detFe )(Fe )−1 T(Fe )−T , and consequently, in the actual configuration the elastic type constitutive equation can be expressed through the elastic strain, be , by ¢ T 1¡ = E[be ] with be = I − (Fe )−T (Fe )−1 . ρˆ 2

(5)

Herein E characterizes the elastic stiffness matrix in the actual configuration and is derived from the matrix of the elastic material, C, with the coefficients given with respect to the lattice configuration by the pushing away procedure, namely £ ¤ E[X] = Fe (C (Fe )T XFe )(Fe )T , ∀ X a symmetric tensor (6) e e e Eijkl = Fim Fjn Fks Flre Cmnsr , in the component representation. In this paper ρ˜, ρˆ and ρˆ0 denote the mass densities with respect to the lattice, current and reference configurations, respectively. 5

Remark. If according to experimental data, C is almost independent of the plastic distortion, then E inevitably depends on Fe as can be seen from (6). In the case of small elastic distortion E is identical to C. Regardless the constitutive description of the material Cauchy stress satisfies, for any motion χ of the body Ω, at any time t (χ : Ω × R → R3 ), the following equilibrium equation in Ωt = χ(Ω, t) divx T(x, t) + ρˆ(x, t)b(x, t) = 0,

x = χ(X, t) ∈ Ωt ,

(7)

u(x, t) |Γ2t = u∗ (x, t).

(8)

and the corresponding boundary condition on ∂Ωt T(x, t)n(x, t) |Γ1t = s∗ (x, t),

Here ∂Ωt ≡ Γ1t ∪ Γ2t denotes the boundary of Ωt , Γ1t ∩ Γ2t = ∅, n(x, t) is the unit outward normal at Γ1t , u(x, t) is the displacement vector. The surface loading s∗ and the displacement vector u∗ are prescribed functions with respect to the configuration at time t. At an arbitrary time τ, say at τ = t + δt, the Cauchy stress T(y, τ ) satisfies the equilibrium equation and appropriate boundary conditions. Due to the fact that the large deformation framework is adopted, the position occupied by the body and its boundary vary in time.

2.2

Incremental (rate) representation of the equilibrium equation

With a view towards the incremental (rate) representation of the equilibrium equation and the appropriate boundary value problem, and the update description of the mechanical state of the body, we estimate the deformation and the stress state of the body Ω at time τ, but with respect to the configuration of the body at time t, in an incremental process. An update Lagrangian formalism or a relative description of the motion χ, (see Truesdell and Noll, 2004) can be used. We use the relative description of the motion following Cleja-T ¸ igoiu (2000), ClejaT ¸ igoiu and Matei (2012). The relative motion from configuration Ωt ≡ χ(Ω, t) to configuration Ωτ ≡ χ(Ω, τ ) is defined by χt (·, τ ) : Ωt −→ Ωτ y = χt (x, τ ) := χ (χ−1 (x, t) , τ ) ,

∀ x ∈ Ωt ,

y ∈ Ωτ

(9)

where τ is the current moment of time. The relative gradient of the motion is given by Ft (x, τ ) = ∇x χt (x, τ ) ,

and Ft (x, τ ) = F (X, τ ) (F (X, t))−1 .

(10)

The spatial velocity can be represented in terms of the relative motion as v (y, τ ) := ∂τ χt (x, τ ) ,

and v (x, t) = ∂τ χt (x, τ ) |τ =t . 6

(11)

Thus as a direct consequence of relationships (9) and (10) ∂ ˙ Ft (x, τ ) |τ =t = F(X, t) (F(X, t))−1 ≡ L(x, t), ∂τ

where L(x, t) = ∇v(x, t)

(12)

is the spatial gradient of the velocity. Conclusion: The advantage of this description is that the position at time t is fixed and τ is the current time. Then, the following estimations hold at time τ = t + δt χt (x, τ ) = x + v (x, t) δt,

Ft (x, τ ) = I + L(x, t) δt,

(13)

as a direct consequence of formulae (11)2 and (12). We can evaluate the new position of the body as a change from configuration at time t to configuration at time τ, through the displacement ut (x, τ ) = χt (x, τ ) − x.

(14)

The nominal stress (i.e. the first Piola-Kirchhoff stress tensor) with respect to configuration at time t is given by St (x, τ ) = (detFt (x, τ )) T (y, τ ) (Ft (x, τ ))−T ,

detFt (x, τ ) =

ρˆ(x, t) . ρˆ(y, τ )

(15)

Proposition 2. 1. The nominal stress satisfies, at moment t, the following relations St (x, t) = T (x, t) d S˙ t (x, t) = ρˆ (x, t) dt

and µ

T (x, t) ρˆ (x, t)

¶

(16) − T (x, t) LT (x, t) .

2. The contact (surface) force acting on the boundary of the part P, at time τ, can be expressed in terms of St (x, τ ) as follows Z Z fc (P, τ ) = T (y, τ ) n (y, τ ) da = St (x, τ ) n (x, t) da, (17) ∂P(τ )

∂P(t)

with n(y, τ ) the unit outward normal at point y on the boundary of part P, at time τ, i.e. P(τ ) = χ(P, τ ). 3. The Piola - Kirchhoff (nominal stress) St satisfies the balance equation with respect to the configuration at time t taken as reference configuration divx St (x, τ ) + ρˆ (x, τ ) bt (x, τ ) = 0,

with bt (x, τ ) = b (χt (x, τ ) , τ ) .

(18)

The boundary value conditions on ∂Ωt , considered to be the reference configuration, are associated with the balance equation (18) St (x, τ ) n (x, t) |Γ1t = s∗t (x, τ ) ,

(χt (x, τ ) − x) |Γ2t = u∗t (x, τ ) . 7

(19)

The surface loading s∗t and the displacement vector u∗t are τ dependent prescribed functions with respect to configuration at time t. We take the time derivative of the balance equation (18) and the boundary conditions (19), at time τ = t, and obtain the rate (incremental) boundary value problem at time t divS˙ t (x, t) + ρˆ (x, t) b˙ t (x, t) = 0, S˙ t (x, t) n (x, t) |Γ1t = s˙ ∗t (x, t) , v (x, t) |Γ2t = u˙ ∗t (x, t) .

2.3

(20)

Evolution equations for plastic distortion and hardening parameters

In the proposed model, the evolution in time of the plastic distortion is described by multislips in the appropriate crystallographic system (i.e. isoclinic configuration) F˙ p (Fp )−1 =

N X

¯ α ), γ˙ α (¯sα ⊗ m

(21)

α=1

where γ˙ α is the plastic shear rate in the slip system α. The slip system initially given ¯ α is the normal to the slip plane and ¯sα is the slip in the lattice configuration, where m direction, is further deformed due to the presence of the elastic distortion Fe . In the actual configuration the slip system, (mα , sα ), is defined by the following formulae sα = Fe¯sα ,

¯ α. mα = (Fe )−T m

(22)

The orthogonality condition, sα · mα = 0, obviously holds. By differentiating the multiplicative relationships (1) with respect to time, the velocity gradient can be expressed as L = F˙ e (Fe )−1 + Fe F˙ p (Fp )−1 (Fe )−1 .

(23)

Hence the rate of elastic distortion can be expressed, from (23) together with (21) and (22), in the actual configuration, in terms of the velocity gradient L as F˙ e (Fe )−1 = L −

N X

γ˙ α (sα ⊗ mα ) ,

L = F˙ (F)−1 ,

(24)

α=1

and we note that F F˙ p (Fp )−1 (Fe )−1 = e

N X

γ˙ α (sα ⊗ mα ) .

(25)

α=1

The following internal variables are used herein: the dislocation densities ρα and hardening variables (or slip resistances) ζ α in the α−slip system, which are described 8

by appropriate evolution equations. The model will be strongly related to the presence, production and motion of dislocations inside the body. The activation condition is formulated in terms of the Schmid law |τ α − τbα | ≥ ζ α

F α ≥ 0 where F α := |τ α − τbα | − ζ α ,

⇐⇒

(26)

where τ α is the reduced shear stress in the α− slip system or the resolved shear stress τ α = τ mα · sα , τ = JT,

where J = detF ≡

ρˆ0 . ρˆ

(27)

In the above formulae, τbα is the back stress on the slip system α and is defined by ¯ α · ∇ρα ), τbα = κ2 (¯sα · ∇ρα )(m

(28)

with κ2 involving a length scale parameter. The expression of the back stress depends on the gradient of the dislocation density in the slip system α. A viscoplastic flow rule associated with the deformation process is given in the form introduced by Teodosiu and Sidoroff (1976) ¯ α ¯ α ¯n ¯ α α ¯ τ − τb ¯ γ˙ = γ˙ 0 ¯ sign(τ α − τbα )H(F α ), ∀α = 1, ..., N, (29) ζα ¯ where the back stress is not included. H denotes the Heaviside function, being defined as H(x) = 0 for any x < 0 and H(x) = 1 for any x ≥ 0. The slip on each system introduces a hardening on all slip systems. The hardening law is described in this work by an evolution law proposed - like in crystal plasticity - , in terms of plastic shear rates ζ˙ α =

N X

¯ ¯ hαβ ¯γ˙ β ¯ .

(30)

β=1

Here hαβ = hαβ (ρq ) are the components of the hardening matrix and they depend on the dislocation density, see Franciosi (1985), Quiti˜ no and Ortiz (1993). The diagonal components hαβ correspond to the self hardening, while the non-diagonal components of hαβ describe the cross hardening effects, and the latent hardening characterizes the effect of the slip systems which are not yet activated. The evolution laws of hardening are identified from single crystal experiments and the matrix (hαβ ) has been represented by Teodosiu et al. (1993) as  Ã !−1/2  Ã !1/2   X X µ αβ 1 αβ αq q q β h = a a ρ ρ − 2yc ρ , (31) K  2 q

q6=β

9

where K is a material parameter, yc denotes a characteristic length associated with the annihilation process of dislocation dipoles, ¡ αβ ¢ µ is the elastic shear modulus and b is the magnitude of the Burgers vector. a is a matrix which takes into account various types of dislocation interactions. Teodosiu et al. (1993) state that the flow rule (29) and the hardening law (30) together with (31) can describe the stages of the hardening curve of copper single crystals and satisfactory correlate multislip and latent hardening data.

2.4

Non-local evolution equation for dislocation densities

In this paper, the evolution in time of the dislocation densities is described by a nonlocal law, which account for the size effect. The non-local evolution equation for the scalar dislocation density, i.e. a diffusion-type constitutive equation, is represented as in Cleja-T ¸ igoiu (2013) in the form β2

d α 1 ρ = div(ˆ ρ0 κ2 ∇ ρα ) − ∂ρα ψ, dt ρˆ0

(32)

where β2 is a positive function and κ2 , involving a length scale parameter and ψ represents the defect energy or the so-called energy density. One could formally incorporate a local model, namely a differential type equation of the form à !−1/2 µ ¶ X 1 1 ρ˙ α = − 2yc ρα |ν α | with Lα = K ρq , (33) b Lα q6=α into the non-local model if the potential ψ in (32) is identified by considering equations (33) and (32) with κ2 = 0. Thus we have (32) with 1 = D | ν α |, β2 1 ψ = yc (ρ ) − K α 2

1 with D = , b à X

(34)

!1/2 ρβ

ρα .

β6=α

Here Lα is the mean free path of mobile dislocations on the system α. A peculiar choice κ2 = k, k a constant length scale parameter, when the initial mass density is homogeneous, leads to the non-local evolution for scalar dislocation densities µ ¶ ∂ψ α α ρ˙ = D k∆ρ − α |ν α | , α = 1, . . . , N. (35) ∂ρ k is called diffusion parameter. 10

Remark. The local evolution equation (33) is given by Mecking and Kocks (1981), and Teodosiu et al. (1993). A diffusive evolution equation non-linear and parabolic has been considered in Bortoloni and Cermelli (2004), where the energy density ψ is chosen to be a non-convex polynomial. In Cleja T ¸ igoiu and Pa¸scan (2013), the expression of the back stress involved in the activation condition for the α− slip system is given by ¯ α × (αN )T ¯sα ) + κ2 (¯sα · ∇ρα )(m ¯ α · ∇ρα ), with two scale parameters, τbα = div(κ3 m p −1 p κ2 , κ3 and αN = (F ) curlF the so-called Noll’s dislocation density. As a particular case, namely when κ3 = 0, the expression of the back stress (28) follows.

3

Rate form of the constitutive model and variational equality

We propose here the procedure to solve the initial and boundary value problem. First we remark that the evolution in time of plastic shears, hardening parameters, slip systems and so on, are introduced by appropriate differential type equations in the subsections 2.3 and 2.4, while the elastic constitutive equation is written by (5). By differentiating the constitutive equation (5) with respect to time, we obtain " N # µ ¶ ³ ´T T³ ´T X d T −1 −1 − F˙ e (Fe ) − F˙ e (Fe ) = E [D] − E γ˙ α {sα ⊗ mα }S , (36) dt ρˆ ρˆ ρˆ α=1 1 (L + LT ) is the rate of strain. Eliminating the rate of elastic distortion in 2 (36) via (24), the rate type constitutive equation which giving rise to the Cauchy stress T is given by where D =

d dt

µ ¶ N h i X T T T T = L + L + E [D] − γ˙ α E {sα ⊗ mα }S ρˆ ρˆ ρˆ α=1 −

N X α=1

γ˙ α (sα ⊗ mα )

T T − ρˆ ρˆ

N X

(37)

γ˙ α (mα ⊗ sα ) .

α=1

Second we summarize the constitutive and evolution equations introduced in Sections 2.3, 2.4 and (37) defining the model in an incremental, or rate form.

11

3.1

Rate form of the adopted constitutive model

For a given history of the deformation gradient, t → F(t), t ∈ [t0 , t∗ ), the differential system defining the fields T, Fe , γ α , sα , mα , ζ α , ρˆ, and ρα is given in the following form d dt

µ

T ρˆ

¶

h i X T T T + L + E [D] − ν α E {sα ⊗ mα }S − ρˆ ρˆ N

=L

α=1

N X

−

N

ν α (sα ⊗ mα )

α=1

˙ e (Fe )−1 F

=L−

T TX α − ν (mα ⊗ sα ) ρˆ ρˆ α=1

N X

(38)

ν α (sα ⊗ mα )

α=1

γ˙ α = ν α , α = 1, . . . , N

s˙ α

= Lsα −

N X

³ ´ ν β sβ ⊗ mβ sα ,

α = 1, . . . , N

β=1

˙α m

=

−LT mα

+

N X

³ ´ ν β mβ ⊗ sβ mα , α = 1, . . . , N (39)

β=1

ζ˙ α

=

N X

¯ ¯ ¯ ¯ hαβ ¯ν β ¯ ,

α = 1, . . . , N

β=1

ρˆ˙ where να

=

¯ ¯

= −ˆ ρ tr L,

1 ρˆ0 γ˙ 0α ¯¯ α ( Tmα ζ ρˆ

α

·s −

¯n µ ¯ ρˆ0 α ¯ τb )¯ sign Tmα ρˆ

¶ α

·s −

τbα

H (F α ) , , α = 1, . . . , N, (40)

¯ α · ∇ρα ) τbα = κ2 (¯sα · ∇ρα )(m The non-local evolution equation for the dislocation density is given by µ ¶ ∂ψ α α α ρ˙ = D |ν | k∆ρ − α , α = 1, . . . , N. ∂ρ

(41)

Initial conditions need to be prescribed for the differential system, as well as a boundary value problem has to be defined in connection with the partial differential equation (41). For instance, we obtain the following boundary condition k

∂ρα = iα (ρq ) ∂n

12

(42)

were iα (ρq ) given scalar functions In formula (42) the derivative with respect to the normal direction for the density of dislo∂ρα cations, = ∇ρα · n, is given on the boundary of the deformed domain. ∂n Remark. In the right-hand side of the differential type equations, ν α replaces γ˙ α inside the evolution equations for the other variables to emphasize the model as standard (normal) differential systems.

3.2

Variational equality

As we have already remarked, the Cauchy stress satisfies the equilibrium equation at every moment t and, consequently, the rate of the nominal stress satisfies the incremental boundary value problem (20). The weak formulation of the rate boundary value problem given by (20) is the following Z Z Z ˙St · ∇wd V = ˙St n · wda + ρˆb˙ t · wd V. (43) Ωt ∂Ωt Ωt for any kinematically admissible field w. We express the rate of the nominal stress, S˙ t , from (16) together with (38)1 as S˙ t (x, t) = L T + ρˆE [D] − ρˆ

N X

h i ν α E {sα ⊗ mα }S −

α=1

−

N X

ν α (sα ⊗ mα )T − T

α=1

(44) N X

ν α (mα ⊗ sα )

α=1

If S˙ t is replaced in (43) by (44), we formulate the following result: Theorem 1. If the activation condition is formulated in terms of Schmid’s law, the rate type boundary value problem at time t leads to an appropriate variational equality to be satisfied by the velocity field, v, when the current state of the body, namely the Cauchy stress, T, the elastic distortion, Fe , the position of the slip systems, (mα , sα ), the dislocation densities, ρα , the mass density ρˆ, and the hardening variables, ζ α , are all known. The corresponding variational equality is given by Z Z (∇v)T · ∇wd V + ρˆE[D] · ∇wd V − Ωt

−

Ωt

N Z X

h i ν α {ˆ ρE {sα ⊗ mα }S + (sα ⊗ mα ) T + T (mα ⊗ sα )} · ∇wd V =

α=1 Ωt

Z = Γ1t

Z s˙ ∗t · wda +

Ωt

ρˆb˙ t · wd V,

0 ∀w ∈ Vad ,

with ν α given by (40).

13

(45)

0 denotes and is given by the set of kinematically admissible velocity fields at time t Here Vad 0 Vad = {v : Ωt −→ R3 |

v |Γ2t = 0},

for the homogeneous condition on Γ2t . v0

Remark. In the case of non-homogeneous condition for the velocity on Γ2t , if a function 0 , where ˆ = v − v0 ∈ Vad ∈ Vad is built, then ∀v ∈ Vad , the function v Vad = {v : Ωt −→ R3 |

v |Γ2t = u˙ ∗t (x, t)} .

We remark that ρˆE has the dimension of the stress as a consequence of the constitutive equation written in (5).

3.3

How to solve the initial and boundary value problem

The initial and boundary value problem The Cauchy stress, T, must be such that it satisfies the equilibrium equation (7), with the boundary conditions given by (8) and (42). The behaviour of the elasto-plastic body Ω be described at any time t in Ωt by the rate type constitutive equation (38)-(40) and the partial differential equation (41) describing the evolution of the dislocation density. The appropriate initial conditions are known. To solve (45), (38)-(40) and (41), the following steps are followed: I Solve an elastic problem on a certain time interval [0, t0 ), for all material points in the current domain, Ωt , for all t ∈ [0, t0 ) : Fp = I,

T be such that | τ α |≤ ζ0α , ∀ α ∈ {1, N }. ρˆ

II Assume that, at time t0 , the plastic behaviour is initiated, i.e. (∃) x0 ∈ Ωt0 | τ α (x0 , t0 ) |≥ ζ0α (x0 , t0 ), for a certain α. At time t0 , for all point in the deformed domain Ωt0 , the deformation and stress state of the material are known, namely T(t), F(t), Fe (t), γ α (t), mα (t), sα (t), ρα (t), ζ α (t), ρˆ(t). They have the meaning of initial data in the behaviour of the material following the moment t0 . We propose a novel algorithm to describe the behaviour of the elasto-plastic body, for a generic time t : • At time t, for any x ∈ Ωt , let us consider that T(t), F(t), Fe (t), γ α (t), mα (t), sα (t), ρα (t), ζ α (t), ρˆ(t) are known. Find the velocity field v(t) : Ωt → R3 as solution of the variational equality (45) together with (40) using the (FEM);

14

• If the velocity field is known, find the solution of the differential equations (38)-(40) and of the non-local evolution equation for the dislocation densities (41), by applying explicit Euler methods at every material point in the current configuration of the body for the differential equation. We also use Crank-Nicolson’s method for the non-local (diffusiontype) evolution equation for the dislocation densities. We consider ζ α (t0 ) = ζ0 and γ˙ α (t0 ) = γ˙ 0 . Note that the Newton-Raphson procedure is not used in the present numerical algorithms.

4

In-plane stress state

In the following, we consider the fcc-crystallographic slip system characterized by the elements given in Table 1, see Cleja-T ¸ igoiu (1996), Bortoloni and Cermelli (2004). Generally, 24 slip systems can be considered taking also into account the slip system which corresponds to −¯sα .

¯ α in fcc-crystallographic Table 1: – The slip directions ¯sα and normals to the slip planes m slip system written with respect of a given Cartesian basis 1 ¯s1 = √ 2 1 2 ¯s = √ 2 1 ¯s3 = √ 2 1 4 ¯s = √ 2 1 5 √ ¯s = 2 1 ¯s6 = √ 2

4.1

£ ¤ 1, 1, 0 £ ¤ 1, 0, 1 £ ¤ 0, 1, 1 [1, 0, 1] [0, 1, 1] [1, 1, 0]

1 ¯1 = √ m 3 1 2 ¯ =√ m 3 1 ¯3 = √ m 3 1 4 ¯ =√ m 3 1 5 √ ¯ = m 3 1 ¯6 = √ m 3

(1, 1, 1) (1, 1, 1) (1, 1, 1) ¡ ¢ 1, 1, 1 ¡ ¢ 1, 1, 1 ¡ ¢ 1, 1, 1

¤ 1 £ ¯s7 = √ 1, 1, 0 2 ¤ 1 £ 8 ¯s = √ 1, 0, 1 2 ¤ 1 £ ¯s9 = √ 0, 1, 1 2 1 10 ¯s = √ [1, 0, 1] 2 1 11 √ ¯s = [0, 1, 1] 2 1 ¯s12 = √ [1, 1, 0] 2

¢ 1 ¡ ¯ 7 = √ 1, 1, 1 m 3 ¢ 1 ¡ 8 ¯ = √ 1, 1, 1 m 3 ¢ 1 ¡ ¯ 9 = √ 1, 1, 1 m 3 ¢ 1 ¡ 10 ¯ = √ 1, 1, 1 m 3 ¢ 1 ¡ 11 √ ¯ = m 1, 1, 1 3 ¢ 1 ¡ ¯ 12 = √ 1, 1, 1 m 3

In-plane stress and generalized plane deformation states

We derive the necessary conditions for a plane stress state   T11 T12 0 T =  T12 T22 0  , 0 0 0

15

(46)

together with the deformation gradient, F = I + ∇u, describing a so-called generalized plane deformation state, namely   F11 F12 0 F =  F21 F22 0  . (47) 0 0 F33 The back stress, τbα , is assumed to be zero. Proposition 3. Let us assume that an in plane stress state (46), combined with a generalized plane deformation state (47), is realized at every moment. Then the plastic and elastic distortions have the same form as the deformation gradient (47). Moreover, the slip systems become active in pairs at the same time. Proof. The first Piola-Kirchhoff stress tensor, S, has namely  S11 S = (detF)TF−T , S =  S21 0

the same form as the Cauchy stress,  S12 0 S22 0  . 0 0

(48)

We look at the reduced shear stress defined by (27) together with (22) written in terms of the first Piola-Kirchhoff stress tensor ¯ α · ¯sα , τ α = FeT SFpT m

or (49)

¯ α · ¯sα , τ α = Mm

for M = FeT SFpT .

We proceed to analyze the possibility of having an activated slip system starting from the expression of the initial values of the reduced shear stress, i.e. Fp = I, Fe = F. We obtain from (49)   M11 M12 0 M =  M21 M22 0  . (50) 0 0 0 The following relationships are obtained from (49) calculated for all slip systems given in Table 1.: τ1 =

√1 6

(M11 + M12 − M21 − M22 ) ≡ −τ 7 , τ 2 =

√1 6

(−M11 − M12 ) ≡ τ 4 ,

τ3 =

√1 6

(−M21 − M22 ) ≡ τ 5 ,

τ6 =

√1 6

(M11 − M12 + M21 − M22 ) ≡ −τ 12 , (51)

τ9 =

√1 6

(M21 − M22 ) ≡ τ 11 ,

τ8 =

√1 6

(−M11 + M12 ) ≡ τ 10 .

Note that only four reduced shears are independent in (51). We remark that τ 1 + τ 7 = 0, τ 2 = τ 4 , τ 3 = τ 5 , τ 6 + τ 12 = 0, τ 8 = τ 10 , τ 9 = τ 11 ,

16

(52)

which proves that the pairs of slip systems become active at the same time. If the plastic shear rates are calculated by the following expression ¯ α ¯n ¯ ¯ α α ¯τ ¯ ν = ν0 ¯ α ¯ sign(τ α )H(|τ α | − ζ α ), (53) ζ as we can prove that ζ α have equal modulus values for the appropriate pairs (this is the case here), then we can conclude that ν 1 = −ν 7 , ν 6 = −ν 12 , ν 2 = ν 4 , ν 3 = ν 5 , ν 8 = ν 10 , ν 9 = ν 11 .

(54)

Thus the evolution equation for the plastic distortion is given by ˙ p (Fp )−1 = ν 1 (¯s1 ⊗ m ¯ 1 − ¯s7 ⊗ m ¯ 7 ) + ν 2 (¯s2 ⊗ m ¯ 2 + ¯s4 ⊗ m ¯ 4) F ¯ 3 + ¯s5 ⊗ m ¯ 5 ) + ν 6 (¯s6 ⊗ m ¯ 6 − ¯s12 ⊗ m ¯ 12 ) +ν 3 (¯s3 ⊗ m

(55)

¯ 8 + ¯s10 ⊗ m ¯ 10 ) + ν 9 (¯s9 ⊗ m ¯ 9 + ¯s11 ⊗ m ¯ 11 ). +ν 8 (¯s8 ⊗ m Consequently, the expression of the rate of plastic distortion in a matrix form is given by:  1  ν − ν2 + ν6 − ν8 ν1 − ν2 − ν6 + ν8 0 ˙ p (Fp )−1 = √2  −ν 1 − ν 3 + ν 6 + ν 9 −ν 1 − ν 3 − ν 6 − ν 9 . 0 F (56) 6 2 3 8 9 0 0 ν +ν +ν +ν The system of differential equation for the components Fijp is derived from equation (56) and contains separately the following equations p F˙3k =

√2 (ν 2 6

p + ν 3 + ν 8 + ν 9 )F3k ,

k = 1, 2, 3,

p F˙13 =

√2 (ν 1 6

p − ν 2 + ν 6 − ν 8 )F13 +

√2 (ν 1 6

p F˙23 =

√2 (−ν 1 6

p − ν 3 + ν 6 + ν 8 )F13 −

p − ν 2 − ν 6 + ν 8 )F23 ,

√2 (ν 1 6

(57)

p + ν 2 + ν 6 + ν 8 )F23 .

p p p p Since at the initial moment t0 , Fp (t0 ) = I, F13 , F23 , F31 , and F32 vanish, while p F33 =e

√2 (γ 2 +γ 3 +γ 8 +γ 9 ) 6

on a small time-interval. Thus, we obtain  p  p F11 F12 0 p p F22 0 . Fp =  F21 p 0 0 F33

(58)

To complete the proof, let us consider the maximal time interval on which Fp takes the form (58), say [t0 , t1 ). Using the multiplicative decomposition of the deformation gradient, it follows that the elastic distortion Fe can be expressed in a similar form. Under the assumption that the in-plane stress state is developed µ ¶ during the deformation d T3j considered process, the components T3j , j = 1, 2, 3 vanish or = 0 for j = 1, 2, 3. dt ρˆ

17

4.2

Realization of an in-plane stress state

We assume that C characterizes the appropriate cubic symmetry of fcc-crystal, with respect to the isoclinic (lattice) configuration. Then, matrix C is characterized by three elastic material constants only, C1111 = C2222 = C3333 ≡ c¯11 , C1122 = C1133 = C2233 ≡ c¯12 , C1212 = C1313 = C2323 ≡ c¯44 , see Ting (1996). Setting the following correspondence between the pairs of indices {(i, j)} for i, j ∈ {1, 2, 3} and the set {k}, k ∈ {1, . . . , 6}, defined as follows: (1, 1) −→ (1), (2, 2) −→ (2), (3, 3) −→ (3), (1, 2) −→ (4), (1, 3) −→ (5), (2, 3) −→ (6), the standard notation will be used for the matrices C, and E, respectively. Proposition 4. In the hypotheses that F and Fe characterize an mation state, namely    e e F11 F12 0 F11 F12 e e F22 F =  F21 F22 0  , Fe =  F21 0 0 F33 0 0

appropriate generalized defor 0 0 , e F33

(59)

the following statements hold: (i) The fourth-order elastic tensor in the following symmetric matrix  E 11 E 12  E 22   (E)i,j=1,6 =    

actual configuration (6) is characterized by the E 13 E 14 E 23 E 24 E 33 E 34 E 44

0 0 0 0 0 0 0 0 E 55 E 56 E 66

    ,   

(60)

Explicit expressions for E ij are given in Appendix B.

³ ´ ˙ p (Fp )−1 (ii) The components (1, 3), (3, 1), (2, 3) and (3, 2) of the rate of plastic distortion F vanish. (iii) The shear stress components are vanishing, i.e. T13 = T23 = 0, provided that they vanish at the initial moment. (iv) The time derivative of the stress component T33 is given by µ ¶ ³ ´ ¢T ¡ d T33 33 ˙ p (Fp )−1 }(Fe )−1 = 2 L33 − Fe {F + dt ρˆ ρˆ 33

(61)

˙ p (Fp )−1 )(Fe )−1 }S )kl . +E33kl Dkl − E33kl ({Fe (F Proof. Starting from (6), we obtain directly the non-vanishing elastic coefficients, which depend on the cubic elastic coefficients and the components of the elastic distortion. The issue e = F e , i=1, 2, concerning the fact that the elastic distortion has zero shear components, Fi3 3i becomes essential for finding the elastic coefficients mentioned above. As a direct consequence of the hypotheses concerning the expression of F and Fe , and of the multiplicative decomposition of the deformation gradient, it follows that Fp , the rate of plastic

18

˙ p (Fp )−1 , and its symmetric part have the same form. Moreover, the rate of plastic distortion, F distortion (21) is necessarily restricted to satisfy the following relationships: ν 1 − ν 2 + ν 4 + ν 7 = 0, ν 8 − ν 9 − ν 10 + ν 11 = 0,

ν 6 − ν 8 + ν 10 + ν 12 = 0 ν 2 + ν 3 − ν 4 − ν 5 = 0.

(62)

Consequently, the following shear components of the tensorial fields, and their symmetrical parts vanish ˙ p (Fp )−1 )kl = 0, (Fe {F ˙ p (Fp )−1 }(Fe )−1 )kl = 0, Lkl = 0, (F

(63)

for (k, l) ∈ {(1, 3), (2, 3), (3, 1), (3, 2)}. We look at the rate type constitutive equation for the stress written inµ(381 )¶togetherµwith¶(25) and projected on i3 . First we consider the equations d T13 d T23 giving rise to and see (A4). Under the hypothesis that at the initial moment dt ρˆ dt ρˆ (T13 (t0 ), T23 (t0 )) = 0, only the trivial solution of the system is obtained (T13 = 0, T23 = 0). Next, we derive the differential type equation which describes the evolution in time of T33 µ ¶ d T33 Tk3 T3k ˙ p (Fp )−1 (Fe )−1 }S )kl − = L3k + L3k + E33kl Dkl − E33kl ({Fe F dt ρˆ ρˆ ρˆ (64) ³ ´ T ³ ´ T k3 3k ˙ p (Fp )−1 }(Fe )−1 ˙ p (Fp )−1 }(Fe )−1 − Fe {F − Fe {F . ˆ ρˆ 3k ρ 3k We evaluate the terms in equation (64) using (63) and hence we obtain equation (61). Let us introduce the following notation N X

ν α sα ⊗ m α =

α=1

N X

¯ α )(Fe )−1 ≡ G. ν α Fe (¯sα ⊗ m

(65)

α=1

for the expression obtained from (25) together with (21) and (22). We can now formulate the main result that follows directly from the proposition above: Theorem 2. Let us assume that the material has cubic elastic symmetry with respect to the isoclinic configuration. In the hypotheses that F and Fe characterize an appropriate generalized deformation state (59) and D33 is given by D33 =

1 (E D11 − E 32 D22 − 2E 34 D12 E 33 31

(66)

+(E 31 − E 33 )GS11 + (E 32 − E 33 )GS22 + 2E 34 GS12 ), then an in-plane stress state is realized in the elasto-plastic material if the components of the Cauchy stress, at time t0 , T13 (t0 ), T23 (t0 ) and T33 (t0 ), vanish. Proof. As a consequence of Proposition 4, the time derivative of the stress component T33 is given by (61). Let us remark that T33 = 0 on a certain time interval, say [t0 , t1 ), if and only if ˙ p (Fp )−1 )(Fe )−1 }S )kl = 0, E33kl Dkl − E33kl ({Fe (F

19

(67)

on [t0 , t1 ). Using the component representation of the expression written in the brackets in relation (67), formula (67) together with (65) is equivalently given by E 31 D11 + E 32 D22 + E 33 D33 + 2E 34 D12 − E 31 GS11 − E 32 GS22 − E 33 GS33 − 2E 34 GS12 = 0.

(68)

Since the plastic incompressibility condition is satisfied by the multislip flow rule written either with respect to reference configuration or to the actual configuration, namely trG = 0, we can eliminate the component G33 . Finally, we write (68) in the form (66).

5

Finite element discretization

The finite element method (FEM) is applied for solving the variational problem at time t defining the velocity field in the actual configuration, v, together with a temporal discretization of the rate type constitutive equations to update the current state in the sheet. The dislocation densities should be treated as field variables which evolve in time and space. The solution of the variational problem (45) formulated at moment t is obtained using the FEM. Let us consider (tn )n=1,N ∗ a partition of the time interval [0, T ] with tn+1 = dt + tn . We assume that the current values of the fields are known at moment tn , i.e. xn , Tn , Fn , Fen , sαn , mαn , γnα , ραn , ζnα , ρˆn .

(69)

Let Ωtn be the domain occupied by the body at moment tn . We require the variational problem at tn be satisfied by the velocity field v : Z Z n £ ¤ o (∇v)T · ∇w dV + ρˆE [D] · ∇w dV − ρˆE GS + GT + TGT · ∇w dV = ΩZ Ωtn Ωtn tn Z (70) ∗ ˙ = s˙ t · w da + ρˆbt · w dV, ∀w ∈ Vad , Z

Γ1tn

Ωtn

where the notation introduced in (65) has been used. The FEs are chosen to be triangles and the shape functions are considered to be linear on each element. By applying a standard FEM, the discretized variational equality can be written in matrix form as: ˜ v = −K¯ ¯ v + f, K˜ ˜ K

=

X e

where

˜ e )T Ke A ˜ e, (A

¯ = K

X ˜ e )T Ke A ¯ e, (A

f=

e

X ˜ e )T f e . (A

(71)

e

˜ is the vector built with the components of the velocities at the FEM nodes, except those Here v that are lying on the boundary where velocities have been imposed. On the contrary, the vector ¯ contains all FEM nodal values on the boundary where velocities are given. Consequently, we v obtain: ¡ ¢ Ke = Be2 T Be1 + ρˆBe4 T EBe3 Ae , Z

fe

¡ ¢ = ρˆBe4 T EGS − ρˆBe4 T Eq + Be4 T Q Ae +

Be1

= ∆1 NeT , Be2 = ∆2 NeT , Be3 = ∆3 NeT , Be4 = ∆4 NeT ,

20

∂Ωetn ∩∂Ωtn

Ne s˙ ∗t da,

(72)

· NeT

=

N1e 0 N2e 0 N3e 0 0 N1e 0 N2e 0 N3e



E

q

E 11 E 12 E 13 E 14 =  E 21 E 22 E 23 E 24  , E 41 E 42 E 43 E 44

=

,

(73)

 S G    11 GS22 GS = GS    33 GS12



     

¸

0 0

      

,

ª 1 © E 31 GS11 + E 32 GS22 + E 33 GS33 + 2E 34 GS12    E   33 0

          

(74)

,

 

 2(G11 T11 + G12 T21 )  2(G21 T12 + G22 T22 ) Q = .   G11 T12 + G12 T22 + G21 T11 + G22 T21 Here Ae denotes the measure of the surface associated with the finite element e in the deformed configuration at moment tn , while ∆1 , ∆2 , ∆3 , ∆4 are appropriate differential operators:    ∂  ∂ ∂ T11

+ T12

∂x1 ∂x2    0  ∆1 =   T12 ∂ + T22 ∂  ∂x1 ∂x2  0



∂ ∂x1

   0  ¶ ∆3 =  1 µ ∂ ∂  −E 31 − E 34  E ∂x1 ∂x2  33 ∂ ∂x2

6

0

  ∂x1  ∂ ∂    0 T12 + T22  ∂x1 ∂x2   , ∆2 =  ∂   0   ∂x 2   T11

∂ ∂ + T12 ∂x1 ∂x2

0

0

∂ ∂x2 0 ∂ ∂x1

    ,   

(75)



 ∂   ∂  ∂x 1   ∂x2 µ ¶   0 , ∆ =  4  1 ∂ ∂  −E 32 − E 34  ∂ E 33 ∂x2 ∂x1  



0

∂ ∂x1

∂x2

0 ∂ ∂x2 ∂ ∂x1

   . (76)  

Time discretization

The update algorithm is related to the temporal discretization of the rate type constitutive equations, see Section 3.1. The first update algorithm has a local characteristic, using the differential-like, evolution equation for the dislocation density, (35) with k=0, i.e. (33).

21

6.1

Temporal discretization of the differential system

Let us introduce the initial conditions at time t0 , corresponding to the differential system (38)¯ α , γ α (t0 ) = γ0 , (40) and (33): T(t0 ) = 0, F(t0 ) = I, Fe (t0 ) = I, sα (t0 ) = ¯sα , mα (t0 ) = m α α ρ (t0 ) = ρ0 , ζ (t0 ) = ζ0 , ρˆ(t0 ) = ρˆ0 . When the velocity field v has been numerically found at time tn , the velocity gradient, L, and the rate of strain, D, follow at once. Further, by applying an update procedure to the appropriate differential-like relationships associated with the constitutive equations describing the model, i.e. an explicit Euler algorithm, we provide the values, at time tn+1 , for the fields listed in (69). We describe the procedure to update the current state variables in the sheet, namely we compute, at time tn+1 , the following values: α α xn+1 , Tn+1 , Fen+1 , sαn+1 , mαn+1 , γn+1 , ραn+1 , ζn+1 , ρˆn+1 .

(77)

We take τ = tn+1 and t = tn , in the formulae (13) and one obtains xn+1 = xn + vn dt,

Ft(n+1) = I + Ln dt.

(78)

Using (44) and (65), the numerical value of the nominal stress St (x, τ ) can be calculated at time τ = tn+1 as St (x, τ ) ≈ T(x, t) + (L(x, t)T(x, t) + ρˆ(x, t)E[D(x, t)]− −ˆ ρ(x, t)E[GS (x, t)] − 2 {G(x, t)T(x, t)}S )(τ − t), or, equivalently, ¡ ¢ St(n+1) = (I + Ln dt)Tn + ρˆn En [Dn ] − ρˆn En [GSn ] − 2{Gn Tn }S dt.

(79)

As a direct consequence of definitions (15) and (10) together with (38)2 and (65) the updated values of the Cauchy stress, mass density, deformation gradient and elastic distortion are obtained

Fn+1

1 S FT , detFt(n+1) t(n+1) t(n+1) = Ft(n+1) Fn ,

Fen+1

= Fen + (Ln − Gn )Fen dt

Tn+1 =

22

ρˆn+1 =

ρˆn , detFt(n+1) (80)

In view of (39), (38)3 and (33), the update procedure is completed by considering sαn+1 = sαn + Ln sαn dt − Gn sαn dt, α = 1, . . . , N, mαn+1 = mαn − LTn mαn dt + GTn mαn dt, α = 1, . . . , N, N ¯ ¯ X ¯ β¯ α α ζn+1 = hαβ n ¯νn ¯ dt + ζn , α = 1, . . . , N, α γn+1

ραn+1 Lαn

β=1

= νnα dt + µ γnα , ¶ 1 1 α α = ρn + − 2yc ρn |νnα | dt, b Lαn  −1/2 X =K ρqn  , α = 1, . . . , N.

(81)

q6=α

As (80)–(81) are calculated, the following approximations for ν α , G and F α can be derived from relations (40) (with n replaced by 1/m), (65) and (26): ¯ ¯ ¯ ρˆ0 Tn+1 mαn+1 · sαn+1 ¯1/m ¡ ¢ ¡ α ¢ α α ¯ ¯ νn+1 = γ˙ 0 ¯ sign Tn+1 mαn+1 · sαn+1 H Fn+1 , α = 1, . . . , N, α ¯ ρˆn+1 ζn+1 Gn+1 =

N X

³ ´ β νn+1 sβn+1 ⊗ mβn+1 ,

(82)

β=1 α Fn+1

6.2

α | − ζα α α α = |τn+1 n+1 = |det(Fn+1 )Tn+1 mn+1 · sn+1 | − ζn+1 .

Resolution of the non-linear parabolic equations

The second update algorithm involves the diffusion-like evolution equation for dislocation densities, (41), written under the form ¯ ¯ 1 ¯¯ ζ α ¯¯n α ∂ψ ρ˙ = k∆ρα − α , α = 1, . . . , N. (83) Dν0 ¯ τ α ¯ ∂ρ for an active α− slip system. We consider the homogeneous boundary condition (42), i.e. iα (ρq ) = 0. We multiply equation (83) by an arbitrary weight function ξ, integrate over the domain Ωt , and using the divergence theorem, obtain the following equation 1 Dν0

¯ α ¯n Z Z ¯ζ ¯ α ∂ψ α ¯ ¯ ρ˙ ξ dV = − k∇ρ · ∇ξdV − ξdV, ¯τα ¯ α Ωt Ωt Ωt ∂ρ

Z

α = 1, . . . , N.

(84)

Let ραh = Neρ T Ae ρα be an approximation of ρα on the element e, ρα represents a global vector containing the nodal values of the dislocation densities, i.e. ρα = [ρα1 , ..., ραN odes ]T where N odes denotes the number of the nodes of the network, and Neρ are the shape functions given by £ ¤ Neρ T = N1e N2e N3e

23

The discretized form of (83) is given by à ! Z ¯ α ¯n X ¯ ¯ 1 ζ ¯ ¯ Neρ Neρ T dV Ae ρ˙ α = AeT ¯ ¯ e τα Dν 0 Ω t à e ! Z Z X X T eT eT eT e α eT − A k(∆Nρ ) (∆Nρ )dV A ρ − A e

or

Ã

X

Ωet

! AeT Meα (ρq )Ae

à ρ˙ α = −

e

where

Ωet

e

X

! AeT Ke Ae

ρα +

X

e

Neρ

∂ψ dV, ∂ρα

AeT fαe (ρq ),

(85)

(86)

e

¯ α ¯n Z ¯ζ ¯ 1 ¯ ¯ Neρ Neρ T dV, = Dν0 ∂Ωet ¯ τ α ¯ Z T e K = k(∆Neρ T ) (∆Neρ T )dV, Ωe Z t ∂ψ fαe (ρq ) = Neρ α dV. ∂ρ Ωet

Meα (ρq )

· We mention that ∆ is the differential operator ∆ = a differential system of the form

∂ ∂ ∂x1 ∂x2

(87) (88) (89)

¸T . The relation (86) is just

Mα (ρq )ρ˙α + Kρα = Fα (ρq ).

(90)

In order to numerically solve of the differential system (97), Crank-Nicolson’s method is applied, see Zl´amal (1977). We proceed as follows. We assume that the solution of the equation (90) is given at times t1 and t2 . For n = 2, N ∗ , ραn+1 is given by µ ¶−1 ·µ ¶ ¸ dt dt α α q α q α α q ρn+1 = M (¯ ρn+1 ) + K M (¯ ρn+1 ) − K ρn + dtF (¯ ρn+1 ) α = 1, ..., N, (91) 2 2 3 q ρ − 2 n unchanged during the ¯ qn+1 = where ρ

1 q ρ , q = 1, ..., N, while the integrals are taken on Ωt , which remains 2 n−1 time interval [tn , tn+1 ], i.e. Ωt = Ωtn

Let us introduce for any α the set of all nodes which satisfy the activation condition at time tn Mαn = {nod/Fnα = |τ α (nod)| − fc ζnα (nod) ≥ 0}

(92)

We check the activation condition and we select the solution following the procedure ½ (ραn+1 )j

=

(ραn+1 )j , if j ∈ Mαact (ραn )j , if j ∈ / Mαact

24

(93)

Finally, the update values for the Cauchy stress, mass density, deformation gradient and elastic distortion are given by (80)–(82) in which ραn+1 is computed by the procedure provided in (91). Remark. The second update algorithm allows the case k = 0, which corresponds to the diffusionless models.

7

Numerical results

As we have already remarked, in a plane stress state the slip systems become active in pairs at the same time, see formulae (51), (53) together with Table 1. Note that the slip systems considered in the examples are from the Table 1, and the normals to the planes and the slip directions could never be in the same plane. Note that the sheet and the Cartesian basis attached to the sheet have the same orientation as the crystalline lattice. The effect of the initial orientations of the slip system has not been investigated. However we do investigate the influence of the pair of slip systems potentially activated. Solving the equilibrium problem for viscoplastic crystalline materials, containing dislocations, the macro boundary conditions which traditionally describe the stress and velocity would have to be completed with micro boundary conditions. We consider that the boundary of the sheet is impenetrable to dislocations, and that we have a plane stress state during the deformation process. Compression/tension of an fcc-single crystal. In our example, a rectangular-shaped single crystal is subject to a simple compression under plane stress conditions as shown in Fig.1. In order to eliminate a rigid rotation of the sheet, a point on the sheet, say (0, 0, 0), is fixed. The macroscopic boundary conditions are assumed to be given as in Kuroda (2011)  v · i1 = 0, s˙ ∗t · i2 = 0 for X1 = 0, X2 ∈ [0, l0 ]    ∗ ∗ v · i1 = −v1 , s˙ t · i2 = 0 for X1 = L0 , X2 ∈ [0, l0 ] (94) ∗ =0 ˙ s for X2 = 0, X1 ∈ [0, L0 ]    t∗ s˙ t = 0 for X2 = l0 , X1 ∈ [0, L0 ] ∂ρα = 0 on the boundary of the sheet and a homogeneous initial distri∂n bution of dislocations, ρα (x, t0 ) = ρ0 , is considered. Moreover, we consider

Here L0 and l0 are the initial length and width of the specimen, respectively, and v1∗ is a prescribed end-displacement rate taken to be a constant and for compression v1∗ > 0 while for the tension v1∗ < 0. In the numerical application we have considered L0 = 50mm, l0 = 20mm and v1∗ = 5 · 10−3 mm/s. The material parameter values for copper are given by Teodosiu et al. (1993): ρˆ0 c¯11 = 166.1 GPa, ρˆ0 c¯12 = 119.9 GPa, ρˆ0 c¯44 = 75.6 GPa, µ = 45 GPa, γ˙ 0 = 10−3 s−1 , n = 20, ρ0 = 2730 mm−2 , b = 2.57 · 10−7 mm, yc = 0.5 · 10−6 mm, K = 75, aαβ = 0.42 if α = β, 0.52 otherwise, ρˆ0 = 8.96 · 10−3 g/mm, ζ0α = ζ0 , ζ0 = 37.7 · 10−3 GP a, i.e. for the initial critical reduced shear, ζ α (t0 ). In the numerical examples the back stress is taken to be zero. Since the axial forces applied to the specimen and the length of the specimen, L, during the deformation process can be measured, in the numerical simulation we consider curves rep-

25

Figure 1: The compression problem of a single crystal under plane stress conditions. resenting the axial stress versus the total strain, with the total strain defined by E11 =

L − L0 . L0

A convergence study has been performed in terms of the time step size, namely for dt = 0.01s and dt = 0.001s. The equivalent steps of strain at the end of the sheet are 10−6 and 10−7 , respectively. The axial stress versus the total strain is computed in the sheet, during a tensile test, when ¯ 1 ) and (¯s7 , m ¯ 7 ), using the two time steps. The effect of time two slip systems are active, (¯s1 , m steps on the axial stress- total axial strain response can be seen in the Fig.2.

80 70 dt=0.01

Axial stress (MPa)

60 dt=0.001

77.1 50 77 40

76.9 76.8

30 8.9

8.95

9 −3

20

x 10

10 0

0

0.002

0.004 0.006 0.008 Total strain ((L−L0)/L0)

0.01

0.012

Figure 2:

The effect of time steps on the stress-strain response, i.e. the axial stress versus L − L0 ¯ 1 ) and (¯s7 , m ¯ 7) total strain , in the tensile problem for k=0, when the slip systems (¯s1 , m L0 are considered. Activation conditions. From the theoretical point of view, the activation condition τ α ≥ ζ α allows us to emphasize the domains in which the plastic deformation occurs at a given moment. In the beginning of the deformation process, the sheet is entirely in an elastic state since no

26

system has been activated, namely τ α < ζ α for all finite elements. At a certain moment, say t∗ , the plastic deformation starts from the points in which the activation condition τ α ≥ ζ α is fulfilled. For the numerical computations performed the activation condition τ α ≥ fc ζ α with a certain correction factor is considered in order to ”filter out the most inactive slip systems,” see Teodosiu et al. (1993). Only if the factor fc is close to 1, or the time step is smaller (say 10−3 ), the sequential spread of the plastic zone can be put into evidence, otherwise the full sheet instantaneously enters in plastic state. In Fig.3 we illustrate the sequential spread of the incipient plastic zone which is continuously developed for the case k=0, when a non-homogeneous initial distribution of dislocation is considered, and fc = 0.96. The distribution of the dislocation density ρ1 is strongly localized in the vicinity of the source of the initial non-homogeneity, as can be seen in Figs.3(b). In this example for the sake of numerical simplicity only two slip systems are considered to be active, ¯ 1 ) and (¯s7 , m ¯ 7 ). they are (¯s1 , m ρ1 at total strain of

The plastic domain at total strain of 0.001344

0.001344

20

20

10

10

0

0

10

20

30

40

50

0

0

10

1500 20

10

10 0 20

40

0

10

2600 20

10

10

20

20

30

40

50

3400

3800

0.002315

20

0

50

2500

3000

0.002315

0

40

0.001355

20

0

30

2000

0.001355

0

20

40

0

0

10 1

20

30 2

40

50

3 4

x 10

Figure 3: The sequential spread of plastically deformed zones in tension and the distribution of ¯ 1 ) and (¯s7 , m ¯ 7 ), at various incipient the dislocation density ρ1 , for the active slip systems (¯s1 , m total strains which are mentioned below, when k = 0 (local model). Idealized slip systems. In the literature, see for instance Asaro (1979), Chang and Asaro (1981), Kuroda (2011), idealized slip systems have been considered although they are not physically possible. Asaro (1979) and Chang and Asaro (1981) introduced an idealized plane version of a double fcc- crystallographic slip systems, where both slip directions and slip plane normals lie in the plane containing the tensile axis. The authors state that ”this geometry incorporates the essential kinematical features” of actual geometry, and the incremental constitutive

27

law for the crystal has been reduced to an appropriate plane strain state. We analyze the assumptions made by Asaro (1979). The double slip mode is characterized by the fcc- crystallographic slip systems ¤ 1 £ 1 s2 = √ 1, 0, 1 , m2 = √ (1, 1, 1) , 2 3

¢ 1 1 ¡ and s5 = √ [0, 1, 1] , m5 = √ 1, 1, 1 , 2 3

(95)

in the notation adopted in the paper, with respect to the crystallographic axes, denoted by (i1 , i2 , i3 ). Let us assume that the unit vectors j1 , j2 are given in the plane generated by (s2 , s5 ), with j2 oriented along the bisectrix of the angle between the two slip directions. j3 = s2 × s5 is normal to the plane. The double slip system is characterized with respect to the axes {jk }k=1,2,3 by √ 1 3 1 2 2 1 s2 = √ j1 + j2 , m2 = √ (− √ j1 + √ j2 − √ j3 ), 2 2 3 2 6 3 √ (96) 1 3 1 2 2 1 5 5 s = − √ j1 + j2 , m = √ ( √ j1 + √ j2 − √ j3 ). 2 2 3 2 6 3 1. Thus the normals m2 and m5 have non-zero components on the direction perpendicular to the plane (j1 , j2 ), which means that the normal and slip directions of the considered slip systems (95) cannot be in the same plane. 2. Under the hypothesis that the stress state is characterized by an uniaxial tensile stress along the axis j2 , we prove that only the two slip systems (s2 , m2 ) and (s5 , m5 ) can be activated. Proof. The reduced shear stress τ α in α− slip system is calculated for the considered here case 1 by τ α = Tmα · sα , T = Tˆ22 j2 ⊗ j2 , where j2 = √ (−i1 + i2 + 2i3 ), and 6 1 1 τ 1 = τ 7 = −τ 9 = −τ 10 = − Tˆ22 , τ 2 = τ 5 = Tˆ22 , 9 6 1 τ 3 = τ 4 = −τ 6 = −τ 11 = Tˆ22 , τ 8 = τ 12 = 0. 18

(97)

Consequently, maxi∈{1,12} | τ i | is taken for i = 2, 5. 3. The rate of plastic distortion corresponding to double slip system is given by L = ν 2 s2 ⊗ m2 + ν 5 s5 ⊗ m5 , and can be expressed, for instance with respect to the basis {jk }k=1,2,3 by   1 2 1 − √ (ν 2 + ν 5 ) √ (ν 2 − ν 5 ) − √ (ν 2 − ν 5 )   6 6 3   1 Lp =  −√6(ν 2 − ν 5 ) √1 (ν 2 + ν 5 ) − √ 2 5 . (ν + ν )   6 4 3 0 0 0

(98)

(99)

We conclude that the plastic strain state does not occur in the plane (j1 , j2 ), and the shear component Lp23 has the same order in magnitude as the normal components Lp22 and Lp11 . In conclusion: a plane which contains the normal and slip directions mentioned in (95) does not exist, and if the uniaxial tensile activates the slip systems (95) then the plastic flow rule does not induce a plastic plane deformation.

28

7.1

Numerical solutions for sheets in-plane stress state

Under a homogeneous compressive (or tensile) stress along the crystalline axis i1 , the initial critical value of the reduced shear stress is achieved in eight slip systems. These activated slip systems are (i)

¯ 1 ), (¯s2 , m ¯ 2 ), (¯s4 , m ¯ 4 ), (¯s6 , m ¯ 6 ), (¯s7 , m ¯ 7 ), (¯s8 , m ¯ 8 ), (¯s10 , m ¯ 10 ), (¯s12 , m ¯ 12 ). (¯s1 , m

(100)

The result is similar with that already proved by Cleja-T ¸ igoiu (1996). Let us remark that the ¯ 2 ) is activated, while (¯s5 , m ¯ 5 ) is not. slip system (¯s2 , m As we have already noticed the slip systems become active in pairs in a plane stress state at the same time. We consider all sets of potentially active pairs (100). The rate of plastic distortion is expressed in terms of the rate of plastic shears by  1  ν − ν2 + ν6 − ν8 ν1 − ν2 − ν6 + ν8 0 ˙ p (Fp )−1 =  −ν 1 + ν 6 −ν 1 − ν 6 0 . F (101) 2 0 0 ν + ν8 Compression and tensile problems with all eight activated slip systems. We apply the above numerical algorithms to the compression and tensile tests of a crystal specimen for all eight slip systems activated in pairs, and with the dislocation densities described by local and non-local laws, respectively. The finite element mesh used in this simulation consists of 1066 elements and 585 nodes in tension, and of 954 elements, and 516 nodes in compression. The shape functions are taken to be linear. Let us denote by ρtot total dislocation density accumulated during the elasto-plastic process. The diagonal of the sheet which contains the fixed point is denoted by dr and the other diagonal is denoted by dl . We analyze the peculiar aspects of the behaviour of the sheet at the various values of the total strains E11 = 0.02, 0.03, 0.04 and 0.05. The diffusion parameter k takes the values: k = 0, and k = 5. In the tensile test the total dislocation density ρtot increases with increasing E11 from 0.02 to 0.05, and has the maximum values located on the diagonal dr for k = 0. For k = 5 the maximum values of ρtot are on the other diagonal. The dislocation density is nearly homogeneously distributed in the central zone of the sheet, keeping the same aspects for values of the total strains considered in Figs.4(a), (b) and Figs.4(c) and (d), respectively. The contour legend is different for Figs.4(c) and (d) (k = 5) since the minimum value of ρtot at E11 = 0.05 is grater than the maximum of ρtot at the previous total strain E11 = 0.04. In the tensile test simulated here, all the graphs show a skew-symmetry with respect to the diagonal, dr , in the distribution of the various fields represented. For the numerical computation performed in the tensile test, the activation condition τ α ≥ fc ζ α is considered with a correction factor fc = 0.82 to avoid the broken aspect of the finite elements which become active. p The incipient shear zones, corresponding to a maximum of F12 , start from the opposite corners of the diagonal dr in the area of the largest values of ρtot . The plastic shear zones are spread toward the central zone, being very well delimited by the lenticular shape zones with low ρtot which include the free boundaries of the sheet, see Figs.5(b) and (e). p The variation of the plastic distortion component F11 − 1 reaches an amount of 6-7 percents p p at the total axial strain of 0.05, around the corners of the diagonal dr . F12 and | F22 − 1 | remain p smaller than F11 − 1 as it can be seen in Figs.5(d), (e) and (f). The maximum variation of the

29

ρtot 20

20

10

10

0

0

0

10

4

6

20 8

30 10

40 12

50 14

0

10

4

6

20

30

8

10

40 12

50 14

6

6

x 10

(a) 20

20

10

10

0 0

10 6

20

30

6.2

40 6.4

0

50 6.6

0

10

8.8

9

20 9.2

6

(c)

x 10

(b)

30

40 9.4

50 9.6 6

x 10

(d)

x 10

Figure 4: Distribution of the total dislocation density ρtot plotted in tension for fc = 0.82, for k = 0 in (a) and (b), and for k = 5 in (c) and (d), at the total strain of 0.04 and 0.05 in (a), (c) and (b), (d), respectively. p p plastic component | F22 − 1 |, with F22 < 1, is expanded from the active corners of the diagonal dr toward the central zone, while the minimum is concentrated near the corners of the diagonal p dl as for F11 − 1. p p The similarity can be observed in the distributions of T11 and F11 − 1 and T12 and F12 , p p respectively. Although F12 and | F22 − 1 | increase significantly with the axial total strain, the stresses T12 and T22 remain very small relative to T11 , see Figs.6(c) and (f), and Figs.6(b) and (e), compared with Fig.6(a) and (d).

The maximum values of stress components T12 and T22 (that correspond to the plastic shear zones shown in Fig.5(e)) can be observed in Figs.6(b) and (e), and Figs.6(c) and (f), respectively. p Compressive tests are first performed for k = 0, here F11 < 1 and T11 < 0. For the numerical computation performed in compression, the activation condition is considered with a correction factor fc = 0.9 If the correction factor is considered to be fc = 0.82 the finite elements are smooth and activated earlier than for fc = 0.9, at the same total strain 0.03 and 0.05, compare e.g. Figs.11(c) and (d), and Figs.11(a) and (b), respectively. Moreover, the total dislocation densities remain smaller than those reached for fc = 0.9, compare Figs.8(c) and (d), and Figs.8(a) and (b), respectively, which have been plotted for k = 0 and k = 5, at the total strain 0.02.

The distribution of ρtot attains its maximum values in four extended lenticular zones, localized near the corners of the diagonal dl and in the central part of the sheet, which are separated by a zone around the diagonal dr corresponding to the minimum values of ρtot , at the total strain of 0.02 (see Fig.8(a)) and 0.03 (see Fig.7(a)). The evolution of the distribution of ρtot is

30

p F11 −1

p F12

p | F22 −1|

20

20

20

10

10

10 0

0

0 0

10

20

0.03

0.04

30

40

0.05

0.06

0

50 0.07

10

0.01

20

30

0.02

(a)

40

0

50

0.03

0.04

0.015

20

10

10

10

0.03

0.04

30

40

0.05

0.06

0

50 0.07

0.025

40

50

0.03

0.035

0

0

0 20

0.02

30

(c)

20

10

20

(b)

20

0

10

10

0.01

20

30

0.02

40 0.03

(d)

0

50 0.04

10

0.015

20 0.02

(e)

30 0.025

40

50

0.03

0.035

(f)

p p Figure 5: Variation of the components of plastic distortion, F11 − 1 ((a) and (d)), F12 ((b) and p p p (e)) and | F22 − 1 | ((c) and (f)) at the total strain 0.04, 0.05. Here F22 < 1 and F11 > 1 in tension.

T11 (M P a)

T12 (M P a)

T22 (M P a)

20

20

20

10

10

10 0

0

0 0

10

100

20

110

30 120

40 130

50

0

140

10 0

20

30

1

(a)

2

40 3

0

50 4

−1.5

5

20

10

10

10

0

0

100

110

20

30 120

(d)

40 130

50 140

20

30

−0.5

0

40 0.5

50 1

(c)

20

10

−1

(b)

20

0

10

0 0

10 0

20 1

30 2

(e)

40 3

50 4

0 5

10 −1.5

−1

20 −0.5

30

40

0

0.5

50 1

(f)

Figure 6: Distribution of the Cauchy stress components (measured in MPa), T11 (in (a) and (d)), T12 (in (b) and (e)) and T22 (in (c) and (f)) at the total tensile strain of 0.04, 0.05 and for k=0.

31

represented in Figs.7(a)-(c) for k = 0 and Figs.7(d)-(f) for k = 5. During their evolutions, two large triangular zones that correspond to the maximum values of ρtot occur at the total strain 0.05. A diagonal band with a low intensity of dislocation density is observed for k=0 at the total strain 0.05, see Fig.7(c). 20

10

0

0

10

20

1

30

2

40

3

50

20

20

10

10

0

0

10

20

1

4

30 2

40 3

0

0

1

(b)

10

10

10

1

20

30 2

40

50

0

0

10 1

3

20

30 2

(d)

40

50

0

0

10 1

3

20

3

(e)

x 10

50 4 x 10

30 2

40

50

3 7

7

7

x 10

2

40

7

20

10

30

(c)

20

0

20

7

x 10

20

0

10

4

7

x 10

(a)

50

(f)

x 10

Figure 7: Distribution of the total dislocation density ρtot is shown at the total strain of 0.03, 0.04 and 0.05 in (a), (b), (c), for k = 0 and in (d), (e) and (f) for k = 5 and fc = 0.9, in compression. p When k=0, both | F11 − 1 | and | T11 | have maximum values inside the zones corresponding to the maximum values for ρtot . These two large triangular zones are separated by a zone around p the diagonal dr , with smaller values for | F11 − 1 | and | T11 | as is shown in Figs.9(a), (d) and (g) and the corresponding Figs.10(a), (d) and (g). We remark that | T11 | is with one order grater than T12 and | T22 |, which have the same order of magnitude. The bands with alternating intensities of T22 cover nearly the full central zone and the contour plot is shown in Fig.10(c), (f) and (i). The multiple shear bands in stress with the maximum intensity of T12 appear in large zones during the compressive process at the total strain of 0.05, as shown in Fig.10(h). p The four areas with largest amount of F12 can be seen in the vicinity of the corners on the diagonal dl and in the central zone at the total strain of 0.03, in Fig.9(b), and they develop p and delimit a zone with low values of F12 and largest values of shear stress T12 . Multiple shear p zones can be seen in Fig.10(h). The variation in F12 follows the variation of ρtot , and reaches the largest magnitude of about 0.06 at the total axial strain of 0.05, during the compression process, see Fig.9(h). p The two parallel zones of maximum values of the plastic component, F22 − 1, occur in the p p zones corresponding to the relative minima of F11 − 1 and F12 , see for instance Figs.9(i), (g) and (h), respectively, and these are similarly associated with the zones of minimum values of

32

ρtot 20

20

10

10

0

0

10 1

20

30

2

40

3

0

50

4

0

5

10 1

20 2

30

40

3

4

6

20

10

10

0

10

20

1

30

2

40 3

0

50 4

0

10 1.7

20 1.8

30 1.9

6

(c)

x 10

(b)

20

0

5 6

x 10

(a)

50

40 2

50 2.1 6

x 10

(d)

x 10

Figure 8:

Distributions of the total dislocation density ρtot at the total strain 0.02, (a) and (c) for k = 0, (b) and (d) k = 5, (a) and (b) for fc = 0.9, and (c) and (d) for fc = 0.82, in compression. the stress T22 , see Fig.10(i). During the compression and tensile test various slip systems could be active or inactive (this is modeled in terms of the Heaviside function which enters the formula (40)). The distribution of the active zones, corresponding to the pair of slip systems, eight and ten, can be seen in Fig.11 at the axial total strain of 0.03 and 0.05, for k = 0. The inactive elements are characterized by shear rates ν 8 = ν 10 = 0 (the blue zones). In compression, the active elements are localized p in the four large zones with largest values of the plastic shear F12 , see Figs.11(a) and (b), for fc = 0.9, and Figs.11(c) and (d), for fc = 0.82. In the simulation of the tensile test, the two active areas are located in the lenticular zones shown in Figs.11(e) and (f) at a total train of 0.03 and 0.05, respectively.

33

p | F11 −1|

p F12

p F22 −1

20

20

20

10

10

10

0

0

10 0.02

20

30

0.04

40 0.06

50

0

0.08

0 0

10

0

20

30

0.02

(a)

0.04

40

50

0

0.06

0.01

(b) 20

20

10

10

10

0

10 0.02

20

30

0.04

40 0.06

50

0

0

10

0

0.08

20

30

0.02

(d)

0.04

40

50

0

10

10

10

0.02

20

30

0.04

40 0.06

50 0.08

0

0 0

10

20 0.02

(g)

30 0.04

40

50

0

0

0.06

0.01

20

10

10

10

10 0.02

20

30

0.04

(j)

40 0.06

50 0.08

0

0 0

10

20 0.02

0.02

40 0.03

50 0.04

20

30

0.02

40 0.03

50 0.04

30 0.04

(k)

20

30

0.02

40 0.03

50 0.04

(i)

20

0

10

(h)

20

0

30

(f) 20

10

10 0.01

0.06

20

0

0

(e)

20

0

20

(c)

20

0

10

40 0.06

50

0

0

10 0.01

20 0.02

30

40 0.03

50 0.04

(l)

p p Figure 9: Components of the plastic distortion, | F11 − 1 | ((a), (d), (g) and (j)), F12 ((b), (e),

p (h) and (k)) and F22 − 1 ( (c), (f), (i) and l) at the total strain of 0.03, 0.04, 0.05 are shown on the first three lines for k = 0, and on the last line for k = 5.

34

T11 (M P a)

T12 (M P a)

T22 (M P a)

20

20

20

10

10

10

0

0

0

10

−150

−100

20

30

−50

40 0

0

10

20

30

40

50

−20

0

20

40

60

80

50 50

(a)

0

0

−60

20

10

10

10

10

−150

−100

20

30

−50

40 0

50 50

0

0

10

20

30

40

50

−20

0

20

40

60

80

(d)

0

10

10

10

−150

−100

20

30

−50

40 0

50 50

0

0

10

20

30

40

50

−20

0

20

40

60

80

(g)

0

0

−60

20

10

10

10

10

−150

−100

20

30

−50

40 0

(j)

50 50

0

0

10

20

30

40

50

−20

0

20

40

60

80

(k)

Figure 10:

20

40

50 60

−40

20 −20

30 0

40 20

40

50 60

10 −40

20 −20

30 0

40 20

40

50 60

(i)

20

0

10

(h)

20

0

40

(f) 20

10

0

−60

20

0

0

(e)

20

0

−20

30

(c)

20

0

−40

20

(b)

20

0

10

0

0

−60

10 −40

20 −20

30 0

40 20

40

50 60

(l)

Distribution of Cauchy stress components (measured in MPa), T11 ((a), (d), (g) and (j)), T12 ((b), (e), (h) and (k)) and T22 ((c), (f), (i) and (l)) at the total strain of 0.03, 0.04, 0.05 is plotted on the first three lines for k = 0 and on the last line for k = 5.

35

20

20

10

10

0

0

10

20

30

40

0

0.5

1

1.5

2

0

50 2.5

0

10

20

30

40

0

0.5

1

1.5

2 x 10

(b)

20

20

10

0

2.5 −4

−4

x 10

(a)

50

10

0

10

20

0

30

40

1

0

50 2

0

10

0

−4

30

40

1

x 10

(c)

20

2 −4

(d)

20

20

10

10

50

x 10

0

0 0 0

10

20

30

1

40

0

50

2

0

3

10

20 1

30

40

50

2

3 −5

−5

x 10

x 10

(e)

(f)

Figure 11: Distributions of the active and inactive (i.e. for the plastic velocities ν 8 = ν 10 = 0) finite elements concerning the eighth and tenth slip systems at the total strain of 0.03 and 0.05, (a) and (b) for fc = 0.9, (c) and (d) for fc = 0.9 in compression, and (e) and (f) for fc = 0.82 in tension.

7.2

A comparison between local and non-local solutions

One of our conclusions from this numerical simulation is that in order to emphasize the differences (over 10−4 ) between the components of the plastic distortion computed for a diffusion parameter k 6= 0 and for k=0, respectively, it is necessary to have a sufficiently large diffusion parameter k (larger than 5) and the total axial strain to be at least E11 = 0.05. We discuss the role of the diffusion for k=5 at the total strain of 0.05. The distribution of the dislocation density. • The distributions of ρtot are looking completely different for k = 5 and k = 0 at various total strains as can be seen in Fig.7 and Fig.8.

36

• The dislocation density tends to become homogeneous (uniform) on large zones in the tension test, having well delimited borders for k = 5, in contrast with diffuse border of ρtot for k = 0, by comparing Figs.4(b) and (d). The maximum values of ρtot are localized in concentric zones at the upper corner on the left, and at the lower corner on the right, respectively, for k = 5, in the above mentioned figures. • In compression, the large band-zones of relative minimum and maximum values of the total dislocation density can be seen for k=0 in Fig.7(c). This dislocation density changes for the case k=5, as we can see from Fig.7(f). Two large and uniform triangular zones that correspond to the maximum values ρtot for k = 5, but smaller than the maximum for k=0, are delimited by the minimum zone displayed on diagonal dr . The magnitude of ρtot is increasing, while the diagonal zone corresponding to the minimum values of ρtot is narrowing (compare Fig.8(b) and Fig.7(d) for k=5). During the compressive deformation process, the maximum values of ρtot are significantly increasing for k = 5, being organized in plateaus, and simultaneously the small isolated zones of minimum values become joined together and develop narrow bands. • The distribution of the total dislocation density is plotted for k=5 in Fig.4(c) and (d) in tension, and in Fig.8(b) and Figs.7(d), (e) and (f) in compression. The difference of the maximum and minimum values of the total dislocation density divided to the maximum value (for k = 5) is observed to become much smaller than its value at k = 0, and corresponding to the total strain of 0.02, 0.03, 0.04 and 0.05. The variation is of 53% versus 95%, as can be seen in Table 2. The minimum value for k=5 is greater than the corresponding minimum value for k=0, while the maximum value for k=5 is smaller than the maximum value for k=0, for the total strain values grater than 0.02, see Table 2.

ρmax − ρmin , for k = 0 and k = 5 and ρmax corresponding to entire range of the total strain values. Table 2:

Variation of ρδ , defined by ρδ =

The plastic distortion and stress. Although a significant variation of the dislocation density occurs in tension and in compression (see Figs.4(c) and (d), Figs.7(e) and (f), respectively), the variations of the stress and plastic distortion remain relatively small, at the total strain 0.05 when the correction factor is 0.82. To study the influence of the diffusion parameter, k, on the plastic distortion and on the

37

stress components for a total strain of about 0.05, we considered the correction factor to be fc = 0.9 in compression and fc = 0.82 in tension. For the numerical simulations corresponding to a tension state, it was used a finer mesh than that corresponding to a compression state in order to control as much as possible the numerical perturbations generated by the material point that has been fixed, whilst the correction factor was taken to be 0.82 to ensure the continuity of the finite elements which become simultaneously active. • The local zones with maximum values of the stress components T11 and T12 tend to unify in a larger zone for k = 5, in compression. The non-homogeneous zones with the alternating tension and compression values for the stress component T22 , see Figs.10(l) and (i), and the shear band-zones corresponding to T12 , see Figs.10(k) and (h), emphasize the tendency to join together and develop well organized large bands. To emphasize the influence of the diffusion parameter we crop small regions from the main contour plots, Figs.9(g) and (h) for k = 5, and Figs.9(j) and (k) for k = 0, respectively. We show the crops in Fig.12. Zones corresponding p to the local maximum values of | F11 − 1 | tend to unify in a larger zone and simultaneously p the difference between the local maximum values is decreasing.The minimum value of F12 is p increasing for k = 5 and a zone with larger values for F12 than previous one is developing. The p zone of local maximum values for F12 is spread toward the zone with minimum value and the curves that separate the regions become smooth.

8

6

4

2 20

8

0.08

22

24

26

28

30

0.08

0.06

6

0.06

0.04

4

0.04

0.02

2 20

22

24

(a)

26

28

30

(b)

8

8

0.06 0.05

6

0.06 0.05

6

0.04

0.04

0.03

4

0.03

4

0.02 2 20

22

24

26

0.02

28

0.02 2 20

30

(c)

22

24

26

28

30

(d)

p p Figure 12: Crops of the contour plots for |F11 − 1| ((a) and (b)) and F12 ((c) and (d)) at the

total strain of 0.05; (a) and (c) for k = 0 and (b) and (d) for k = 5, compare with Fig.9.

We crop the main contour plots of the stress components T11 and T12 , Figs.10(g) and (h) for k=0, and Figs.10(j) and (k) for k=5. In Figs.12(a) and (b) we show that the alternating bands

38

with tensile and compressive values for T11 tend to join together and develop a large zone with smaller difference between the local maximum values of T11 . The band-zone of minimum values for the shear stress T12 , and which separates the zones of relative maximum values, becomes smaller for k = 5, than the corresponding zone for the diffusion-less case, as can be seen in Figs.13(d) and (c).

18

0

16

−50

14 −100 12

18

0

16

−50

14 −100 12

30

35

40

−150

30

(a)

35

40

(b)

80

20

20

80

60

18 16

20

14

0 30

35

60

18

40

25

−150

40 16

20 0

14

−20

25

(c)

30

35

−20

(d)

Figure 13: Crops of the contour plots for T11 ((a) and (b)) and T12 ((c) and (d)) at the total strain of 0.05; (a) and (c) for k = 0 and (b) and (d) for k = 5, compare with Fig.10. • To emphasize the role of the diffusion in the tensile test, we have computed the algebraic differences between the components of the plastic distortion and stress for k = 5 and k = 0. Starting from these computed algebraic differences, in Fig.14 we present the values of the rescaled differences superposed on the values of the field for k = 0. The following conclusions can be drawn: Well marked zones can be observed in Figs.14(d), (e) and (f), in the distribution of Cauchy stress components. The local zones with maximum values of T11 tend to unify in a larger zone for k = 5. The shape of the zone that corresponds to the maximum intensity of T11 is perturbed toward the zone that corresponds to the maximum of ρtot , see Figs.14(d) and 4(d). Three well localized zones can be seen in the distribution of T12 , Fig.14(e), a central shear zone and two ones containing obliquely oriented shear bands with relatively low and high stress intensities (k=5). The contour plots for T22 emphasize the tendency to form non-homogeneous zones with the alternating tension and compression values. The stress components T12 and T22 remain very small comparing with T11 , and their distribution has a broken aspect. • The most significant influence of the diffusion can be seen in the shear component of the

39

20

20

20

10

10

10 0

0

0 0

10 0.03

20 0.04

30 0.05

40 0.06

0

50 0.07

10

20

0.01

30

0.02

(a)

40

0

50

0.03

0.04

0.015

20

10

10

10

0

0

100

110

20

30 120

(d)

0.02

40 130

50 140

30 0.025

40

50

0.03

0.035

(c)

20

10

20

(b)

20

0

10

0 0

10 0

20 1

30

40

50

3

4

5

2

(e)

0

10 −1.5

−1

20 −0.5

30

40

0

0.5

50 1

(f)

p p p Figure 14: The values of plastic distortion components F11 −1, F12 and | F22 −1 | with rescaling

difference and associated with k = 5 are plotted at the total strain of 0.05 in (a), (b) and (c), p p respectively. Here F22 < 1 and F11 > 1 in tension. The values of stress components T11 , T12 and T22 , calculated for the rescaling difference and associated with k = 5 are plotted at the total strain of 0.05 in (d), (e) and (f), respectively. p plastic distortion F12 , in tensile test. A new branch is developed from the zone of maximum plastic shear in the direction of the maximum of ρtot , see Fig.14(b) and Fig.4(d). A succession of zones with low and relatively large intensity occurs from the maximum zone for F22 in the direction of maximum for ρtot .

Tensile problem with non-homogeneous initial dislocation density. Consider now the non-local evolution equation (41) with the boundary conditions considered to be the impenetrability of the sheet borders to dislocations, when a non-homogeneous initial distribution of dislocations is given by the function  ¡ ¡ ¢¢ 1  − 2 ln(ρmin /ρmax ) (X1 − X10 )2 + (X2 − X20 )2 − R2 ,  ρmax · exp R ρα0 (X1 , X2 ) = (102) for (X1 − X10 )2 + (X2 − X20 )2 ≤ R2 ,   ρmax , for (X1 − X10 )2 + (X2 − X20 )2 > R2 , ¯ 1 ), (¯s7 , m ¯ 7 ), are considered to be active. For when only two slip systems, namely (ii) (¯s1 , m the sake of simplicity we restricted the analysis to this case (with only two slip systems). The initial dislocation density ρα0 is considered constant outside the circle centered at (X10 , X20 ) and of radius R. The minimum dislocation density occurs in the center of afore mentioned circle, i.e. ρ10 (X10 , X20 ) = ρmin . In the numerical applications, X10 = L0 /2, X20 = l0 /2, ρmax = ρ0 , ρmin = 1400mm−2 , R = 1.5mm and R/L0 = 0.03. In addition, the diffusion constant k is chosen to be 0, 0.0005, 0.005 and 0.05, respectively. The results show that the small initial non-homogeneity in the distribution of the dislocation density disappears in the solution to the diffusion-like evolution equation of dislocation density

40

5

ρ1

x 10

k=0.05

1.9

5

x 10

k=0.005

1.85

1.85

1.8 1.8 k=0.0005

1.75

1.75

1.7 k=0

1.65

1.7

1.6 20

10

1.65 0

10

20

30

40

50

Figure 15: Distribution of the dislocation density ρ1 at the total strain of 0.01 for the numerical values of the diffusion parameter k = 0, k = 0.0005, k = 0.005 and k = 0.05. (for k non-zero). Moreover, for these cases, a nearly homogeneous distribution of the density of dislocation can be observed in Fig.15, at a total strain of E11 = 0.01. By contrast, for k =0, the diffusionless evolution takes place, and the distribution of the dislocation density remains nonhomogeneous, with a large amount in the vicinity of the initial non-homogeneity, still visible (like a small crater )(see also Fig.3). The shape of the plot for a small values of k, (e.g.k = 0.0005,) shows am evolution of the distribution of dislocation density from the case k = 0 to larger values, say k = 0.005 or k = 0.05, at small values of the total strain. The aspect of the graph for larger values of the total strain (E11 = 0.01), even for k small (k=0.0005) is drastically changed and becomes similar to the plots obtained for small values of the total strain but bigger values of k (k = 0.005 and k = 0.05), see Fig.15. Otherwise said, the effect of the diffusion is clearly showed, as for larger values of the total strain the response of the sheet ” forgets” the existence of the initial inhomogeneity in the distribution of dislocation density. We remark the concordance of the Figs.15 and 3. Conclusions In this work, a novel numerical algorithm has been developed for the stable and accurate numerical solutions of the initial and boundary value problems, resulting from the simulation of in compression/tension experiments of a sheet of metal. The numerical simulations prove the efficiency of the proposed numerical algorithms for the analysis of the deformation under a plane stress state, of the sheet made up from an fcc-single crystal, in which different fcc-slip systems are simultaneously activated. Compressive and tensile tests were solved by considering all eight activated slip systems, together with the activation conditions. We do not considered that once a slip system was activated it remains active for the rest of the simulation. Specifically, we focus on the mobility of the passage from active to inactive slip systems in diffusion and non-diffusion evolution equation for dislocation density. The diffusion parameter k involves a material length scale. As the physical interpretation of the k- parameter is still an open problem, a numerical study of the influence of k on the flow patterns has been conducted this. The influence of the diffusion parameter k is major

41

on the distribution of the total dislocation density. The distribution of the total density ρtot of dislocations tends to become uniform in certain localized zones, the plateaus of maximum values are criss-crossed by thin bands of minimum values. The influence of the diffusion is more significant in compression than in tension. The local zones with maximum values of plastic distortion components tend to unify in a larger zone for k = 5, in compression. The nonhomogeneous band-zones with the alternating maximum and minimum values for the stress components emphasize the tendency to join together and develop well organized large bands. The relative differences between local maximum and minimum values of the fields are decreasing with increasing k. The curves that separate the various regions become more smooth and the regions become more diffuse with increasing k. In particular, we analyzed a simple example in which only two active slip systems are considered, and when an initial non-homogeneity of the dislocation density is considered. For this case, we characterize the influence of the diffusion parameter k, on the distribution of the dislocation density. Future work would have to include a systematic numerical study for in compressive/ tensile tests considering all activated slip systems, and with the non-homogeneous initial distribution of the dislocation density, as well as very large total strain, (say of 0.1 -0.2). Acknowledgments: The authors acknowledge the support from the Ministry of Education Research and Innovation under CNSIS PN2 Programme Idei, PCCE, Contract No. 100/2009. The very constructive comments and suggestions made by the anonymous reviewers are gratefully acknowledged. References Acharya, A., Bassani, J.L., 2000. Lattice incompatibility and gradient theory of crystal plasticity. J. Mech. Phys. Solids 48, 1565–1595. Arsenlis, A., Parks, D.M. 1999, Crystallographic aspects of geometrically-necessary and statistically-stored dislocation density. Acta Mater. 47, 1597–1611. Arsenlis, A., Parks, D.M., Becker, R., Bulatov, V.V., 2004. On the evolution of crystallographic dislocation density in non-homogeneously deforming crystals. J. Mech. Phys. Solids 52, 1213–1246. Asaro, R.J. 1979, Geometrical effects in inhomogeneous deformations of ductile single crystal, Acta Matallurgica, 27, 445-453. Asaro, R.J., Needleman, A., 1985. Texture development and strain hardening in rate dependent polycrystals. Acta Metall. 33, 923–953. Bortoloni, L., Cermelli, P., 2004. Dislocation Patterns and Work-Hardening in Crystalline Plasticity. J. Elasticity, 76, 113-138. Clayton, J.D., Hartley, C., 1985. Texture development and strain hardening in rate dependent polycrystals. Acta Metall. 33, 923–953. Clayton, J.D., Hartley. C.S., McDowell, D.L., 2014. The missing term in the decomposition of finite deformation Int. J. Plasticity, 52, 51-76

42

Chang, Y.W., Asaro, R.J., 1981. An experimental study of shear localization in aluminumcopper single crystals. Acta Metallurgica, 29, 241–257. Cleja-T ¸ igoiu, S., So´os, E., 1990. Elastoplastic models with relaxed configurations and internal state variables. Applied Mechanics Reviews, 131–151. Cleja-T ¸ igoiu, S., 1996. Bifurcations of homogeneous Deformations of the bar in finite elastoplasticity. European Journal of Mechanics - A/Solids, 15, 761–786. Cleja-T ¸ igoiu, S., 2000. Anisotropic and dissipative finite elasto-plasticity. Rendiconti del seminario matematico, Univ. Pol. Torino, 58, 1, 69–82. Cleja-T ¸ igoiu, S, Matei A., 2012. Rate boundary value problems and variational inequalities in rate-independent finite elasto-plasticity. Mathematics and Mechanics of Solids, 17 (6), 557-586. Cleja-T ¸ igoiu, S, 2013. Non-Local elasto-viscoplastic models with dislocations in finite elastoplasticity. Part I: Constitutive framework, Mathematics and Mechanics of Solids, 18, 349–372. Cleja-T ¸ igoiu, S, Pa¸scan, R., 2013. Non-Local elasto-viscoplastic models with dislocations in finite elasto-plasticity. Part II: Influence of dislocations in crystal plasticity, Mathematics and Mechanics of Solids, 18, 373–396. Cermelli, P., and Gurtin, M.E.,2001. On the characterization of the geometrically necessary dislocations in finite plasticity. J. Mech. Phys., 49, 1539–1568. Ortiz Cuiti˜ no, A., Ortiz, M., 1993. The hardening of single crystal in Large Plastic Deformations, Fundamental Aspects and Applications to Metal Forming, eds. Teodosiu, Raphanel and Sidoroff, Brook-fiels, Rotterdam, 39–51. Delaire, F., Raphanel, J.L., Rey, C, 2000, Plastic heterogeneities of a copper multicrystal deformed in uniaxial tension: experimental study and finite element simulations. Acta Materialia 48, 1075–1087. Evers L.P., Parks, D.M., Brekelmans, W.A.M., Geers, M.G.D., 2002. Crystal plasticity model with enhanced hardening by geometrically necessary dislocation accumulation. J. Mech. Phys. Solids, 50, 2403–2424. Franciosi, P., 1985, The concepts of latent hardening and strain hardening in metallic single crystals, Acta Metallurgica, 9, 1601-1612. Gurtin, M.E., 2002. A gradient theory of single-crystal viscoplasticity that account for geometrically necessary dislocations. J. Phys. Solids, 50, 5–32. Gurtin, M.E., 2003. On a framework for small-deformation viscoplasticity: free energy, microforces, strain gradients. Int. J. Plasticity 19 47–90. Inal, K., Wu, P.D., Neale, K.W., Instability and localized deformation in polycrystalline solids under plane-strain tension, 2002. Int. J. Solid Struct. 39, 983–1002. Ha, S., Kim, K., 2011. Heterogeneous deformation of Al single crystal: Experiments and finite element analysis. Mathematics and Mechanics of Solids, 16, 651-661.

43

Homayonifar, M., Mosler, J., 2011. On the coupling of plastic slip and deformation-induced twinning in magnesium: A variationally consistent approach based on energy minimization. Int. J. Plasticity 27, 983–1003. Hughes, T.J.R., 1987. The Finite Element Method. Linear Static and Dynamic Finite Element Analysis, Prentice-Hall, Inc., Englewood Cliffs, New Jersey 07632. Keller, C., Hug, E., Habraken, A.M., Duchene, L., 2012. Finite element analysis of the free surface effects on the mechanical behavior of thin nickel polycrystals. Int. J. Plasticity 29, 155-172 Klusemann, B., Svendsen, B., Vehoff, H., 2013. Modeling and simulation of deformation behavior, orientation gradient development and heterogeneous hardening in thin sheets with coarse texture. Int. J. Plasticity 50, 109–126. Klusemann, B., Yal¸cinkaya, T., 2013. Plastic deformation induced microstructure evolution through gradient enhanced crystal based on a non-convex Helmholtz energy. Int. J. Plasticity 48, 168-188. Van Houtte, P., Li, S., Seefeldt, M., Delannay, L., 2005, Deformation texture prediction: from the Taylor model to advanced Lamel model. Int. J. Plasticity 21, 589–624. Kuroda, M., 2004, A phenomenological plasticity model accounting for hydrostatic stresssensitivity and vertex-type of effect. Mechanics of Materials 36, 285–297. Kuroda, M., 2011. On large-strain finite element solutions of higher-order gradient crystal plasticity, 2011. Int. J. Solid Struct. 48, 3382–3394. Kuroda M., Tvergaard, V., 2008. A finite deformation theory of higher-order gradient crystal plasticity. J. Mech Phys. Solids 56, 2573–2584. Ling X., Horstemeyer M.F., Potirniche, G.P., 2005. On the numerical implementation of 3D rate-dependent single crystal plasticity formulations, Int. J. Numer. Meth. Engng. 63, 548–568. Man, C-S, Gao, X., Godefroy, S., Kenik, E.A., 2010. Estimating geometric dislocation densities in polycrystalline materials from orientation imaging microscopy. Int. J. Plasticity 26, 423–440. Mandel, J., 1972. Plasticit´e classique et viscoplasticit´e, CISM-1971, Udine, Springer-Verlag, Vienna, New York. Mayeur, J.R., McDowell, D.L., D.J. Bammann, D.J., 2011. Dislocation-based micropolar single crystal plasticity: Comparison of multi and single criterion theories. J. Mech. Phys. Solids 59, 398–422. McMeeking, R.M., Rice, J.R., 1975, Finite-Element Formulations for Problems of Large Elastoplastic Deformations, International Journal of Solids and Structures 11, 601-616. Mecking, H., Kocks,U., 1981. Kinetics flow, strain hardening. Acta Metallurgica 29, 1865-1875. Needleman, A., Tvergaard, V., 1984. Finite element analysis of localization in plasticity. In: Oden, J.T., Carrey, G.F. (Eds.), Finite elements- Special Problems in Solid Mechanics.

44

Prentice Hall, Englewood Cliffs, pp. 94–157. Ohashi, T., 2005, Crystal plasticity analysis of dislocation emission from micro voids. Int. J. Plasticity, 21, 2071-2088. Ortiz, M., Stainier, L., 1999. The variational formulation of viscoplastic constitutive updates. Comput. Methods Appl. Mech. Eng. 171, 419-444. ¨ Oztop, M.S., Niordson, C.F., Kysar, J.W., 2013, Length-scale effect due to periodic variation of geometrically necessary dislocation densities. Int.J. Plasticity 41, 189–201. Peirce, D., Asaro R.J., Needleman, A., 1983. Material rate dependence and localized deformation in crystalline solids. Acta metall. 31 1951–1976. Svendsen, B., 2002, Continuum thermodynamic models for crystal plasticity including the effects of geometrically-necessary dislocations, J. Mech. Phys. Solids, 50, 1297–1329. Teodosiu, C., Raphanel, J.L., Tabourot, L., 1993. Finite element simulation of the large elastoplastic deformation of multicrystals, in Large Plastic Deformations, Fundamental Aspects and Applications to Metal Forming, eds.Teodosiu, Raphanel and Sidoroff, Brook-fiels, Rotterdam, 153–168. Teodosiu, C., Sidoroff, F., 1976. A theory of finite elastoviscoplasticity of single crystals. Int. J. Engn. Sci. 14 (1976), p. 165–176. Ting, T.C.T., 1996, Anisotropic Elasticity, Oxford University Press. The MathWorks, Inc., 3 Apple Hill Dr., Natick, MA, Matlab. Triantafyllidis, N., Needleman, A., Tvergaard, V., 1982, On the development of shear bands in pure bending. Int. J. Solids Struct. 18, 121-138. Truesdell, C., Noll, W., 2004, The Non-Linear Field Theories of Mechanics, ed. Antman, S.S., Third Edition, Springer. Segerlind, L.J., 1984. Applied Finite Element Analysis, John Wiley and Sons, New York. Simo, J.C., Hughes. T.J.R., 2000, Computational Inelasticity, Springer. Zl´amal, M., 1977. Finite element methods for nonlinear parabolic equations, RAIRO-Analyse num´erique, tome 11, n0 1, p.93–107.

Appendix A. Proof in the Theorem 2. We look at the rate type constitutive equation written in (23) projected on i3 . First we consider the following expression µ ¶ d T13 Tk3 T1k ˙ p (Fp )−1 )(Fe )−1 }S )kl A.1 = L1k + L3k + E13kl Dkl − E13kl ({Fe (F dt ρˆ ρˆ ρˆ ³ ´ − Fe {F˙ p (Fp )−1 }(Fe )−1

1k

´ Tk3 T1k ³ e ˙ p p −1 − F {F (F ) }(Fe )−1 . ρˆ ρˆ 3k

45

We analyze, term by term, the expressions in equation (A.1) L1k

Tk3 T13 T23 = L11 + L12 , ρˆ ρˆ ρˆ

T1k T13 L3k = L33 ρˆ ρˆ

(A.2)

E13kl Dkl = E1313 D13 + E1323 D23 = 0 ˙ p (Fp )−1 }(Fe )−1 )kl = E13kl (Fe {F ˙ p (Fp )−1 }(Fe )−1 )13 + E1323 (Fe {F˙ p (Fp )−1 }(Fe )−1 )23 = 0. E1313 (Fe {F The last two terms in equation (A.1) can be written as ³ ´ ˙ p (Fp )−1 }(Fe )−1 Fe {F

1k

Tk3 ρˆ

³ ´ ˙ p (Fp )−1 }(Fe )−1 = Fe {F

11

T13 + ρˆ

³ ´ ˙ p (Fp )−1 }(Fe )−1 + Fe {F

12

T23 , (A.3) ρˆ

´ ´ T1k ³ e ˙ p p −1 T13 ³ e ˙ p p −1 F {F (F ) }(Fe )−1 = F {F (F ) }(Fe )−1 ρˆ ρˆ 3k 33 µ ¶ d T13 Similar transformations can be derived if we look at expression of . dt ρˆ Thus we can conclude that the linear differential system which allows us to derive the expression of (T13 , T23 ) has the following form µ ¶ d T13 T13 T23 = B11 + B12 (A.4) dt ρˆ ρˆ ρˆ d dt

µ

T23 ρˆ

¶ = B21

T13 T23 + B22 , ρˆ ρˆ

where {Bij }ij are supposed to be continuous functions, dependent of the state of the deformation of the material. In hypothesis that at the initial moment (T13 (t0 ), T23 (t0 )) = 0, only the trivial solution is obtained, namely (T13 = 0, T23 = 0). Appendix B. Elastic matrix in the deformed configuration The elastic moduli are expressed in terms of the appropriate cubic elastic constants, c¯11 , c¯44 , c¯12 , and the components of elastic distortion as follows

46

e )4 + (F e )4 )¯ e )2 (F e )2 c e 2 e 2¯ E 11 ≡ E1111 = ((F11 c11 + 2(F11 44 12 12 ¯12 + 4(F11 ) (F12 ) c

(B.1)

e )4 + (F e )4 )¯ e )2 (F e )2 c e 2 e 2¯ , E 22 ≡ E2222 = ((F21 c11 + 2(F21 44 22 22 ¯12 + 4(F21 ) (F22 ) c e )2 (F e )2 + (F e )2 (F e )2 )¯ e F e F e F e )¯ E 12 ≡ E1122 = ((F11 c11 + (4F11 21 12 22 12 21 22 c44 e )2 (F e )2 + (F e )2 (F e )2 )¯ +((F11 c12 . 22 12 21

e )2 ((F e )2 + (F e )2 )¯ E 13 ≡ E1133 = (F33 c12 , 11 12

(B.2)

e )3 F e + (F e )3 F e )¯ e 2 e e e 2 e e c E 14 ≡ E1112 = ((F11 12 21 12 22 c11 + ((F11 ) F12 F22 + (F12 ) F11 F21 )¯ e )2 F e F e + F e (F e )2 F e )¯ +2((F11 12 22 11 12 21 c44 .

e )2 ((F e )2 + (F e )2 )¯ E 23 ≡ E2233 = (F33 c12 , 21 22

(B.3)

e )3 F e + (F e )3 F e )¯ e 2 e e c E 24 ≡ E2212 = ((F21 44 11 22 12 c11 + (F21 ) F12 F22 )¯ e )2 F e F e + (F e )2 F e F e )¯ e e 2 e +((F21 12 22 22 11 21 c12 + 2(F21 (F22 ) F11 , e )4 c E 33 ≡ E3333 = (F33 ¯11 , e )2 (F e F e + F e F e )¯ E 34 ≡ E3312 = (F33 11 21 12 22 c12 , e )2 (F e )2 + (F e )2 (F e )2 )¯ e Fe Fe Fe c E 44 ≡ E1212 = ((F11 c11 + 2F11 21 12 22 21 12 22 ¯12 e )2 (F e )2 + 2F e F e F e F e + (F e )2 (F e )2 c +((F11 22 11 22 12 21 12 21 ¯44 .

e )2 ((F e )2 + (F e )2 )¯ E 55 ≡ E1313 = (F33 c44 (B.4) 11 12 e )2 ((F e )2 + (F e )2 )¯ E 66 ≡ E2323 = (F33 c44 21 22 e )2 ((F e F e + F e F e )¯ E 56 ≡ E1323 = (F33 11 21 12 22 c44 .

47

List of the figures: 1. The compression problem of a single crystal under plane stress conditions. 2. The effect of time steps on the stress-strain response, i.e. the axial stress versus total L − L0 ¯ 1 ) and (¯s7 , m ¯ 7 ) are strain in the tensile problem for k=0, when the slip systems (¯s1 , m L0 considered. 3. The sequential spread of plastically deformed zones in tension and the distribution of the ¯ 8 ) and (¯s10 , m ¯ 10 ) at various incipient dislocation density ρ1 , for the active slip systems (¯s8 , m total strains which are mentioned below, when k=0 (local model). 4. Distribution of the total dislocation density ρtot plotted in tension for k = 0 and k = 5, at the total strain of 0.04 and 0.05 in (a), (b) and (c), (d), respectively, for fc = 0.82. p p 5. Variation of the components of plastic distortion, F11 − 1 ((a) and (d)), F12 ((b) and (e)) p p p and | F22 − 1 | ((c) and (f)) at the total strain 0.04, 0.05. Here F22 < 1 and F11 > 1 in tension.

6. Distribution of the Cauchy stress components (measured in MPa), T11 ((a) and (d)), T12 ((b) and (e)) and T22 ( (c) and (f)) at the total tensile strain of 0.04, 0.05, and for k=0. 7. Distribution of the total dislocation density ρtot is shown at the total strain of 0.03 and 0.05 in (a), (b), (c), for k = 0 and in (d), (e) and (f) for k = 5 and fc = 0.9, in compression. 8. Distributions of the total dislocation density ρtot at the total strain 0.02, (a) and (c) for k = 0, (b) and (d) k = 5, (a) and (b) for fc = 0.9, and (c) and (d) for fc = 0.82, in compression. p p 9. Components of the plastic distortion, | F11 − 1 | ((a), (d), (g) and (j)), F12 ((b), (e), (h) p and (k)) and F22 − 1 ( (c), (f), (i) and l) at the total strain of 0.03, 0.04, 0.05 are shown on the first three lines for k = 0, and on the last line for k = 5.

10. Distribution of Cauchy stress components (measured in MPa), T11 ((a), (d), (g) and (j)), T12 ((b), (e), (h) and (k)) and T22 ((c), (f), (i) and (l)) at the total strain of 0.03, 0.04, 0.05 is plotted on the first three lines for k = 0 and on the last line for k = 5. 11. Distributions of the active and inactive (i.e. for the plastic velocities ν 8 = ν 10 = 0) finite elements concerning the eighth and tenth slip systems at the total strain of 0.03 and 0.05, (a) and (b) for fc = 0.9, (c) and (d) for fc = 0.9 in compression, and (e) and (f) for fc = 0.82 in tension. p p 12. Crops of the contour plots for |F11 − 1| ((a) and (b)) and F12 ((c) and (d)) at the total strain of 0.05; (a) and (c) for k = 0 and (b) and (d) for k = 5, compare with Fig.9.

13. Crops of the contour plots for T11 ((a) and (b)) and T12 ((c) and (d)) at the total strain of 0.05; (a) and (c) for k = 0 and (b) and (d) for k = 5, compare with Fig.10.

48

p p p 14. The values of plastic distortion components F11 − 1, F12 and | F22 − 1 | with rescalling difference and associated with k = 5 are plotted at the total strain of 0.05 in (a), (b) and (c), p p respectively. Here F22 < 1 and F11 > 1 in tension. The values of stress components T11 , T12 and T22 , calculated for the rescalling difference and associated with k=5 are plotted at the total strain of 0.05 in (d), (e) and (f), respectively.

15. Distribution of the dislocation density ρ1 at the total strain of 0.01 for the numerical values of the diffusion parameter k=0, k=0.0005, k=0.005 and k=0.05.

49

A rate-type viscoplastic model with non-local evolution equations for dislocations and multislip flow rules is considered. Conditions to ensure the in-plane stress for material with cubic elastic symmetry. The variational equality for the velocity field is associated with the incremental boundary value problem at time t. The finite element method for solving the variational problem, coupled with an update algorithms. Tests accounting for all eight activated slip systems, the influence of the diffusion equation for dislocations are analyzed.