ARTICLE IN PRESS Journal of the Mechanics and Physics of Solids 58 (2010) 844–859
Contents lists available at ScienceDirect
Journal of the Mechanics and Physics of Solids journal homepage: www.elsevier.com/locate/jmps
Dynamic crystal plasticity: An Eulerian approach Oana Cazacu a,, Ioan R. Ionescu b a b
Department of Mechanical and Aerospace Engineering, University of Florida, REEF, 1350 N. Poquito Rd, Shalimar, FL 32579, USA LPMTM, University Paris 13, 99 Av. J.-B. Clement, 93430 Villetaneuse, France
a r t i c l e i n f o
abstract
Article history: Received 22 August 2009 Received in revised form 16 February 2010 Accepted 3 April 2010
In this paper an Eulerian rate-dependent single crystal model that accounts for highstrain rates, large strains and rotations is developed. The viscoplastic law as well as the evolution equations for the lattice are written in terms of vectorial and tensorial quantities associated with the current configuration. The viscoplastic law is obtained from Schmid law using an overstress approach. Such an expression for the viscoplastic law is motivated by the microdynamics of crystal defects. A general analysis of the plane-strain response of the proposed rigid-viscoplastic single crystal model is presented. It is shown that only one differential equation, involving the orientation of one composite in-plane slip system, is necessary to describe the lattice evolution. Several two-dimensional boundary value problems, such as equal-channel die extrusion and channel die compression are selected to illustrate the predictive capabilities of the model. The results show that even at relatively low strain rates the viscosity plays an important role in the development of localized deformation modes. At high crosshead velocity, the plastic properties and crystal anisotropy are less important while inertia effects are dominant. Finally, the grains interaction is investigated by analyzing the compression of a grains multicrystal. & 2010 Elsevier Ltd. All rights reserved.
Keywords: Rate dependent crystal plasticity Eulerian crystal lattice kinematics Equal-channel die extrusion Channel die compression
1. Introduction The fundamental importance of single crystal modeling is reflected in the abundance of theories on monocrystalline plasticity. The description of plastic flow by crystallographic slip in metals with cubic crystals structure dates back to the pioneering work of Taylor and Elam (1923, 1925), and Taylor (2003a, b), that established the experimental laws that are the foundation of crystal mechanics. A general kinematics of the finite deformation of elastic–plastic single crystals that includes lattice distortion was first given by Rice (1971) and may be found in equivalent forms in Hill and Rice (1972) and Asaro and Rice (1977). An overview of constitutive models for single crystals and polycrystals subjected to arbitrarily large strains can be found in Asaro and Needleman (1985), Hafner (1992), Schmidt-Baldassari (2003), Dawson (2004), Van Houtte et al. (2004), etc. In some processes, such as metal cutting, extrusion or penetration, large deformations, high pressures, and high strain rates may occur. That is why in the mechanical and numerical modeling, an Eulerian description seems to be most suitable. The main ingredients of an Eulerian model for large strains dynamic deformation of single crystals are the Eulerian equations describing the lattice rotations and the rate-dependent flow rule. However, in current rate-dependent models an finite element (FE) Lagrangean framework is generally adopted (e.g. Ling et al., 2005; Rousselier and Leclerc, 2006). Thus,
Corresponding author.
E-mail addresses:
[email protected]fl.edu (O. Cazacu),
[email protected] (I.R. Ionescu). 0022-5096/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.jmps.2010.04.001
ARTICLE IN PRESS O. Cazacu, I.R. Ionescu / J. Mech. Phys. Solids 58 (2010) 844–859
845
the large lattice rotations that occur in the dynamic flow regime may not be accurately captured. As concerns the flow rule, a Norton-type power-law is generally used. This law is adequate for the description of low strain rate behavior. Because Norton law is very stiff, current rate-dependent formulations may predict unrealistic slip rates (see for example, McGinty and McDowell, 2006; Busso and Cailletaud, 2005, and further discussion in Section 3 of this paper). The main goal of this paper is to develop an Eulerian rate-dependent single crystal model that accounts for high-strain rates, large strains, and rotations of the material and the lattice. Specific objectives include the Eulerian description of the crystal lattice kinematics and a viscopastic extension of the Schmidt law. Although the proposed model and associated numerical strategy are fully 3-D, of particular interest is the in-plane response. Since for in-plane deformation, the lattice orientation can be described by one angle, interpretation of the numerical results is largely facilitated. Furthermore, analytical treatments based on composite slip systems allow obtaining a reduced, 2-D form of the model. Besides providing a deeper understanding of the orientation flow, the 2-D form is computationally much less costly than the full 3-D. The outline of the paper is as follows. First, the viscoplastic law as well as the evolution equations for the lattice are written in terms of vectorial and tensorial quantities associated with the current configuration (Sections 2 and 3). Specifically, evolution equations for the slip directions and slip planes, involve material time derivatives. The viscoplastic law is obtained from Schmid law using an overstress approach. By adopting such a law, the numerical instabilities that arise in the integration of the classical Norton law are eliminated. The expression of the proposed viscoplastic law is motivated by the microdynamics of crystal defects (see the analysis of Teodosiu and Sidoroff, 1976). Next, Eulerian numerical methods, developed in Cazacu and Ionescu (2010) to solve boundary value problems, are briefly presented (Section 4). The numerical strategy consists in using the finite element method for solving for the velocity field and finite volume method to solve for the hyperbolic equations that describe the evolution of the orientations of the slip systems. In Section 5, we present a detailed analysis of the plane-strain response of the proposed rigid-viscoplastic single crystal model. It is shown that only one differential equation, involving the orientation of one composite in-plane slip system, is sufficient to fully describe the lattice evolution. Furthermore, for an FCC single crystal, analytic expressions for the slip rates on individual in-plane composite systems are obtained (Section 6). Several two-dimensional boundary value problems are selected to illustrate the predictive capabilities of the mechanical model. Specifically, in Section 7, is presented the application of the model to equal-channel die extrusion. We analyze the stationary flow of a FCC single crystal and the extrusion of a grain embedded in a parent crystal. Section 8 concerns channel die compression of an FCC single crystal. The influence of viscosity and of the inertia are explored. It is found that even at relatively low strain rates the viscosity plays an important role in the development of localized deformation modes. At high crosshead velocity, the plastic properties and crystal anisotropy are less important while inertia effects are dominant. Finally, we investigate the grains interaction by analyzing the compression of a four grains multicrystal. 2. Eulerian description of the lattice rotations Consider a single crystal at time t =0, free of any surface tractions and body forces and let choose this configuration, say K0 as reference configuration of the crystal. Let K ¼ KðtÞ denote the current configuration. In Fig. 1 are shown the geometric decomposition and mapping of the crystal deformation and rotation of the lattice axes. Since in applications involving large deformations and high strain rates, the elastic component of the deformation is small with respect to the inelastic one, it can be neglected and a rigid-viscoplastic approach will be adopted (such a hypothesis is generally used, e.g. Hutchinson, 1976; Lebensohn and Tome´, 1993; Kok et al., 2002, etc.). The gradient P ¼ PðtÞ, due to the dislocation glide on the slip ~ with no change in systems, maps the reference crystal neighborhood K0 into an imaginary intermediate configuration K volume. P is called the viscoplastic deformation with respect to the reference configuration of a material neighborhood of
Fig. 1. Various configurations associated with the lattice rotation of a viscoplastic crystal.
ARTICLE IN PRESS 846
O. Cazacu, I.R. Ionescu / J. Mech. Phys. Solids 58 (2010) 844–859
the material point X at time t. Crystal slip systems are labeled by integers s= 1,y,N, with N denoting the number of slip 0 0 systems. Each slip system s is specified by the unit vectors ðbs ,m0s Þ, where bs is in the slip direction and m0s is normal to the slip plane in the perfect undeformed lattice. Since the viscoplastic deformation does not produce distortion or rotation of the lattice, the mean lattice orientation is the same in the reference and intermediate configurations and is specified by 0 ðbs ,m0s Þ, s= 1,y,N. Let R denote the rotation of the crystal axes with respect to its isoclinic orientation. It maps the intermediate configuration into the current configuration. This leads to the following decomposition for the deformation gradient F (see for example Kok et al., 2002): F ¼ RP. This decomposition implies that we neglect the elastic lattice strains. Such a hypothesis is valid since during forming, the elastic component of deformation is negligibly small (typically 10 3) in comparison to the plastic component (typically 4101 , see also the channel die extrusion results of Fig. 9 where plastic deformations are up to 2). It also to be noted that once the elastic/plastic transition is over the stress evolution in the grains is controlled by plastic relaxation (see Tome and Lebensohn, 2003). Let bs ¼ bs ðtÞ and ms ¼ ms ðtÞ, the glide direction and glide plane normal, respectively, in the deformed configuration, i.e. 0 bs ð0Þ ¼ bs and ms ðtÞ ¼ m0s . Since elastic effects are neglected 0
ms ¼ Rm0s ,
bs ¼ Rbs ,
ð1Þ
Note that according to (1), bs and ms are unit vectors. Furthermore 0
bs ms ¼ Rðbs m0s ÞRT :
ð2Þ
We seek to express the lattice evolution equations only in terms of vector and tensor fields associated with the current configuration. Let v ¼ vðt,xÞ, the Eulerian velocity field, L the velocity gradient, D the rate of deformation, and W the spin tensor L ¼ LðvÞ ¼ rv,
D ¼ DðvÞ ¼ ðrvÞsymm ,
W ¼ WðvÞ ¼ ðrvÞskew :
ð3Þ
The viscoplastic deformation is due to slip only; the slip contribution to the viscoplastic deformation being (Rice, 1971; Teodosiu and Sidoroff, 1976) X s 0 g_ bs m0s , ð4Þ P_ P 1 ¼ s s
s
where g_ ¼ g_ ðtÞ is the viscoplastic shear rate on the slip system s. If we denote by M s ¼ ðbs ms Þsymm ,
Rs ¼ ðbs ms Þskew ,
ð5Þ
_ T þ RP_ P 1 RT , Eqs. (2) and (4), the rate of deformation D can be written as then, using L ¼ F_ F 1 ¼ RR D¼
N X
g_ s M s :
ð6Þ
s¼1
_ T A, where the skew tensor A is defined as Taking the antisymmetric part of L, we obtain that the spin tensor is W ¼ RR X s A ¼ g_ Rs : ð7Þ s
Next, by taking the material time derivative of (1) we obtain the Eulerian evolution equation for the glide direction vectors bs and for the normal vectors ms as d bs ¼ ðW þ AÞbs , dt
d ms ¼ ðW þ AÞms , dt
ð8Þ
where the material time derivative d=dt is expressed in terms of the velocity field by d @ ¼ þ v r: dt @t The evolution equations (8) describe the evolution of the lattice in terms of vector and tensor fields associated with the current configuration.
3. Overstress visco-plastic flow rule In order to complete the model, we need to provide the constitutive equation for the slip rate g_s as a function of ts , the stress component acting on the slip plane of normal ms in the slip direction bs . In the current configuration, ts is expressed as
ts ¼ r : M s ,
ð9Þ
where r ¼ rðtÞ is the Cauchy stress tensor acting in the current configuration K while M s is defined by (5). Note that if N 45 then fts gs ¼ 1,N are not independent; they belong to a fifth dimensional space of RN .
ARTICLE IN PRESS O. Cazacu, I.R. Ionescu / J. Mech. Phys. Solids 58 (2010) 844–859
847
In rigid-plastic formulations, it is assumed that the onset of plastic flow of a slip system s is governed by Schmid law: the slip system s is active if and only if jts j ¼ ts0 , i.e.
g_s ðjts jts0 Þ ¼ 0, g_s ts Z0, jts jts0 r 0,
ð10Þ
s 0
s 0
where t is the slip resistance (also called critical resolved shear stress or CRSS). For a given time t, the t are material constants. Thus, the planes jts j ¼ ts0 are the facets of the current yield surface of the single crystal in the stress space. For example, for B.C.C. crystals, there are 24 potentially active slip systems (N= 24). Since it is not known in advance which slip system will be active, all N slip rates must be treated as unknowns. It is clear that in order to accommodate arbitrary prescribed strains at least five slip systems need to be simultaneously activated. On the other hand, the solution is not unique. Additional assumptions are needed in order to restrict the number of solutions. One way to overcome this difficulty of determining the active slip systems problem is to adopt a rate-dependent approach for the constitutive response of the single crystal. A widely used rate-dependent (viscoplastic) model is the Norton type model, which relates the shear strain rate g_ s on a slip system s to the resolved shear stress ts through a power-law (see Asaro and Needleman, 1985) s n t g_ s ¼ g_ s0 s signðts Þ, ð11Þ
t0
where g_ 0 is a reference shear strain rate, while the exponent n has a fixed value. Generally, the exponent n has very large values (typically n = 20) so, as noted in the literature (e.g. McGinty and McDowell, 2006), this formulation predicts unrealistic slip rates for jts =ts0 j much larger than unity. Moreover, even for jts =ts0 j close to the unity the associated tangent problem has a huge variation of the slope for a small variation of jts =ts0 j (see Fig. 2). There have been a lot of efforts in order to overcome the severe numerical instabilities that arise in the integration of the model (11). A review of computational strategies can be found, for example, in Delannay et al. (2006) and Kuchnicki et al. (2008). In this paper, we consider a Perzyna-like viscoplastic law of the form (see also Fig. 2)
g_ s ¼
1
Zs
½jts jts0 þ signðts Þ,
ð12Þ
where Zs is the viscosity, which may depend on the slip rate, and ½x þ ¼ ðx þ jxjÞ=2 denotes the positive part of any real number x. Furthermore, substituting (12) in (6), follows that the Perzyna-like visco-plastic law expressing the rate of deformation D as a function of the Cauchy stress r, is N X ts0 1 1 ðr : M s ÞM s : ð13Þ D¼ Z jr : M s j þ s¼1 s Note that the viscoplastic flow rule (12) is the visco-plastic extension of the rigid-plastic Schmid law using an overstress approach. The physical motivation for the dependence of the viscoplastic shear rate on the overstress ðts ts0 Þ was provided by Teodosiu and Sidoroff (1976) based on an analysis of the microdynamics of crystals defects. The main assumption in (12) is that the resolved shear stress ts can be represented as the sum of a viscous contribution tsV ¼ Zs g_ s and a contribution tsP ¼ ts0 sgnðg_ s Þ, related to plasticity effects. The above flow rule can be rewritten as ( s t ¼ Zs g_ s þ ts0 sgn ðg_ s Þ if g_ s a0, ð14Þ if g_ s ¼ 0: jts j r ts0 s Overstress law
s 0
.s Power law s 0 Schmid law
Fig. 2. Schematic representation of two visco-plastic regularizations of the Schmid law (red line): a power law (blue line) and overstress Perzyna-type viscoplastic law (black line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
ARTICLE IN PRESS 848
O. Cazacu, I.R. Ionescu / J. Mech. Phys. Solids 58 (2010) 844–859
Since the resolved shear stresses ts are not independent, the shear rates g_s given by the viscoplastic flow rule (12) are not independent; they have to satisfy the kinematic constraint (6). Given the rate of deformation D, the slip rates g_s can be determined by minimizing the internal power Jðg_ 1 , g_ 2 , . . . , g_ N Þ ¼
N X Z
s
s¼1
2
jg_ s j2 þ ts0 jg_ s j,
ð15Þ
under the constraint (6). Slip on each system produces hardening on all slip systems. This is taken into account by adopting for the slip resistances ts0 an evolution law, which in Eulerian description is of the form: N X @ts d s t ¼ 0 þ v rts0 ¼ hsr jg_ r j, @t dt 0 r¼1
ð16Þ
with the initial condition ts0 ð0Þ ¼ tsc , with tsc constant for g_ s ¼ 0. The matrix hsr, called hardening matrix, describes the deformation resistance on system s which is caused by slip on system r. As shown by Franciosi (1984, 1985), the matrix hsr is not constant. Expressions for the hardening matrix that account for a monotonic dependence on the cumulated shear on all systems, g, defined as N X d @g g ¼ þ v rg ¼ jg_ s j, @t dt s¼1
ð17Þ
are widely used in FE simulations of polycrystals (see Chang and Asaro, 1981; Peirce et al., 1982; Anand and Kothari, 1996, etc.). Hardening laws that use dislocation densities on all slip planes as internal variables have also been proposed (e.g. Teodosiu and Sidoroff, 1976; Tabourot, 1992; Teodosiu et al., 1993, etc.). All the hardening laws mentioned are rate-independent, so they can be identified using quasi-static single crystal experiments. On the other hand, there are limited high-strain rate data from which evolution laws for the viscous coefficients Zs can be determined. If Zs are considered constant, then according to the overstress law (12) ts has a linear dependence on g_ s . Note that if
Zs ¼ A
argsinhðjg_ s j=g_ s0 Þ jg_ s j
the overstress law (12) reduces to the Prandtl-Eyring type model, while for
Zs ¼
ts0 g_ s 1=n1 g_ s0 g_ s0
we recover the Norton model (11).
4. Boundary value problem and numerical strategy We begin by presenting the equations governing the motion in a domain D ¼ DðtÞ of an incompressible rigidviscoplastic crystal. In an Eulerian description of a crystal visco-plasticity theory, the unknowns are: the velocity v : ½0,T D-R3 , the crystal lattice orientation, i.e. the N slip systems characterized by the slip directions bs : ½0,T D-S 3 and the normals to the slip planes ms : ½0,T D-S 3 (here S 3 denotes the unit vectors of R3 ), with bs ms ¼ 0, for all . Let r ¼ r0 þ pI, where r0 is the stress deviator while p : ½0,T D-R is s =1,y,N, and the Cauchy stress r : ½0,T D-R33 S the pressure. The momentum balance in the Eulerian coordinates reads
rð@t v þv = vÞdivr0 þ rp ¼ rf
in D,
ð18Þ
where the mass density r 4 0 and the body forces f are supposed to be known. The incompressibility condition is div v ¼ 0
in D:
ð19Þ
These balance equations are complemented by the constitutive equation (13), which relates the stress tensor r and the rate of deformation tensor DðvÞ, the evolution equations for each slip system s given by (8), and the hardening laws (e.g. (16)). The boundary @D of the domain D is decomposed into two disjoint parts, G0 and G1 , such that the velocity is prescribed on G0 and traction is prescribed on G1 , at any time t: vðtÞ ¼ VðtÞ on G0 ,
rn ¼ FðtÞ on G1 ,
ð20Þ
where n stands for the outward unit normal on @D, V is the imposed velocity and F is the prescribed stress vector. We also consider another partition of @D into @in DðtÞ and @out DðtÞ corresponding to incoming ðv n o 0Þ and outcoming ðv n Z0Þ
ARTICLE IN PRESS O. Cazacu, I.R. Ionescu / J. Mech. Phys. Solids 58 (2010) 844–859
849
flux. The boundary conditions associated to the lattice evolution equations (8) are bs ðtÞ ¼ Bs ðtÞ,
ms ðtÞ ¼ Cs ðtÞ,
s ¼ 1,: :,N
on @in DðtÞ,
ð21Þ
i.e. a given lattice orientation of the crystal is prescribed on @in DðtÞ. To the field equations, we add the initial conditions 0
vð0Þ ¼ v0 ,
bs ð0Þ ¼ bs ,
ms ð0Þ ¼ m0s ,
s ¼ 1,: :,N
in D:
ð22Þ
which express that the slip systems are known at t= 0. Note that in this model, there is no need to prescribe initial conditions for the stress. This is very convenient since the initial stress field is generally not known or it cannot be easily measured. To integrate the governing equations, a mixed finite-element and finite-volume strategy, developed by Cazacu and Ionescu (2010), will be used. We will give here only the principles of the method and present simulation results in the next section. A detailed presentation of the numerical scheme is given in Cazacu and Ionescu (2010). Time implicit (backward) Euler scheme is used for the field equations, which gives a set of nonlinear equations for the velocities v and lattice orientation ðbs ,ms Þ. At each iteration in time, an iterative algorithm is used in order to solve these nonlinear equations. Specifically, the variational formulation (inequality) for the velocity field is discretized using the finite element method, while a finite volume method with an upwind choice of the flux is adopted for solving the hyperbolic equations that describe the evolution of the lattice orientation. One of the advantages of the proposed numerical strategy is that it does not require the consideration of elastic deformation since it is not based on an elastic predictor/plastic corrector scheme (e.g. Simo and Hughes, 1998). It is to be noted that in the case of the proposed rigid-viscoplastic model (13), numerical difficulties arise from the non-differentiability of the viscoplastic terms. That means that one cannot use finite element techniques developed for Navier–Stokes fluids (see for instance Pironneau, 1989; Temam, 1979). To overcome these difficulties the iterative decomposition-coordination formulation coupled with the augmented Lagrangian method of Glowinski and Le Tallec (1989) and Fortin and Glowinski (1982) was modified. The reason for this modification is that in the crystal model (13) there is no co-axiality between the stress deviator and the rate of deformation as it is the case in the Bingham model used by Glowinski and Le Tallec (1989) and Fortin and Glowinski (1982). Note that this iterative decomposition-coordination algorithm permits to solve at each iteration the equations for the unit vectors that define the lattice orientation. Furthermore, if the domain D occupied by the single crystal (or multicrystals) varies in time (for example, see the last application presented in Section 9), then an arbitrary Eulerian–Lagrangean (ALE ) description was adopted. In the numerical simulations presented in the last two sections we use the followings spatial discretization: P2 for the velocity field, P1 for the pressure, P1 discontinuous for the stress field, slip rates, and lattice orientations. 5. In-plane form of the model In this section, we consider the problem of plane strain deformation of a crystal obeying the rigid-viscoplastic Eulerian model described by Eqs. (8) and (13). We show that for a rigid-viscoplastic crystal, the same slip systems as in the rigidplastic case need to be simultaneously activated in order to accommodate plane strain deformation. Furthermore, we prove that the lattice rotation is described by one angle and its evolution by a single differential equation. This analytical treatment gives insight into the physics of the deformation and allows a deeper understanding and interpretation of the numerical results that will be presented in Sections 7 and 8. We begin by briefly recalling the main results of Rice (1973) (see also Crone et al., 2004; Gan et al., 2006) concerning the plane strain, rate-independent deformation of rigid-plastic single crystals obeying Schmid law. Rice (1971) showed that certain pairs of the three-dimensional systems that are potentially active need to combine in order to achieve plane-strain deformation. Let r =1,y,Np be the index of the composite slip system formed by two 3D slip systems k and l. For a pair of slip systems k and l to produce in-plane deformation in the plane (x1x2), their slip directions and normals must satisfy the following conditions (see Rice, 1973; Crone et al., 2004; Gan et al., 2006): bk3 7bl3 ¼ 0,
mk3 7 ml3 ¼ 0,
bk3 mk3 ¼ 0,
bl3 ml3 ¼ 0,
b k b l ¼ 7 1,
m k m l ¼ 71,
ð23Þ
where b k , b l , m k , m l denote the normalized projections of the three-dimensional slip directions bk , bl and normal directions mk , ml onto the x1x2-plane. Note that since (bk3)2 = (bl3)2 and (mk3)2 = (ml3)2 we can put qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi br :¼ 1ðbk3 Þ2 ¼ 1ðbl3 Þ2 , mr :¼ 1ðmk3 Þ2 ¼ 1ðml3 Þ2 , qr :¼ mr br , ð24Þ for r =1,y,Np, thus k bk ba ¼ a ,
br
l bk ba ¼ a ,
br
m ka ¼
mka
mr
,
m la ¼
mla
mr
,
a ¼ 1,2:
Note that the 2-D vectors b k , b l , m k , m l are also orthogonal in the plane (i.e. b k m k ¼ 0, b l m l ¼ 0). Since b k m k ¼ b l m l we can define the in-plane slip directions b r ¼ b k and plane normals m r ¼ m k (or equivalently b r ¼ b l ,m r ¼ m l ). We can introduce now M r , the symmetric part of the two-dimensional tensor product of slip directions
ARTICLE IN PRESS 850
O. Cazacu, I.R. Ionescu / J. Mech. Phys. Solids 58 (2010) 844–859
and normals M r ¼ ðb r m r Þsymm :
ð25Þ
Note that, for any r, the skew part of b r m r has the same expression, denoted by R 1 R r ¼ ðb r m r Þskew ¼ R ¼ ðð1,0Þ ð0,1Þð0,1Þ ð1,0ÞÞ: 2
ð26Þ
In the following, we will deduce the governing equations of the in-plane flow through a domain D ¼ O R, of a crystal described by the proposed rigid-viscoplastic model. These equations take a simpler form if we make use of only twodimensional vector and tensor quantities. For plane-strain conditions the unknowns of the problem are: the in-plane , where s ab ¼ sab , the in-plane velocity v : ½0,T O-R2 ðv ¼ ðv,0Þ and v ¼ vðt,x1 ,x2 Þ), Cauchy stress r : ½0,T O-R22 S and the angle y : ½0,T O-R, that specifies the current orientation of the lattice. Proof that these are the only unknowns of the problem is given in the following. Since the out of plane shear components of the stress tensor are zero (i.e. s13 ¼ s23 ¼ 0) from (23), it follows that r : M k ¼ r : M l . If we assume that the two 3D systems k and l have the same yield limit and the same viscosity, denoted by tr0 and by Zr , respectively
tl0 ¼ tk0 ¼ tr0 , Zl ¼ Zk ¼ Zr ,
ð27Þ
then from the Perzyna-type flow rule (12), we obtain that g_ k ¼ g_ l . This common value will be denoted in what follows by g_ r . Also, since bk3mki + bl3mli =0 for any i= 1,y,3, we obtain ðg_ k M k þ g_ l M l Þ3i ¼ 0,
ðg_ k Rk þ g_ l Rl Þ3i ¼ 0
for any i ¼ 1,2,3:
ð28Þ
This means that the combination of the slip systems k and l satisfying (23) will result in plane strain deformation. Thus, in order to accommodate plane strain deformation in the rigid-viscoplastic case, the same geometrical constraints (23) have to be satisfied as in the rigid-plastic case. Next, we derive the in-plane form of the rigid viscoplastic flow rule (12). Let r =1,y, Np be the index of the composite slip system formed by the three-dimensional slip systems k and l. r Denoting by g_ ¼ 2qr g_ r , using Eqs. (25), we get r
r
ðg_ k M k þ g_ l M l Þab ¼ g_ ðM r Þab ,
ðg_ k Rk þ g_ l Rl Þab ¼ g_ ðRÞab
for all a, b ¼ 1,2:
ð29Þ r
Furthermore, since r : M l ¼ qr r : M l , the overstress based flow rule (12) becomes a relation between g_ and r , the two dimensional tensor obtained from the first two columns and lines of the stress tensor r r 1 t r0 g_ ¼ 1 r : Mr: ð30Þ Zr jr : M r j þ In (30), Z r0 and t r0 denote the ‘‘in-plane’’ viscosity and yield limit, respectively
t r0 ¼
tr0 qr
,
Zr
Zr ¼
2q2r
:
ð31Þ
Note that in order to have a plane strain configuration, all three-dimensional slip systems s that do not form twodimensional composite slip systems should not be active (i.e. for s different than any (k,l) satisfying (23), g_ s ¼ 0). This is the case, since from (28) and (29) it follows that the only non-zero components of the rate of deformation D are Dab , with a, b ¼ 1,2, i.e. Di3 ¼ D3i ¼
N X
s g_ s M3i ¼ 0, Dab ¼
s¼1
N X
g_ s Mas b ¼ 2
s¼1
Np X r¼1
g_ r Mar b ¼
Np X
r
r
g_ M ab ,
ð32Þ
r¼1
for all i= 1,2,3. Let also denote by DðvÞ :¼ ð= vÞsymm and by WðvÞ :¼ ð= vÞskew , the symmetric and antisymmetric parts, respectively, of the two-dimensional velocity gradient = v. Note that Eq. (32) expresses that for plane strain deformation in the Ox1x2 plane, i.e. for velocity v ¼ ðv,0Þ, with vðt,x1 ,x2 Þ ¼ ðv1 ðt,x1 ,x2 Þ; v2 ðt,x1 ,x2 ÞÞ 2 R2 , the two-dimensional rate of deformation tensor DðvÞ can be decomposed into Np composite plane-strain systems specified by ðb r ,m r Þ, the normalized projections of the three-dimensional slip directions and normal directions onto the x1x2-plane. Thus DðvÞ ¼
Np X
r
g_ M
r
ð33Þ
r¼1
Furthermore, the proposed rigid viscoplastic law can be written now as a relation between the two-dimensional stress tensor r and the two-dimensional rate of deformation tensor DðvÞ DðvÞ ¼
Np X 1 r¼1
Zr
1
t r0
jr : M r j
ðr : M r ÞM r : þ
ð34Þ
ARTICLE IN PRESS O. Cazacu, I.R. Ionescu / J. Mech. Phys. Solids 58 (2010) 844–859
851
As already mentioned, for plane strain deformation, important simplifications of the evolution Eqs. (8) that characterize the change in lattice orientations can be deduced. Indeed, from (29) it follows that ! ! Np Np N N X X X X r s s s s r _ _ _ _ g R ¼ 0, g R ¼2 g R ¼ g R , ð35Þ ab
ab
3i
s¼1
r¼1
s¼1
ab
r¼1
for all i= 1,2,3 and a, b ¼ 1,2. From the above equation we get ðAÞ3i ¼ 0 (see Eq. (7) for the definition of the tensor A). Since ðWðvÞÞ3i ¼ 0, from the evolution Eqs. (8) we obtain d s b ¼ 0, dt 3
d s m ¼0 dt 3
for all s ¼ 1, . . . ,N:
In other words, all active three-dimensional slip systems are characterized by slip and normal vectors having the component along the normal to the plane of deformation constant. Moreover for any composite slip system r, the factors br and mr do not depend on time. Thus, the change in lattice orientation can be described by the evolution equations for the projections of the 3-D slip systems vectors onto the Ox1 x2 plane. Again using (8) we obtain ! ! ! ! Np Np X X d r d r b r ¼ WðvÞ g_ R b r , m r ¼ WðvÞ g_ R m r : ð36Þ dt dt r¼1 r¼1 It is convenient to specify the two-dimensional vectors b r and m r by their polar coordinates b r ¼ ðcosyr ,sinyr Þ,
m r ¼ ðsinyr ,cosyr Þ:
ð37Þ
The in-plane traction associated with the composite slip system r is 1 2
t r ¼ r : M r ¼ sinð2yr Þðs22 s11 Þ þcosð2yr Þs12 , while the Eq. (33), expressing the decomposition of the in-plane rate of deformation into Np composite in-plane systems, becomes Np Np r r @v1 1X 1 @v2 @v1 1X ¼ ¼ g_ sinð2yr Þ, D12 ¼ þ g_ cosð2yr Þ: ð38Þ D11 ¼ @x1 @x2 2r ¼1 2 @x1 2r ¼1 r Furthermore, given DðvÞ, the shear rates g_ are determined by minimizing the internal power Np 1 2 X Np Zr _ r 2 r _ r jg j þ t 0 jg j, Jplane g_ , g_ , . . . , g_ ¼ 2 r¼1
ð39Þ
under the constraint (38). Note that the Eqs. (36) can be written as d 1 XNp _ r @v1 @v2 , yr ¼ g r¼1 @x2 @x1 dt 2 for any r = 1,y,Np. It follows that d=dtðyr yq Þ ¼ 0, with raq. Since the angles between any two systems do not change in time, it is sufficient to compute the change in orientation of only one of the composite in plane slip systems, say y ¼ y1 , by using the equation ! Np r d @y 1 X @v1 @v2 _ : ð40Þ þ v ry ¼ y¼ g @t @x2 @x1 dt 2 r¼1 The lattice orientation of all the other Np 1 composite plane strain systems is obtained from the relation
yr ðtÞ ¼ y1 ðtÞ þ yr ð0Þy1 ð0Þ:
ð41Þ
For plane strain conditions in the Ox1x2 plane, the three-dimensional momentum balance law (18) and the incompressibility condition in D ¼ O R read
rð@t v þ v = vÞdivr 0 þ rp ¼ rf , div ðvÞ ¼ 0, in O:
ð42Þ 2
The boundary value problem in the in-plane configuration consists in finding the velocity field v : ½0,T O-R , the , the pressure p : ½0,T O-R and the angle y : ½0,T O-R deviator of the Cauchy stress tensor r 0 : ½0,T O-R22 S solution of Eqs. (34), (40) and (42), with the tensors M r defined by (25) and (37). 6. In-plane configurations of F.C.C. crystals We will further apply the in-plane model to a F.C.C. crystal. Let Ox3 axis be parallel to [1 1 0] in the crystal basis, which means that the plane-strain plane (Ox1x2) is the plane ½1 1 02½0 0 1. For a F.C.C. crystal there are 12 potentially active slip systems. The geometrical constraints (23) are satisfied only by three pairs of these systems (Np = 3). As shown in Fig. 3 the
ARTICLE IN PRESS 852
O. Cazacu, I.R. Ionescu / J. Mech. Phys. Solids 58 (2010) 844–859
c3
~
x2 2
2
~ 1
c3
x1
~
1
~ 3
x2 2
c2
5 ~ 3 46
~ 1
3
x1
c2 _ m = (1 1 1)
m = (1 1 1)
c1
c1
x3
c3
x3
c3
x1
x1
12 ~ 2
7
~ 3
x2
~ 1
8
9
c1
c2
~ 2 x2 ~ 1
_ m = (1 1 1)
11 ~ 3
c2 __ m = (1 1 1)
c1
x3
10
x3
Fig. 3. Geometrical constraints (23) to be satisfied by a F.C.C. crystal deforming in the plane ½110½001. These three (Np = 3) composite in-plane systems are: 1~ ¼ ð3,12Þ, 2~ ¼ ð4,5Þ and 3~ ¼ ð7,8Þ.
in-plane system r ¼ 1~ is formed by (k,l) =(3,12), r ¼ 2~ is formed by (k,l)= (4,5), and r ¼ 3~ is formed by (k,l)= (7,8). We suppose in what follows that the systems 1,2,6,9,10 and 11 are not active. To describe the orientation of the crystal, we denote by y the angle counterclockwise between the Ox1 axis and the ½110 direction. The three active composite in-plane slip systems b 1 , b 2 , b 3 are specified by the angles pffiffiffi ð43Þ y1 ¼ y, y2 ¼ yf, y3 ¼ y þ f with f ¼ arctanð 2Þ, while the corresponding in-plane factors (see the definition (24) for qr) are pffiffiffi 1 3 : q1 ¼ pffiffiffi , q2 ¼ q3 ¼ 2 3
ð44Þ
r The in-plane slip rates g_ (r =1,2,3) are related to the 3-D slip rates g_ s (s =1,y,12) by 1
2
2
pffiffiffi
2
pffiffiffi
3
pffiffiffi
pffiffiffi
g_ ¼ pffiffiffi g_ 3 ¼ pffiffiffi g_ 12 , g_ ¼ 3g_ 4 ¼ 3g_ 5 , g_ ¼ 3g_ 7 ¼ 3g_ 8 , 3
3
while the apparent in-plane viscosities and the apparent in-plane critical resolved shear stresses are 3 2
3 2
2 3
2 3
2 3
2 3
Z 1 ¼ Z3 ¼ Z12 , Z 2 ¼ Z4 ¼ Z5 , Z 3 ¼ Z7 ¼ Z8 , pffiffiffi
pffiffiffi
2
2
2
2
t 10 ¼ 3t30 ¼ 3t12 t 20 ¼ pffiffiffi t40 ¼ pffiffiffi t50 , t 30 ¼ pffiffiffi t70 ¼ pffiffiffi t80 : 0 , 3
3
3
3
We derive here, analytic expressions for the slip rates on individual composite systems. Specifically, for a given lattice orientation, we want to find the decomposition of the in-plane rate of deformation into the three in-plane composite slip r systems. For the 2-D case the slip rates g_ (r = 1,2,3) can be determined by minimizing the internal power (39) under the constraints (38). It is worth noting that these slip rates can be determined in closed form. The proof is given in the following. Let us first decompose the three dimensional vector:
G_ ¼ ðg_ 1 , g_ 2 , g_ 3 Þ ¼ sS y þcC y þ zZ y
y
ð45Þ y
in the basis fS ,C ,Zg given by S ¼ ðsinð2y1 Þ,sinð2y2 Þ,sinð2y3 ÞÞ, C ¼ ðcosð2y1 Þ,cosð2y2 Þ,cosð2y3 ÞÞ and Z ¼ S 4C y ¼ ðsinð4fÞ,sinð2fÞ,sinð2fÞÞ. Substitution of (45) into Eq. (38) leads to sjS y j2 þcS y C y ¼ 2D11 ,
y
sS y C y þcjC y j2 ¼ 2D12 ,
y
ð46Þ
where jS y j2 ¼
3 X
sin2 ð2yr Þ,
Sy C y ¼
r¼1
3 1X sinð4yr Þ, 2r ¼1
y 2
and jC j ¼ 3jS y j2 . Solving the above linear system (46) for the unknowns s and c, we obtain s ¼ 2
jC y j2 D11 þ ðS y C y ÞD12 y 2
y 2
y
y 2
jS j jC j ðS C Þ
,
c¼2
jS y j2 D12 þðS y C y ÞD11 jS y j2 jC y j2 ðS y C y Þ2
:
ð47Þ
ARTICLE IN PRESS O. Cazacu, I.R. Ionescu / J. Mech. Phys. Solids 58 (2010) 844–859
853
1 2 3 Thus, the in-plane internal power is a function of z only, i.e. Jplane g_ , g_ , g_ ¼: JðzÞ, with 3 X Z
JðzÞ ¼
r¼1
r
2
jsSyr þ cC yr þ zZ r j2 þ t r0 jsSyr þ cC yr þzZ r j:
ð48Þ
The minimization problem to be solved reduces to J0 ðzÞ ¼ 0. Note that J(z) is differentiable everywhere with the exception of three points, zy1 , zy2 and zy3 , given by zy1 ¼
ssinð2y1 Þ þ ccosð2y1 Þ , sinð4fÞ
zy2 ¼
ssinð2y2 Þ þ ccosð2y2 Þ , sinð2fÞ
zy3 ¼
ssinð2y3 Þ þ ccosð2y3 Þ : sinð2fÞ
However, the left and right derivatives J 0 ðzyr Þ,J0 ðzyr þÞ exist and can be obtained from the general expression of J 0 ðzÞ as: J0 ðzÞ ¼ Ay þ Bzþ
3 X
t r0 jZr j
r¼1
zzyr jzzyr j
with Ay ¼
3 X
Z r ðsSyr þ cC yr ÞZr , B ¼
r¼1
3 X
Z r Zr2 :
r¼1
We can now solve the equation J 0 ðzÞ ¼ 0 to get the analytical expression of z 8 y if J 0 ðzyr ÞJ 0 ðzyr þÞ r0, z > > < r" # 3 X J 0 ðzy þÞ z¼ 1 > otherwise t r0 jZr j 0 ry Ay þ > :B jJ ðzr þÞj r¼1
ð49Þ
Further substitution of (47) and (49) into (45) gives the analytical decomposition of the strain rate DðvÞ into the three slip rates g_ 1 , g_ 2 , g_ 3 according to the constitutive equation (34). 7. Extrusion in an equal channel angular die In this section, we apply the in-plane model, described by Eqs. (34) and (40), to the analysis of the equal channel angular die extrusion of an F.C.C single crystal. The equal channel angular extrusion (ECAE) was first developed by Segal (1981) and intensely studied thereafter (see for instance Luri et al., 2006). The main advantage of using ECAE processes rather than traditional manufacturing processes, such as rolling, extrusion, or conventional forging is the slight variation on the cross section of the processed billet. This allows to process the material several times in order to achieve the desired accumulated deformation. The material is extruded through a die formed by two channels that intersect at a 903 angle (see Fig. 4 left for a schematic representation of the die geometry ). As the billet crosses the intersection of these channels, the material is severely deformed in shear to a very high level of accumulated plastic strains. Submicrometric or even nanometric grain size can be achieved by using ECAE. This grain size reduction improves mechanical properties such as yield stress and ultimate stress. The deformation field in the ECAE process can deviate significantly from idealized shear deformation. Even with models that do not account for hardening or consider only isotropic hardening (Rosochowski and Olejnik, 2002; Li et al., 2006), the predicted deformation history in the ECAE processed sample is highly heterogeneous, influenced strongly by the geometry of the ECAE dies. However, the predicted textures obtained using a Taylor-type crystal plasticity model with hardening are in better agreement with the experiments (Toth et al., 2004; Li et al., 2005; Kalidindi et al., 2008). The main problem we focus here is: given the initial orientation of an FCC single crystal, determine the main orientation of the lattice after deformation in an equal channel extruder i.e. its orientation yf at the exit of the extruder (see Fig. 4 left), r determine the slip systems which are active and the distribution of the slip rates g_ (r =1,2,3) in the extruder. The die
b2
IsoValue 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
b1
θ0
x2 V
[110]
V
x3
b2 [001]
[110] b3 b1 θf
x1
Fig. 4. Left: extruder geometry and the orientation of the crystal given by the angles y0 (at the entrance) and yf (at the exit). Right: a typical distribution of the modulus of the rate of deformation (stationary) in the extruder, i.e. ðx1 ,x2 Þ-jDðv Þjðx1 ,x2 Þ.
ARTICLE IN PRESS 854
O. Cazacu, I.R. Ionescu / J. Mech. Phys. Solids 58 (2010) 844–859
geometry and the coordinate system for the analysis are depicted in Fig. 4 left, with Ox3 axis chosen to be parallel to [1 1 0]. Since the flow takes place in the plane ½1 1 02½0 0 1, the mechanical response is governed by the constitutive relations (34) and (40) and the balance equations (42). As already shown in the previous section, there are only three composite in-plane systems (Np = 3). At any point (x1,x2) of the extruder, the orientation of the crystal (i.e. the three active composite in-plane slip systems specified by (43)) is described by the angle y counterclockwise between ½1 1 0 and the Ox1 axis of the extruder; the corresponding in-plane factors qr, r= 1yNp, are given by (44). In all the simulations, the entrance velocity was v ¼ ð0,VÞ with V =10 m/s while the numerical values for the material coefficients were: r ¼ 3,000 kg=m3 , same slip resistance for all 3-D slip systems ts0 ¼ 10 MPa (thus, t 10 ¼ 17:32 MPa, t 20 ¼ t 30 ¼ 11:55 MPa, while the viscosity was considered constant Zs ¼ 2 kPa s (hence Z 1 ¼ 3 kPa s, Z 2 ¼ Z 3 ¼ 1:33 kPa s). Hardening was neglected. The characteristic length of the process was considered to be L= 10 2 m, hence the Reynolds and the Bingham numbers, are Re ¼
rVL t0 L ¼ 5: ¼ 0:15, Bi ¼ Z ZV
ð50Þ
Let us first analyze the steady-state behavior of a single crystal. It is worth noting that in an Eulerian setting, it is sufficient to do the computations only in a region around the extruder angle, whereas in a Lagrangian description the computations need to be done in all the material points of the sample. We have conducted simulations for two different values of the crystal orientation at the entrance of the extruder ye ¼ 0 and ye ¼ p=2, respectively. In Fig. 4 (right) is presented the distribution of the rate of deformation for ye ¼ p=2. Note that deformation occurs only in a region around the 903 angle of the extruder while at the entrance and exit of the extruder the crystal is rigid. Similar distributions of the modulus of the rate of deformation tensor are obtained irrespective of the initial orientation of the crystal. According to our model, when the crystal is subjected to in-plane deformation the lattice orientation can be described by a unique scalar evolution Eq. (40). Fig. 5 shows the stationary distribution ystat of the lattice orientation for the two initial (entrance) orientations of the crystal. Note that the distributions of ystat at the exit of the extruder is not exactly homogeneous (the inhomogeneity being more pronounced for ye ¼ p=2) thus, we cannot have a single value for the exit orientation yf . However, since y is almost constant on a large zone, the mean-value yf of ystat ðx1 ,x2 Þ at the exit of the extruder gives an estimate of the crystal rotation. If the entrance orientation of the crystal is ye ¼ 0 then yf ¼ 0:9123, which means that the crystal rotated with a1 ¼ 0:9123 ¼ 52:263 . In the second case, we obtained yf ¼ 1:13108 which corresponds to a rotation of a2 ¼ 0:43977 ¼ 25:23 . This effect of the initial orientation on the rotation of the lattice can be r explained by analyzing the slip activity in each case. Fig. 6 shows the distribution of the slip rates g_ for each of the three potentially active composite systems for the case when the initial orientation of the crystal ye ¼ 0. Note that for ye ¼ 0 the active composite slip systems are r = 1 and 3, while for an entrance orientation ye ¼ p=2 the active composite slip systems are r = 2 and 3. Next, we analyze the evolution of the lattice orientation of a square grain embedded in another grain. Initially, the ‘‘parent grain’’ has the orientation ystat ðx1 ,x2 Þ, the stationary distribution corresponding to an entrance orientation ye , (see the previous stationary simulation), while the ‘‘inserted grain’’, located at the entrance of the extruder, has an initial orientation ygrain aye . In other words, yjt ¼ 0 ¼ ystat ðx1 ,x2 Þ for all points (x1,x2) in the extruder (see Fig. 5 for the spatial distribution of ystat ), with the exception of a zone occupied by the square grain, where yjt ¼ 0 ¼ ygrain . We have considered two cases. In the first case, we set ye ¼ 0 and ygrain ¼ p=2, while in the second case: ye ¼ p=2 and ygrain ¼ 0. The evolution in time of the lattice orientations during the extrusion process is plotted in Fig. 7 for both cases. Intuitively, we expect the ‘‘inserted grain’’ to accommodate its orientation with the ‘‘parent crystal’’ as soon as it enters the deformation zone of the IsoValue -1.9 -1.65 -1.4 -1.15 -0.9 -0.65 -0.4 -0.15 0.1 0.35 0.6 0.85 1.1 1.35 1.6 1.85
IsoValue -1.9 -1.65 -1.4 -1.15 -0.9 -0.65 -0.4 -0.15 0.1 0.35 0.6 0.85 1.1 1.35 1.6 1.85
Fig. 5. Stationary distribution of the crystal orientation ðx1 ,x2 Þ-ystat ðx1 ,x2 Þ for two entrance orientations: ye ¼ 0 (left) and ye ¼ p=2 (right).
ARTICLE IN PRESS O. Cazacu, I.R. Ionescu / J. Mech. Phys. Solids 58 (2010) 844–859 IsoValue -0.45 -0.35 -0.25 -0.15 -0.05 0.05 0.15 0.25 0.35 0.45
855 IsoValue -0.45 -0.35 -0.25 -0.15 -0.05 0.05 0.15 0.25 0.35 0.45
IsoValue -0.45 -0.35 -0.25 -0.15 -0.05 0.05 0.15 0.25 0.35 0.45
r Fig. 6. The stationary distribution ðx1 ,x2 Þ-g_ stat ðx1 ,x2 Þ (left figure corresponds to r= 1, middle figure corresponds to r =2, while r =3 corresponds to the right figure) for the entrance orientation ye ¼ 0. Blue corresponds to zero slip rate. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
IsoValue -1.9 -1.65 -1.4 -1.15 -0.9 -0.65 -0.4 -0.15 0.1 0.35 0.6 0.85 1.1 1.35 1.6 1.85
IsoValue -1.9 -1.65 -1.4 -1.15 -0.9 -0.65 -0.4 -0.15 0.1 0.35 0.6 0.85 1.1 1.35 1.6 1.85
IsoValue -1.9 -1.65 -1.4 -1.15 -0.9 -0.65 -0.4 -0.15 0.1 0.35 0.6 0.85 1.1 1.35 1.6 1.85
IsoValue -1.9 -1.65 -1.4 -1.15 -0.9 -0.65 -0.4 -0.15 0.1 0.35 0.6 0.85 1.1 1.35 1.6 1.85
IsoValue -1.9 -1.65 -1.4 -1.15 -0.9 -0.65 -0.4 -0.15 0.1 0.35 0.6 0.85 1.1 1.35 1.6 1.85
IsoValue -1.9 -1.65 -1.4 -1.15 -0.9 -0.65 -0.4 -0.15 0.1 0.35 0.6 0.85 1.1 1.35 1.6 1.85
IsoValue -1.9 -1.65 -1.4 -1.15 -0.9 -0.65 -0.4 -0.15 0.1 0.35 0.6 0.85 1.1 1.35 1.6 1.85
IsoValue -1.9 -1.65 -1.4 -1.15 -0.9 -0.65 -0.4 -0.15 0.1 0.35 0.6 0.85 1.1 1.35 1.6 1.85
Fig. 7. The distribution ðx1 ,x2 Þ-yðt,x1 ,x2 Þ of the crystal lattice orientation during the extrusion process (t = 0, t = T/3, t= 2T/3, t= T from the left to the right) for the cases: ye ¼ 0, ygrain ¼ p=2 (top) and ye ¼ p=2, ygrain ¼ 0 (bottom).
extruder (i.e. where jDðvÞjðx1 ,x2 Þa0). However, analysis of Fig. 7 shows that at the exit of the extruder, the inserted square grain has orientation a ¼ a1 ¼ 0:9123 ¼ 52:263 , if ygrain ¼ 0 while a ¼ a2 ¼ 0:43977 ¼ 25:23 if ygrain ¼ p=2, i.e. in each case the inserted grain has the same exit orientation as in the stationary case. This means that the ‘‘parent crystal’s’’ initial orientation and rotation do not influence the lattice rotation of the ‘‘inserted grain’’. Yet, the ‘‘parent grain’’ influences the deformation of the inserted grain (note that the shape of inserted grain at the exit of the extruder is different in the two cases). While in both cases, the sides of the grains that initially were vertical have a p=2 rotation (a little less in the second case), the horizontal sides have a 241 rotation in the first case, and a 271 rotation in the second case (see Fig. 7).
8. Channel-die compression In this section, we analyze the flow of an F.C.C. crystal compressed in a channel die. In channel-die compressions tests (see for instance Maurice and Driver, 1993 for the device description) if friction is small, almost perfect plane strain conditions can be achieved. Assume that [1 1 0] axis of the crystal is along Ox3 axis of the die (see Fig. 8 for a schematic representation of the die geometry). The crystal flow is in the plane ½1 1 02½0 0 1. According to the proposed in-plane
ARTICLE IN PRESS 856
O. Cazacu, I.R. Ionescu / J. Mech. Phys. Solids 58 (2010) 844–859
[110]
x2
x2
V V
_ b2 _ b1
θ _ b3
[110] x3
x3
V
x1
x1 Fig. 8. The geometry of the channel-die compression flow.
model (see Eqs. (34) and (40)) there are only three active composite in-plane systems (Np = 3) and the orientation of the crystal at any point (x1,x2) and at any time t is described through the angle yðt,x1 ,x2 Þ between ½1 1 0 and the Ox1 axis of the channel-die. The three active composite in-plane slip systems are specified by (43), while the corresponding in-plane factors are given by (44). Let first analyze the influence of the viscosity and of the inertial forces. The material coefficients considered correspond to copper (Cu) of density r ¼ 8:940 kg=m3 , slip resistance for all systems equal to ts0 ¼ 60:8 MPa, which means that for each in-plane system t 10 ¼ 105:84 MPa, t 20 ¼ t 30 ¼ 69:92 MPa; no hardening is taken into consideration. In its initial configuration the single crystal sample has a rectangular section of Lc 0:5Lc , with Lc = 0.1 m denoting the characteristic length, and an initial orientation defined by the angle yð0,x1 ,x2 Þ ¼ 23:93 ¼ 0:417 rad between ½110 and the Ox1 axis. The time interval [0,T] was chosen such that the final height corresponds to half of the initial height of the sample, i.e. a 50% reduction. To examine the influence of the viscosity and of the inertial forces three simulations were conducted for two different crosshead velocities V and for two different viscosities Zs . The reference crosshead velocity on the top of the sample v ¼ ð0,VÞ was chosen to be V0 =1 m/s, while the reference viscosity was Z0 ¼ 20 kPa s (i.e. Z 01 ¼ 1:5Z0 , Z 20 ¼ Z 03 ¼ 0:66Z0 ). The associated reference non-dimensional numbers, given by (50) are Re = 0.0447 and Bi= 304. In the first simulation Zs ¼ Z0 ,V ¼ V0 , in the second one Zs ¼ 250Z0 ,V ¼ V0 (Re = 0.0001788, Bi = 1.216) and in the third one, we set Zs ¼ Z0 ,V ¼ 250V0 (Re = 11.175, Bi= 1.216). In Fig. 9 (bottom) are plotted for all these three cases, the distributions of the cumulated shear deformation ðx1 ,x2 Þ-gðt,x1 ,x2 Þ (given by (17)), while on the top of Fig. 9 are plotted the corresponding lattice rotations ðx1 ,x2 Þ-yðt,x1 ,x2 Þ given by the angle between ½1 1 0 and the Ox1 axis. To analyze the influence of the viscous effects we fixed the crosshead velocity and we increased the viscosity by a factor of 250. Note that even at relatively small strain rates (10 s 1), the viscosity plays an important role in strain localization. For smaller viscosity (Fig. 9 left), we remark the occurrence of shear bands parallel to the slipping direction b 3 . At the beginning of the deformation process there are two shear bands. However, to accommodate the large deformation (50% in height reduction) at the end of the process, a third shear band appears in the middle of the sample. Note that the lattice has an important rotation only in these shear bands. For larger viscosity, localization takes place but is much more diffuse (see Fig. 9 middle). Only the upper left corner and the bottom right corner of the sample have an important lattice rotation. These two zones were initially at the lateral free surface and during the deformation they come into contact with the crosshead and with the bottom plate. Even if the lattice rotations are quite different, the shapes of the sample at the final stage of deformation are similar in both cases (compare Fig. 9 left and middle). The effect of the inertial forces can be observed by comparing the results of the second computation (see Fig. 9 middle) and the third one (see Fig. 9 right). Indeed, the viscous contributions tsV ¼ Zs g_ s and the plastic stress contributions tsP ¼ ts0 sgn ðg_ s Þ are the same (see Eq. (14)) hence, the Cauchy stress tensor is of the same order. From the momentum balance equation, it follows that the difference between these two cases consists only in the inertial effects. We remark that in the third case, the shear bands have completely disappeared, the lattice has an almost uniform small rotation, and the final shape of the sample is completely different from the other cases. To analyze the combined effects of viscosity and inertia, we kept the same material parameters and we increased the crosshead velocity by a factor of 250. Comparison of the results of the first computation (see Fig. 9 left) and third one (see Fig. 9 right) shows that the final sample shapes and the final lattice orientations are very different. At high crosshead velocity, the plastic properties and crystal anisotropy are less important while inertia effects are dominant. Next, we analyze the interaction between grains in a multicrystal subjected to channel die compression. The model is applied to a tantalum (Ta) multicrystal of 4 grains. The material coefficients are: density r ¼ 16,650 kg=m3 , slip resistance for all systems equal to ts0 ¼ 66 MPa (i.e. t 10 ¼ 114:31 MPa, t 20 ¼ t 30 ¼ 76:23 MPa) with no hardening. The viscosity was Zs ¼ 44 kPa s (i.e. Z 10 ¼ 66 kPa s, Z 20 ¼ Z 30 ¼ 29:33 kPa s). The crosshead velocity was set to V =1 m/s; the associated non-dimensional numbers, given by (50) are Re= 0.0378 and Bi =150.
ARTICLE IN PRESS O. Cazacu, I.R. Ionescu / J. Mech. Phys. Solids 58 (2010) 844–859
857
Fig. 9. The distribution of the lattice rotation ðx1 ,x2 Þ-yðt,x1 ,x2 Þ (top) and of the cumulated shear deformation on all systems ðx1 ,x2 Þ-gðt,x1 ,x2 Þ (bottom) at different stages of the deformation (t = 0, t = T/3, t = 2T/3, t =T) for two different crosshead velocities V and for two different viscosities Zs . Left: the reference setting Zs ¼ Z0 ,V ¼ V0 , middle: Zs ¼ 250Z0 ,V ¼ V0 , right: Zs ¼ Z0 ,V ¼ 250V0 .
Fig. 10. The distribution of the lattice rotation ðx1 ,x2 Þ-yðt,x1 ,x2 Þ (left) and of the strain rate ðx1 ,x2 Þ-jDðvÞjðt,x1 ,x2 Þ (right) at different stages of the deformation (t = 0, t = T/3, t =2T/3, t= T).
In Fig. 10 (left) is plotted the distribution of the initial orientation of the four grains multicrystal ðx1 ,x2 Þ-yð0,x1 ,x2 Þ, given by the angle between ½1 1 0 and the Ox1 axis. Specifically, the bottom left grain, denoted by 1, is orientated at y ¼ p=8, the bottom right grain, denoted by 2, is orientated with y ¼ p=8, top right grain, denoted by 3 is orientated at y ¼ 0 and the left top grain, denoted by 4, is orientated with y ¼ p=4. The evolution of the lattice orientation of the multicrystal is plotted in Fig. 10 (left). Note that at the final stage of the deformation (t= T), in the grains 1 and 2 the orientation is almost homogeneous while each of the grains 3 and 4 contain two regions with different orientations. It means that during deformation, two subgrains have been created. Moreover, a part of the boundary between the grains 2 and 3 is no longer a discontinuity line since the two grains have the same the lattice orientation. Thus, at the end of the deformation process, the multicrystal is made of six grains. Some of the boundaries between these six grains are the initial grain boundaries (plotted in red in Fig. 10) but the other ones are new interfaces created during the deformation process. To explain this phenomenon, in Fig. 11 are plotted the evolution in time of the three slip rates ðx1 ,x2 Þ-g_ r ðt,x1 ,x2 Þ distributions in the multicrystal. In general, the shear bands corresponding to a given slip system do not cross the grain boundaries. The only exception is the slip system r = 2, since some shear bands cross a part of the boundary between the
ARTICLE IN PRESS 858
O. Cazacu, I.R. Ionescu / J. Mech. Phys. Solids 58 (2010) 844–859
Fig. 11. The distribution of the slip rates ðx1 ,x2 Þ-g_ r ðt,x1 ,x2 Þ (left: r = 1, middle: r =2, right: r= 3) at different stages of the deformation (t= T/3, t = 2T/3, t= T).
grains 2 and 3. The occurrence of these shear bands may explain why the grains 2 and 3 have the same orientation on both sides of their common boundary.
9. Conclusions In this paper, we have developed an Eulerian rate-dependent single crystal model applicable at high-strain rates, large strains, and large material and lattice rotations. Among the novel ingredients of the model is the use of an Eulerian form for the evolution equations of the slip directions and slip planes (i.e. involving material time derivatives). Also, by adopting an overstress based extension of the Schmid law, motivated by the microdynamics of crystal defects, numerical instabilities associated with the classical Norton law are eliminated. Although the proposed model and associated numerical strategy are fully 3-D, of particular interest is the in-plane response. It was shown that for the proposed rigid-viscoplastic crystal model, the same 3-D slip systems as in the inviscid case (see Rice, 1973) need to be simultaneously activated in order to accommodate plane strain deformation. Furthermore, for the plane strain case, the proposed rigid viscoplastic law can be written as a relation between a two-dimensional stress tensor and a two-dimensional rate of deformation tensor (see Eq. (34)). This 2-D flow rule has the same form as in the 3D case, with in-plane CRSS and in-plane viscosity coefficients being obtained by multiplying the respective 3D properties with some factors that depend on the lattice geometry (see Eq. (31)). Also, it was shown that in the plane strain case, the lattice evolution is described by a single differential equation. Furthermore, for an FCC single crystal, analytic expressions for the slip rates on individual composite systems have been obtained. Specifically, for a given lattice orientation, the decomposition of the in-plane rate of deformation into the three in-plane composite slip systems are given by Eqs. (45)–(49). The in-plane model was applied to predict constitutive response in equal die extrusion and channel die compression of single FCC crystals and multicrystals. The two example problems selected illustrate one of the main advantages of using an Eulerian description, which is the evaluation of the basic characteristics of the orientation flow during plastic deformations. Since both problems involve plane strain deformation, the computational cost is low. To integrate the governing equations, a mixed finite-element and finite-volume strategy, developed by Cazacu and Ionescu (2010), was used. One of the main advantages of using an Eulerian setting is that it is sufficient to do the computations of the flow only in a region around the extruder angle, whereas in a Lagrangian description the computations need to be done in all the material points of the sample. Our results have shown that the initial crystal orientation has a significant effect on the lattice rotation. Analysis of the evolution of the lattice orientation of a single square grain embedded in another grain has shown that the parent crystal’s initial orientation and rotation do not influence the lattice rotation of the inserted grain. Yet, the parent grain influences the deformation of the inserted grain. The viscosity plays an important role in strain localization during channel die compression. For low viscosity, shear bands parallel to a slip direction are present from the beginning of the process. The lattice has an important rotation only in these shear bands. For larger viscosity, localization still takes place but is much more diffuse. At high crosshead velocity, the plastic properties and crystal anisotropy are less important while inertia effects are dominant. Finally, grain interactions were analyzed by simulating compression of a four grains multicrystal. It was found that during deformation, sub-grains have been created. Moreover, a part of the boundary between certain grains is no longer a discontinuity line: the two adjacent grains have the same final lattice orientation. This phenomenon can be explained by the formation of shear bands, associated to an individual slip system, that cross the grain boundary. Thus, at the end of the deformation process, some of the boundaries between grains are initial grain boundaries but other ones are new interfaces created during the deformation process. The proposed Eulerian rate-dependent single crystal model is particularly suitable for solving boundary-value problems involving large strains and strain rates. In this paper, the model was applied to FCC single crystals and multicrystals. Future plans include the application of the model to describe the single-grain behavior in a polycrystal along with texture evolution.
ARTICLE IN PRESS O. Cazacu, I.R. Ionescu / J. Mech. Phys. Solids 58 (2010) 844–859
859
Acknowledgment The authors would like to thank Drs. Cristian Teodosiu and Patrick Franciosi for insightful comments on the manuscript. Partial financial support for this work through grant FA8651-08-1-0009 is gratefully acknowledged. References Anand, L., Kothari, M., 1996. A computational procedure for rate-independent crystal plasticity. J. Mech. Phys. Solids 44, 525–558. Asaro, R.J., Rice, J.R., 1977. Strain localization in ductile single crystals. J. Mech. Phys. Solids 25, 309–338. Asaro, R.J., Needleman, A., 1985. Texture development and strain hardening in rate dependent polycrystals. Acta Metall. 33, 923–953. Busso, E.P., Cailletaud, G., 2005. On the selection of active slip systems in crystal plasticity. Int. J. Plasticity 12, 1–28. Cazacu, O., Ionescu, I.R., 2010. An augmented Lagrangian method in Eulerian modeling of visco-plastic crystals. Comput. Methods Appl. Mech. Eng. 199, 689–699. Chang, Y.W., Asaro, R.J., 1981. An experimental study of shear localization in aluminum–copper single crystals. Acta Metall. 29, 241–250. Crone, W.C., Shield, T.W., Creuziger, A., Henneman, B., 2004. Orientation dependence of the plastic slip near notches in ductile FCC single crystals. J. Mech. Phys. Solids 52, 85–112. Dawson, P.R., 2004. Crystal plasticity. In: Raabe, D. Roters, F., Barlat, F., Chen, L.-Q. (Eds.), Continuum Scale Simulations of Engineering Materials. Wiley-VCH, pp. 115–141. Delannay, L., Jacques, P.J., Kalidindi, S.R., 2006. Finite element modeling of crystal plasticity with grains shaped as truncated octahedrons. Int. J. Plast. 22, 1879–1898. Fortin, M., Glowinski, R., 1982. Methodes de Lagrangien augmente´s, application a la re´solution de proble mes aux limites, Dunod. Franciosi, P., 1984. Etude the´orique et expe´rimentale du comportement e´lastoplastique de monocristaux me´talliques se de´formant par glissement: Mode´lisation pour un chargement complexe quasi-statique., Ph.D. Thesis, Paris 13 University, France. Franciosi, P., 1985. The concepts of latent hardening and strain hardening in metallic single crystals. Acta Metall. 33, 1601–1612. Gan, Y.X., Kysar, J.W., Morse, T.L., 2006. Cylindrical void in a rigid-ideally plastic single crystal II: experiments and simulations. Int. J. Plast. 22, 39–72. Glowinski, R., Le Tallec, P., 1989. Augmented Lagrangian and operator splitting method in non-linear mechanics. SIAM Studies in Applied Mathematics. Hafner, K.S., 1992. Finite Plastic Deformation of Crystalline Solids. Cambridge University Press, pp. 110–157. Hill, R., Rice, J.R., 1972. Constitutive analysis of elastic–plastic crystals at arbitrary strain. J. Mech. Phys. Solids 20, 401–413. Hutchinson, J.W., 1976. Bounds and selfconsistent estimates for creep of polycrystalline materials. Proc. R. Soc. Lond. A 348, 101–127. Kalidindi, S.R., Donohue, B.R., Li, S., 2009. Modeling texture evolution in equal channel angular extrusion using crystal plasticity finite element models. Int. J. Plast. 25, 768–779. Kok, S., Beaudoin, A.J., Tortorelli, D.A., 2002. A polycrystal plasticity model based on the mechanical threshold. Int. J. Plast. 18, 715–741. Kuchnicki, S.N., Radovitzky, R.A., Cuitino, A.M., 2008. An explicit formulation for multiscale modeling of bcc metals. Int. J. Plast. 24, 2173–2191. Lebensohn, R.A., Tome´, C.N., 1993. A self-consistent anisotropic approach for the simulation of plastic deformation and texture development of polycrystals: application to zirconium alloys. Acta Metall. Mater. 41, 2611–2624. Ling, X., Horstemeyer, M.F., Potirniche, G.P., 2005. On the numerical implementation of 3D rate-dependent single crystal plasticity formulations. Int. J. Numer. Methods Eng. 63, 548–568. Li, S., Beyerlein, I.J., Alexander, D.J., Vogel, S.C., 2005. Texture evolution during multi-pass equal channel angular extrusion of copper: neutron diffraction characterization and polycrystal modeling. Acta Mater. 53, 2111–2125. Li, S., Beyerlein, I.J., Necker, C.T., 2006. On the development of microstructure and texture heterogeneity in ECAE via Route C. Acta Mater. 54, 1397–1408. Luri, R., Luis, C.J., Leo´n, J., 2006. A new configuration for equal channel angular extrusion dies. Trans. ASME 128, 860–865. Maurice, C.L., Driver, J.H., 1993. High temperature plane strain compression of cube oriented aluminum single crystal. Acta Metall. Mater. 41, 1635–1664. McGinty, R., McDowell, D., 2006. A semi-implicit integration scheme for rate independent finite crystal plasticity. Int. J. Plast. 22, 996–1025. Peirce, D., Asaro, R.J., Needleman, A., 1982. An analysis of nonuniform and localized deformation in ductile single crystals. Acta Metall. 30, 1087–1119. Pironneau, O., 1989. Finite Element Methods for Fluids. John Wiley & Sons Ltd., Chichester. Rice, J.R., 1971. Inelastic constitutive relations for solids: an internal variable theory and its application to metal plasticity. J. Mech. Phys. Solids 19, 433–455. Rice, J.R., 1973. Plane strain slip line theory for anisotropic rigid/plastic materials. J. Mech. Phys. Solids 21, 63–74. Rousselier, G., Leclerc, S., 2006. A simplified polycrystalline model for viscoplastic and damage finite element analyses. International Journal of Plasticity 22, 685–712. Rosochowski, A., Olejnik, L., 2002. Numerical and physical modelling of plastic deformation in 2-turn equal channel angular extrusion. Journal of Materials Processing Technology 125, 309–316. Schmidt-Baldassari, M., 2003. Numerical concepts for rate-independent single crystal plasticity. Comp. Methods Appl. Mech. Engng. 192, 1261–1280. Segal, V.M., 1981. Plastic Working of Metals by Simple Shear. Russ. Metall. 1, 99–105. Simo, J.C., Hughes, T.J.R., 1998. Computational Inelasticity. Springer-Verlag, New York. Tabourot, L., 1992. Loi de comportement e´lastoviscoplastique du monocristal en grandes de´ formations, Ph.D. Thesis, Grenoble University, France. Taylor, G.I., Elam, C.F., 1923. The distortion of an aluminum crystal during a tensile test. Proc. Royal Soc. London A 102, 643–647. Taylor, G.I., Elam, C.F., 1925. The plastic extension and fracture of aluminum single crystals. Proc. Royal Soc. London A 108, 28–51. Taylor, G.I., 1938a. Plastic strains in metals. J. Institute of Metals 62, 307–324. Taylor, G.I., 1938b. Analysis of plastic strain in a cubic crystal. In: Stephen Timoshenko 60th Anniversary Volume. McMillan Co., New York, pp. 218–224. Teodosiu, C., Sidoroff, F., 1976. A theory of finite elastoviscoplasticity of single crystals. Int. J. Engng. Sci. 14, 165–176. Teodosiu, C., Raphanel, J.L., Tabourot, L., 1993. Finite element simulations of the large elastoplastic deformation of multicrystals. In : Teodosiu, C., Raphanel, J.L., Sidoroff, F. (Eds.), Large Plastic Deformations: Fundamental Aspects and Applications to Metal Forming. Balkema, Rotterdam, pp. 153–168. Temam, R., 1979. Navier–Stokes Equations. Theory and Numerical Analysis. North-Holland, Amsterdam. Tome, C.N., Lebensohn, R.A., 2003. Manual for Code Visco-plastic self-consistent (VPSC) Version 6. Toth, L.S., Massion, R.A., Germain, L., Baik, S.C., Suwas, S., 2004. Analysis of texture evolution in equal channel angular extrusion of copper using a new flow field. Acta Mater. 52 (7), 1885–1898. Van Houtte, P., Li, S., Engler, O., 2004. Taylor type homogenization methods for texture and anisotropy. In: Raabe, D., Roters, F., Barlat, F., Chen, L.-Q. (Eds.), Continuum Scale Simulations of Engineering Materials. Wiley-VCH, pp. 459–471.
