Numerical simulation of yield loci for PST crystals of TiAl

Numerical simulation of yield loci for PST crystals of TiAl

Materials Science and Engineering A239 – 240 (1997) 790 – 803 Numerical simulation of yield loci for PST crystals of TiAl S.M. Schlo¨gl a,*, F.D. Fis...

467KB Sizes 0 Downloads 43 Views

Materials Science and Engineering A239 – 240 (1997) 790 – 803

Numerical simulation of yield loci for PST crystals of TiAl S.M. Schlo¨gl a,*, F.D. Fischer b a

Christian Doppler Laboratory for Micromechanics of Materials, Uni6ersity for Mining and Metallurgy Leoben, Franz-Josef-Straße 18, A-8700 Leoben, Austria b Institute of Mechanics, Uni6ersity for Mining and Metallurgy Leoben, Franz-Josef-Straße 18, A-8700 Leoben, Austria

Abstract The yield strength of polysynthetically twinned (PST) crystals of TiAl depends on the orientation angle between the lamellar boundaries and the load axis. It is the aim of this work to model this anisotropic behaviour. A micromechanical model is employed which considers the lamellar microstructure and the deformation mechanisms such as ordinary slip and deformation twinning in g-TiAl. Since the calculated and measured yield strength data show good agreement in the case of uniaxial loading, there is a real motivation to study numerically the behaviour of a PST crystal under biaxial loading. Yield points are calculated for various proportional stress paths defined by the ratio of the principal stresses s1/s2 in relation to the orientations of the load axes towards the lamellae. These predicted yield loci are compared with the von Mises and the Hill yield criteria. © 1997 Elsevier Science S.A. Keywords: g-TiAl; a2-Ti3Al; Crystal plasticity; Deformation twinning; Yield loci; Hill yield criterion

1. Introduction

2. Micromechanical model

Intermetallic titanium aluminides are potential new high temperature structural materials. Unfortunately, they suffer from poor ductility and low toughness at ambient temperatures which appear to be the most serious obstacle to their full utilization. Especially single-phase g-alloys [Ti(50 – 56)Al] show very low room temperature ductility. However, ductility can be improved in two-phase compounds with an Al content of 48%. These materials possess a two-phase microstructure composed of g-TiAl (L10-structure, cubic) and 5 –10% a2-Ti3Al (DO19-structure, hexagonal). By means of various heat treatments it is possible to design different microstructures such as fully-lamellar, nearly-lamellar, duplex and near g [1]. In this paper attention is devoted to the lamellar microstructure, because fine grained fully lamellar g-TiAl based alloys are of main interest technologically. The aim of this work is to create a basis for finding a constitutive law for a lamellar polycrystal.

The lamellar structure of a single grain is formed according to the orientation relationships {111}g (0001)a 2 and Ž110g Ž112( 0a 2. As a consequence, six different types of ordered domains exist in the g-phase in the lamellar structure [2]. The effects of the lamellar structure can be best studied with so called polysynthetically twinned (PST) crystals of TiAl. These are single-crystal-like ingots with a structure composed of a single set of lamellae with a specific orientation [3]. The deformation behaviour of PST crystals at room temperature was investigated systematically [3,4]. It is well known that the plastic behaviour of a PST crystal of TiAl is highly anisotropic: under uniaxial loading at room temperature its yield strength depends on the orientation angle f between the load axis and the

* Corresponding author. Present address: Delft University of Technology, Laboratory for Engineering Mechanics, 2628 CD Delft, The Netherlands. 0921-5093/97/$17.00 © 1997 Elsevier Science S.A. All rights reserved. PII S 0 9 2 1 - 5 0 9 3 ( 9 7 ) 0 0 6 6 8 - 0

Table 1 Slip and twinning systems in the g-phase of TiAl-intermetallics Plane

Slip direction

Twinning direction

(111) (1( 11) (11( 1) (111( )

[1( 10] [110] [110] [1( 10]

[112( ] [1( 12( ] [11( 2( ] [112]

S.M. Schlo¨gl, F.D. Fischer / Materials Science and Engineering A239–240 (1997) 790–803

Fig. 1. Unit cell for biaxial loading of a PST crystal in the (x, z)plane. The angles f1 and f2 determine the orientation of the two principal stress axes in reference to the lamellae (f1 +f2 =90°). Definition of the (x¯, y¯, z¯ )-coordinate system.

lamellar boundaries [5]. The main deformation mechanisms in the g-phase are ordinary slip on the {111} plane along Ž11( 0] and true twinning of the 1 {111} Ž111( ]-type (Table 1). Possible slip systems in the 6 a2-Ti3Al phase are the prismatic {101( 0}Ž12( 10 and pryramidal {112( 1}ŽI( 1( 26 slip systems with very different values of critical shear stress. In this work slip in g-TiAl as well as in a2-Ti3Al is modelled in the framework of crystal plasticity. This concept is based on Schmid’s law and has been adapted for numerical applications by Asaro and colleagues [6 –8]. In this rate dependent formulation the shear rate on the slip system a is determined by the resolved shear stress t (a) and the slip system hardness g (a). Deformation twinning is described with a critical shear stress criterion analogous to Schmid’s law for crystallographic slip. This means that the shear stress t (b) across the twinning plane, resolved in the twinning direction, is the driving force for the nucleation and growth of the twin b. The constitutive equations used for slip and deformation twinning and the chosen material parameters are described in detail in [9,10]. Our micromechanical model is based on the unit cell technique using the finite element method ABAQUS [11]. A three-dimensional unit cell with periodic boundary conditions reflects the lamellar microstructure (Fig. 1). An amount of the a2-Ti3Al phase of j= 10% is studied. All six ordered domains of the g-phase are taken into account, namely the matrix domains I (g 1), II (g 2) and III (g 3), and the twin domains I, II and III (g 4, g 5 and g 6). Crystal plasticity is implemented into the finite element code by the special user supplied material routine UMATCRYSPL initially developed by Huang [12] for cubic crystals. We

791

extended this subroutine to describe crystallographic slip in hexagonal crystals and deformation twinning in the cubic crystal system [10]. This approach was applied to simulate room temperature uniaxial tests for orientation angles f varying in the range of 05 f590°. Since this simulation has already been described in detail in the paper by Schlo¨gl and Fischer [9], only the results are presented here. Fig. 2a shows the calculated stress–strain curves for the selected orientations. From these curves the stress at 0.2% offset strain, s0.002, was determined as a function of f, and compared with experimentally obtained values [5]. Fig. 2b confirms the excellent agreement between experiment and numerical simulation. Hence, the simulation also reflects the anisotropic behaviour of the PST crystal of TiAl. Therefore, we can rely on this model to predict the behaviour of a PST crystal under biaxial loading. It is quite interesting to point out the influence of the orientation of the lamellae on the initial yield locus under proportional loading, since the behaviour of PST crystals of TiAl under multiaxial loading has not been investigated experimentally yet.

Fig. 2. Influence of the angle f on the calculated flow curves of a PST crystal of TiAl. (a) Computed stress – strain curves in tension ( + ) and compression ( − ) [9]. (b) s0.002 As a function of f as obtained from experiments [5] and numerical simulation [9].

792

S.M. Schlo¨gl, F.D. Fischer / Materials Science and Engineering A239–240 (1997) 790–803

process time ti, corresponding to a (s1, s2)-pair, the finite element model supplies the cumulative shear strain g on the slip and twinning systems a, which occurs on an average in the whole unit cell: g=% a

Fig. 3. Applied principal stresses s1 and s2 vs. process time for a r = 0.5 biaxial test. The vertical lines t1, t2 and t3 refer to the times of 0.001, 0.003 and 0.005 cumulative shear strain g.

3. Simulation of biaxial tests Yielding in a homogenous and isotropic material can be characterized by the principal stresses s1, s2 and s3. The most common yield criterion for such a material is the von Mises criterion:

(1/2)[(s1 −s2)2 +(s2 −s3)2 +(s3 −s1)2]= sf

(1)

with sf being the uniaxial yield strength. In the case of an anisotropic material like the investigated PST crystal, it is not sufficient to know the values of the principal stresses. Its response to biaxial loading also depends on the directions of the principal stresses with respect to the lamellae. The first step is to fix the orientations of the load axes with respect to the lamellae. The presented simulations are limited to plane stress problems where one principal stress is zero. The other two principal stress axes lie in the (x, z)-plane of the unit cell, and their orientations are characterized by the angles f1 and f2 (Fig. 1) whose sum (f1 +f2) must be equal to 90°. The micromechanical model and the material parameter set 1 presented in [9] are taken to determine the plane stress yield locus of a PST crystal of TiAl. The first calculations are performed for the orientation f1 = 0°, f2 = 90°. The unit cell of Fig. 1 is deformed along a proportional stress path defined by the ratio r of the principal stresses s1 and s2 (r = s1/s2). This unit cell is employed for all calculations, no matter which orientation (characterized by f1 and f2) is investigated. Since the periodic boundary conditions only allow for sxx, syy, szz, sxy, sxz and syz, the prescribed (s1, s2)stress state has to be transformed into the (x, y, z)-coordinate system. To predict the yield point along such a proportional stress path, the principal stresses s1 and s2 (in fact, the corresponding sxx, szz and sxz stresses) are increased with time (which is only a ‘process’ time) by keeping the ratio r = s1/s2 constant (Fig. 3). For each

&

t

g; (a) dt,

0

where g; (a) is the shearing rate on the slip (twinning) system a measured relative to the lattice. Since slip and deformation twinning result in a plastic deformation, the cumulative shear strain on the slip and twinning systems is taken to evaluate the yield locus. It is usually common in continuum mechanics to use the concept of effective stress and effective strain to reduce a multiaxial response to a universal stress–strain curve, from which the effective stress at various prescribed offset strains (for example 0.1 or 0.2%) can be determine [13]. However, in the presented case the cumulative shear strain g is calculated in the framework of crystal plasticity. This knowledge is sufficient to compare multiaxial stress states. Various offset-g values for yielding are investigated: g=0.001, g= 0.003 and g= 0.005. Fig. 4 shows the evolution of the computed cumulative shear strain for a r=0.5 biaxial test. The vertical lines in Fig. 4 indicate the times corresponding to 0.001, 0.003 and 0.005 cumulated shear strain. The coordinates of a single point on the yield locus can be determined in Fig. 3 by using the time found in Fig. 4. The yield points, corresponding to these three offset-g values, are now known for the proportional stress path s1/s2 = 0.5 and can be depicted in a (s1, s2)-plane. The whole procedure is repeated for other proportional stress paths: “ First quadrant: r= 0, r= 0.5, r=6/7, r=1, r=7/ 6, r=2, r= 4, r= , “ Fourth quadrant: r= − 0.5, r= −1, r= −2, r= − .

Fig. 4. Computed cumulative shear strain g on the slip and twinning systems averaged over the whole unit cell as a function of process time for a r= 0.5 biaxial test with the orientation f1 = 0°, f2 = 90°.

S.M. Schlo¨gl, F.D. Fischer / Materials Science and Engineering A239–240 (1997) 790–803

793

Fig. 5. Computed yield loci at 0.001, 0.003 and 0.005 cumulative shear strain g on all slip and twinning systems averaged over the whole unit cell for biaxial tests with proportional loading paths and an orientation of f1 =0°, f2 =90°.

794

S.M. Schlo¨gl, F.D. Fischer / Materials Science and Engineering A239–240 (1997) 790–803

Fig. 6. Computed yield loci at 0.001, 0.003 and 0.005 cumulative shear strain g on all slip and twinning systems averaged over the whole unit cell for biaxial tests with proportional loading paths and an orientation of f1 =15°, f2 =75°.

S.M. Schlo¨gl, F.D. Fischer / Materials Science and Engineering A239–240 (1997) 790–803

795

Fig. 7. Computed yield loci at 0.001, 0.003 and 0.005 cumulative shear strain g on all slip and twinning systems averaged over the whole unit cell for biaxial tests with proportional loading paths and an orientation of f1 =30°, f2 =60°.

796

S.M. Schlo¨gl, F.D. Fischer / Materials Science and Engineering A239–240 (1997) 790–803

Fig. 8. Computed yield loci at 0.001, 0.003 and 0.005 cumulative shear strain g on all slip and twinning systems averaged over the whole unit cell for biaxial tests with proportional loading paths and an orientation of f1 =45°, f2 =45°.

S.M. Schlo¨gl, F.D. Fischer / Materials Science and Engineering A239–240 (1997) 790–803

797

Additionally, the yield loci are also calculated in the second and third quadrant of Fig. 5 (f1 = 0°, f2 =90°). The overall results for f1 = 0°, f2 = 90° are illustrated in the so called ‘pi-plane’. This is a plane in the (s1, s2, s3)-space which is normal to the space diagonal. The pi-plane shows the cross-section of the yield locus normal to the space diagonal. Any three-dimensional stress vector can be projected onto two orthogonal axes in the pi-plane. Its coordinates x1 and x2 are given by [14]: x1 = 1/2(s2 − s1), x2 = 1/6(2s3 − s1 −s2). These relations are taken to transfer the yield loci predicted for the stress subspace with s3 = 0 onto the pi-plane. The resulting yield loci for g= 0.001, g=0.003 and g=0.005 are shown in Fig. 9. These curves characterize the shape of the cross-section of the yield loci. Since the von Mises criterion describes a circle in the pi-plane, the anisotropy of a PST crystal becomes evident from Fig. 9. Fig. 9. Computed yield loci at g = 0.001, g = 0.003 and g =0.005 for the orientation of f1 = 0°, f2 = 90° depicted in the pi-plane.

The simulation is confined to the first and fourth quadrant of the (s1, s2)-plane because of the large CPU time to calculate a single yield point. The predicted yield points are plotted in the (s1, s2)-plane, and the point possessing the same cumulative shear strain g are connected by straight lines to construct the yield locus. Fig. 5 shows the computed yield loci of a PST crystal under biaxial testing with principal stresses orientated parallel and perpendicular to the lamellae (f1 = 0°, f2 = 90°). The inner curve refers to a cumulative shear strain g of 0.001, the middle one to a g of 0.003 and the outer one to 0.005. The shapes of these yield loci differ, which indicates that the hardening depends on the chosen stress path. Since a PST crystal of TiAl behaves in an anisotropic way, its response to biaxial loading has to be investigated for other orientations, too. The above mentioned procedure is repeated to evaluate the yield loci for principal stress axes inclined by 15° and 75° to the lamellae (f1 =15°, f2 =75°). Fig. 6 demonstrates the resulting yield loci in the first and fourth quadrant of the (s1, s2)-plane. Fig. 7 shows the results of the calculations for the orientation f1 =30°, f2 =60°. Fig. 8 presents the computed yield loci for the orientation f1 = 45°, f2 = 45°. By comparing Figs. 5 – 8 one can observe that the orientation of the lamellae greatly influences the response of a PST crystal to biaxial loading. As expected, the two extreme orientations are f1 = 0°, f2 = 90° (hard orientation) and f1 =45°, f2 = 45° (soft orientation). It is clear, that in the case of s1 = s2 = 1 (r= 1) different orientations of the lamellae have no effect.

3.1. Comparison with common yield criteria A yield function is now looked for, which can describe the predicted yield loci of Figs. 5–8. First, we compare the yield loci predicted from crystal plasticity with the von Mises criterion. The yield locus at g= 0.005 is defined as the initial one, since it nearly corresponds to the stress at 0.2% offset strain, s0.002, which is usually determined in experiments. So the yield loci at g= 0.005 are taken from Figs. 5–8 and compared with the yield surface one obtains by solving Eq. (1). The measured uniaxial yield strength sf is not a constant; it depends on the orientation of the lamellae towards the load axis (characterized by f). Fig. 10 shows the yield locus for the orientation f1 = 0°, f2 = 90° and the von Mises ellipse with ff = 343 MPa (which corresponds to s0.002 for uniaxial loading with an orientation angle f= 0°). For the orientation f1 = 15°, f2 = 75° Eq. (1) is solved with ff = 232 MPa, and this ellipse is depicted together with the computed yield locus at g=0.005 in Fig. 11. To draw the von Mises ellipses for the orientation f1 = 30°, f2 = 60° and f1 = 45°, f2 = 45°, ff is taken as 135 and 120 MPa, respectively. The results can be seen in Figs. 12 and 13. Table 2 Yield strenghts X, Y, Z calculated with the crystal plasticity model for tension and compression Yield strength

X (MPa)

Y (MPa)

Z (MPa)

In tension In compression Average

343 311 327

287 372 329.5

430 450 440

The values stem from Fig. 9 and are used as material constants for the Hill criterion.

798

S.M. Schlo¨gl, F.D. Fischer / Materials Science and Engineering A239–240 (1997) 790–803

Fig. 10. Crystal plasticity yield locus at g = 0.005, von Mises yield criterion with sf =343 MPa and Hill yield criterion for biaxial tests with proportional loading paths and an orientation of f1 = 0°, f2 =90°.

S.M. Schlo¨gl, F.D. Fischer / Materials Science and Engineering A239–240 (1997) 790–803

799

Fig. 11. Crystal plasticity yield locus at g = 0.005, von Mises yield criterion with sf =232 MPa and Hill yield criterion for biaxial tests with proportional loading paths and an orientation of f1 = 15°, f2 =75°.

800

S.M. Schlo¨gl, F.D. Fischer / Materials Science and Engineering A239–240 (1997) 790–803

Fig. 12. Crystal plasticity yield locus at g = 0.005, von Mises yield criterion with sf =135 MPa and Hill yield criterion for biaxial tests with proportional loading paths and an orientation of f1 = 30°, f2 =60°.

S.M. Schlo¨gl, F.D. Fischer / Materials Science and Engineering A239–240 (1997) 790–803

801

Fig. 13. Crystal plasticity yield locus at g = 0.005, von Mises yield criterion with ff =120 MPa and Hill yield criterion for biaxial tests with proportional loading paths and an orientation of f1 = 45°, f2 =45°.

802

S.M. Schlo¨gl, F.D. Fischer / Materials Science and Engineering A239–240 (1997) 790–803

It turns out, that the von Mises criterion is a rather unsatisfying approximation for the response of a lamellar grain to multiaxial loading. Hence, for future simulations of polycrystalline lamellar g-TiAl based alloys it is necessary to look for a better approach. Hill [15,16] proposed a yield criterion for anisotropic materials which is a generalization of the von Mises yield criterion, F(sy¯y¯ − sz¯z¯ )2 +G(sz¯z¯ −sx¯x¯ )2 +H(sx¯x¯ −sy¯y¯ )2 +2Ls 2y¯z¯ + 2Ms 2x¯z¯ +2Nt 2y¯z¯ =1.

(2)

F, G, H, L, M and N are material constants. The x¯ -, y¯ - and z¯ -axes, which are called the principal axes of anisotropy, are assumed to be axes of two-fold rotational symmetry. The parameters F, G and H are related to the yield strength in the x¯ -direction (X), the yield strength in y¯ -direction (Y) and the yield strength in z¯ -direction (Z), respectively, and the other parameters L, M and N can be determined from shear tests [17]. It is interesting to figure out which yield loci are predicted by the Hill criterion. First, the x¯ -, y¯ -, z¯ axes have to be fixed in the unit cell (Fig. 1). x¯ is assumed to be parallel to the lamellae of the PST crystals; it corresponds to the [112( ] direction of the matrix domain g 1-TiAl. y¯ also lies in the plane of the lamellae but is parallel to the [11( 0] direction of g 1TiAl. z¯ is perpendicular to the lamellae (Fig. 1). For plane stress loading (sy¯y¯ =sy¯z¯ =sx¯y¯ =0) Eq. (2) reduces to Fs 2z¯z¯ + G(sz¯z¯ −sx¯x¯ )2 +Hs 2x¯x¯ +2Ms 2x¯z¯ =1.

(3)

To evaluate the material parameters F, G and H, the calculated values for the yield strengths X, Y, Z can be read from the yield loci at g =0.005 in Fig. 9. Since the principal axis 1 of Fig. 9 corresponds to x¯, axis 2 to z¯ and axis 3 to y¯, the yield points of these three axes can be taken as the yield strengths X, Y, Z. The uniaxial yield strengths are listed in Table 2; they differ for tension and compression. Since Eq. (3) only involves quadratic terms, the Hill criterion can only describe yielding that is insensitive to the sign of the applied stress. Therefore, X, Y, Z receive the average value for tension and compression (Table 2), and the following material constants are used: F =2.511 · 10 − 6, G= 2.653 · 10 − 6 and H =6.698 · 10 − 6 MPa − 2. There is still one material constant missing, namely M, which is related to the yield point of the stress state sx¯z¯ " 0, sx¯x¯ =sy¯y¯ =sz¯z¯ =sx¯y¯ =sy¯z¯ = 0 by Ms 2x¯z¯ =1/2. Its principal stresses are s1 = +sx¯y¯, s2 = − sx¯y¯ and s3 = 0 with principal stress axes inclined by 45° to the (x¯, y¯, z¯ )-coordinate system. Such a stress state has already been simulated by a stress path with r= − 1 for the orientation f1 =45°, f2 =45° (Fig. 8). The micromechanical model states that yielding

occurs at s1 = 61 MPa and s2 = − 61 MPa, which leads to an M parameter of 1.34 · 10 − 4 MPa − 2. The chosen material constants fulfil the two inequalities M\ (H+ 2G) and M\ (F+ 2G). This is the prerequisite [17] that the Hill criterion can describe the marked orientation dependence of the uniaxial yield strength: the maximum in yield strength occurs at f= 0° and 90°, the minimum at f= 45°. With these material parameters, Eq. (3) and various transformation matrices for switching between the different coordinate systems (x¯, y¯, z¯ )“ (1, 2, 3), it is possible to calculate the yield loci according to Hill for biaxial loading in various directions. The results are recorded in Figs. 10–13. While it is true that these curves differ a little from the ones found with the more sophisticated micromechanical model, it is also obvious that the Hill criterion is a much better approximation than the von Mises criterion. Because of the long CPU time that exists today it is impractical to model each lamellar grain in a polycrystal within the framework of crystal plasticity. It does make sense, however, to model a lamellar grain with the Hill criterion (Figs. 10–13).

4. Conclusions This paper deals with the micromechanical study of the flow behaviour of PST crystals of TiAl at room temperature. These crystals are single-crystal-like ingots composed of a single set of parallel g-TiAl and a2-Ti3Al lamellae. A micromechanical model is presented which considers crystallographic slip and deformation twinning as deformation mechanisms. This model is based on the unit cell technique using the finite element method. The microstructure and the crystallography are incorporated in the framework crystal plasticity. As a first step, the flow behaviour of a PST crystal of TiAl is simulated for uniaxial loading at room temperature and checked by experimental results known from the literature [5]. The simulation of uniaxial tests reflects the observed anisotropic plastic behaviour of a single lamellar grain. Then the model is employed for biaxial loading. Yield loci are predicted for various orientations. The findings can be stated as follows: (1) The orientations of the two load axes play a significant role in the flow behaviour. There is a great difference between the yield locus of the ‘hard orientation’ (load axes are parallel and perpendicular to the lamellae) and of the ‘soft’ orientation (both load axes are inclined by 45° to the lamellae). (2) The yield loci of different orientations predicted by a more sophisticated crystallographic model can be described approximately by the Hill criterion for anisotropic materials.

S.M. Schlo¨gl, F.D. Fischer / Materials Science and Engineering A239–240 (1997) 790–803

References [1] Y.-W. Kim, J. Met. 41 (1989) 24–30. [2] H. Inui, A. Nakamura, M.H. Oh, M. Yamaguchi, Philos. Mag. A 66 (1992) 557 –573. [3] T. Fujiwara, A. Nakamura, M. Hosomi, S.R. Nishitani, Y. Shirai, M. Yamaguchi, Philos. Mag. A 61 (1994) 591 – 606. [4] H. Inui, M.H. Oh, A. Nakamura, M. Yamaguchi, Acta Metall. Mater. 40 (1992) 3095–3104. [5] M. Yamaguchi, H. Inui, in: C.T. Liu, R.W. Cahn, G. Sauthoff (Eds.), Ordered Intermetallics—Physical Metallurgy and Mechanical Behaviour, Kluwer, Dordrecht, 1992, pp. 217 – 235. [6] R.J. Asaro, Adv. Appl. Mech. 23 (1983) 1–115. [7] R.J. Asaro, A. Needleman, Acta Metall. 33 (1985) 923 – 953. [8] D. Peirce, R.J. Asaro, A. Needleman, Acta Metall. 31 (1983) 1951 – 1976.

.

803

[9] S.M. Schlo¨gl, F.D. Fischer, Philos. Mag. A 75 (1997) 621–636. [10] S.M. Schlo¨gl, Micromechanical modelling of the deformation behaviour of gamma titanium aluminide, Thesis, Montanuniversita¨t Leoben, Austria, 1997. [11] Hibbit, Karlsson and Sorensen, ABAQUS User manual version 5.5, Pawtucket, RI, 1995. [12] Y.G. Huang, Harvard University Report MECH-178, 1991. [13] M.G. Stout, S.S. Hecker, R. Bourcier, ASME J. Eng. Mater. Technol. 105 (1983) 242 – 249. [14] D.C. Stouffer, L.T. Dame, Inelastic Deformation of Metals, Wiley, New York, 1996. [15] R. Hill, Proc. R. Soc. London 193A (1948) 281 – 291. [16] R. Hill, The Mathematical Theory of Plasticity, chap. XII, Clarendon Press, Oxford, 1950. [17] W.F. Hosford, The Mechanics of Crystals and Textured Polycrystals, Oxford University Press, Oxford, 1993.