THE INCORPORATION OF TEXTURE-BASED YIELD LOCI INTO ELASTO-PLASTIC FINITE ELEMENT PROGRAMS

Elasto-plastic t’mite elements (FE) methods are nowadays widely used to simulate complex metal forming processes. It is then useful to generate an anisotropic yield criterion from the crystallographic texture and incorporate it into such model. The theory of dual plastic potentials (one in strain rate space and one in stress space) helps to achieve this. There is however a certain danger of losing the convexity of the yield locus during this procedure. Examples of this phenomenon are given and discussed. It is furthermore explained how the yield locus can be used to generate an elasto-plastic modulus for implementation in the FE code. Finally several examples of successful applications of the anisotropic FE code to metal forming problems are given.


INTRODUCTION
Elastic-plastic finit element models (FEM) can be used to calculate the stress, strain and temperature distributions during complex forming processes such as rolling, forging and deep drawing.A vast literature exists in this field; an extensive reference list is given by Rowe et al. (1990).For many industrial applications, a high level of accuracy is desired.This requires the use of good material models.One of the material properties which cannot be neglected is the plastic anisotropy.By this it is meant that the flow stress of the material depends on the orientation of the axes of the loading system.
In engineering materials such as steel sheet or aluminium alloys (both in annealed condition), this anisotropy is mainly determined by the crystallographic texture in the material.After subsequent plastic deformation without annealing, other factors will also contribute to the anisotropy (Teodosiu, 1992).
Initially, nearly all the FE models used the isotropic von Mises yield criterion.It was however soon realised that engineering materials used for example for can making or for car body parts are always anisotropic.So in many of these FE codes, the von Mises law has now been replaced by the anisotropic Hill criterion.In its most general form, the latter is expressed as follow: eijkl Odij O'kl { (1) in which the 0" are the components of the deviatoric stress tensor, o" 0 is the flow stress in a tensile test.The coefficients Bijkl have to be found by fitting the equation to the results of mechanical tests.There are 81 such coefficients, but not all of them are independent.The 4th order tensor is symmetric: and nijkl-njikl (2) Bun Bkt u (3) This reduces the number of independent coefficients to 21.The fact that for ordinary engineering materials only deviatoric stresses need to be considered further reduces the number of independent coefficients to 15.Finally, for materials of which the plastic properties possess symmetry, some of the coefficients can be set to zero, using the principles explained by Nye (1976).Only 6 independent coefficients are left in case of orthorhombic symmetry.Such symmetry can be expected in sheet material after rolling and/or annealing, but it would be lost as soon a more complex plastic deformation has taken place.But even the determination of 6 coefficients requires tiresome mechanical testing.Moreover, a Hill expression cannot reproduce all the features of the anisotropy of the material.An extreme example is given by Figure 1.It shows the O'11-O'22 section of the yield locus of 2 polycrystals of a f.c.c, material: one with a very sharp texture (almost a (100) [011] single crystal) and one with an orientation distribution function (ODF) which has the shape of a Gaussian distribution.Its centre is the (100) [011] orientation; the spread is characterised by an angle 0 16-5.Cubic- orthorhombic symmetry has been enforced.We have used the formulas given by Bunge (1982) in order to generate this ODF.Both yield loci are "crystallographic" yield loci; they have been generated using the "geometrical" method explained by Van Houtte el al. (1989).The yield locus of the nearly-single crystal is almost a perfect square, Figure 1 (a) t-t22 section of the yield locus of a nearly-single crystal with (100)[011] orientation.
(b) Idem for a polycrystal; the ODF (texture index: 8.9) is a Gaussian spread (D 16.5) around the same orientation.Fcc material.
as would be predicted by standard crystal plasticity theory (Backofen, 1972).Even the other yield locus section looks like a square with rounded comers.Since the Hill criterion is represented by a quadratic expression (eq. ( 1)) the best fit that can ever be expected from it is a circle.This would really be a very bad approximation, especially when one is interested in the plastic strain associated to a particular stress.According to the normality rule, a vector which represents a plastic strain tensor must be normal to the yield locus at the point that represents the stress (see for example Backofen, 1972).
At points in the vicinity of the equibiaxial stress state (O'l tr22), the vectors normal to Hill's yield locus (a circle) would almost make an angle of 45 with the vectors normal to the crystallographic yield locus.This can be improved by using more complicated expressions.Chung and Shah (1991) for example have performed finite element calculations using a more versatile anisotropic yield locus which was proposed by Barlat and Lian (1989).Another yield locus expression which could be used to such purpose has been proposed by Toth et al. (1991).In the same way as Hill's expression, these expressions contain a certain number of unknown parameters which have to be obtained by fitting to mechanical test data.As already mentioned, this often turns out to be very difficult and indeed sometimes impossible.Therefore we prefer to calculate the yield locus from the crystallographic texture of the material.The latter property is comparatively easy to measure.Such yield locus is expected to be a good approximation of the yield locus for annealed metals (which has been confirmed by MacEwen et al. (1992) for an aluminium alloy), but not necessarily for cold worked metals.The method requires a quantitative knowledge of the ODF that describes the texture (Bunge, 1982); it also requires a model for the plastic deformation of a polycrystalline material.In the next chapters, we will explain the method of calculating and representing the yield locus, followed by a description of how it can be incorporated in the FE program.
A different method (which does not require the calculation of the yield locus) has been proposed by Kalidindi and Anand (1993) and Dawson et al. (1994).In this method, the polycrystalline material is not represented by an anisotropic yield locus, but by a discrete set of "crystallites".The crystal plasticity equations are solved for each of these and then averaged each time the FE code needs information about the material behaviour.

CALCULATION OF THE YIELD LOCUS FROM THE TEXTURE
Model for the plastic deformation of polycrystals A theory is needed which allows the construction of a model for the plastic deformation of a polycrystalline material with a known texture.An overview of existing models has been given by Leffers et al. (1988) and by Aemoudt et al. (1993).The well-known theory proposed by Taylor (1938) was found to be comparatively simple and leading to satisfactory results for many engineering applications.The basic ideas of this theory (in a slightly generalised form) are the following: the strain rate tensor tij that describes the macroscopic deformation of the polycrystalline material is supposed to be known.It is the sum of an elastic part and a plastic part: The elastic part is usually neglected in the Taylor theory.In the present paper, the symbol D will also be used to denote the plastic strain rate tensor: (5) d, which is the local strain rate tensor that describes the deformation of the individual crystallites is supposed to be equal to D, and is therefore also known; the local stress acting in the individual crystallites will be calculated for each crystallite separately.The crystallites are characterised by a lattice orientation g and by a volume fraction.The sum of the volume fractions of all crystallites that are within a range dg around an orientation g is equal to f(g) dg.f(g) is a statistical distribution function, called "orientation distribution function" (ODF), which characterises the crystallographic texture of the polycrystal (Bunge, 1982).The calculation then proceeds as follows for each individual crystallite: deformation is achieved by multiple slip.There are n slip systems available.The subscript s will be used to refer to a particular slip system.the slip rates 's must be chosen in such way, that (i) together they achieve the plastic part of the prescribed d D; (ii) the internally dissipated frictional work (per unit time and volume) is minimal: n ff x El'sl Min s=l c is the critical resolved shear stress.
(6) Note that some slip rates ', will be found equal to zero.This means that the corresponding slip systems are not active.These zero slip rates do not contribute to fi in eq. ( 6).Once the values of the slip rates are found, the local stress acting in the crystallite is found by using the requirement, that for active slip systems, the resolved shear stress must be equal to the critical resolved shear stress.Note that the local stresses are functions of g, and that they are usually different in each crystallite.
More detailed descriptions of this theory have for example been given by Kocks (1970) or by Aernoudt et al. (1993).Computer codes are available that can find the values of the ',, the ratio v/r and the local stress, as functions of the crystallite orientation g and the prescribed strain rate D (Van Houtte, 1988).It is then customary to define the Taylor factor as follows: M (7) TCDM DvM is the von Mises equivalent strain rate: DvM -D/.,D (8)   It is a scalar measure of the strain rate.
It is clear that M is a function of D and of g.Suppose now that after having calculated M, one doubles D and calculates M again, fi will have doubled, but so has D,M; this results in an unchanged value of M. Therefore M can be considered a function of the crystallite orientation g and the ratio D/DvM.The latter describes the strain mode.It is then useful to introduce a sort of direction in strain rate space denoted by a: a is a tensor with unit length.It is proportional to D/Dv, so M can be considered a function of a instead of D/D,.For the purpose of the present work, we only need the Taylor factor M, not the local stress acting in the crystallites.The Taylor theory will be used to construct a data base of M-factors: they will be calculated once and for all for a large number of crystal orientations g and for a large number of strain modes.
Calculation of the average Taylor factor of a polycrystal It follows from eq. ( 7) that the plastic work per unit time and volume is given by: fi, Dv M (g, a) (10) This can be averaged over the polycrystal, leading to with f(g) is the ODF that describes the crystallographic texture, as already explained above.
It can be obtained from X-ray or neutron diffraction measurements (pole figures), followed by an appropriate data processing (Bunge, 1982).Hence for or a particular sample with a known ODF, the average Taylor factors M can be obtained from the data base of Taylor factors M.This is done for a discrete set of strain modes characterised by unit tensors a (approx. 70 000for an angular resolution of 10 in strain rate space).The integrals in eq. ( 12) can be calculated very rapidly by a computer if series expansion methods are used.This was first shown by Bunge (1970).Note that eqs.(11-12) together offer a means of calculating the average plastic work dissipated in the polycrystalline material from the knowledge of the strain rate tensor, since Dv and a can be calculated easily from D. The texture acts as a parameter in this calculation, as well as the critical resolved shear stress of the slip systems.It will be shown in the next section how these formulas can be used as a basis for the calculation of a relationship between macroscopic stresses and strains.

Macroscopic Stress-Strain Relationship
Let us now treat a material as a continuum, disregarding its polycrystalline nature.
Stresses will now be "macroscopic", i.e. a sort of average in a volume which is large compared to the size of the crystallites.Assume that for such material I, the plastic work per unit time and per unit volume, can be obtained as a function of the plastic strain rate D. The stress should not appear as an argument or a parameter of this function.
P. VAN HOUT ET AL.
For materials which are plastically incompressible, it can then be shown that the deviatoric stress is given by (Van Houtte, 1994): For strain rate insensitive materials, v turns out to be equal to 1.The function I/ (D) can be regarded as one of a. dual set of plastic potentials (Hill, 1989, Van Houtte,   1994).Eq. ( 11) shows that W can indeed be considered as a function of D, since it essentially depends on DvM and a, both simple functions of D (eqs. (8-9)).It then becomes possible to calculate the deviatoric stress from the plastic strain rate.Work hardening can be incorporated by making " in eq. ( 11) a function of the equivalent strain.

Interpolation in strain rate space
There is however a serious practical problem.The expression for lcontains the texture- dependent function M (a).The latter is only known for a finite number of strain modes a belonging to a discrete set.An interpolation for intermediate strain modes is needed if eq. ( 13) must be evaluated.It must be such that at least the first order partial derivatives are continuous, otherwise the deviatoric stress cannot be obtained as a continuous function, which would almost certainly lead to unstable behaviour of the FE code.One way of turning the discrete set of data for M into a continuous function of a is by fitting an analytical function G(a) to it (Van Houtte et al., 1989; Van Bael, 1994).Let G for example be a 6th-order series expansion.For simplicity, a contracted index notation is used for the components of a (see e.g.Appendix 1 in Van Houtte et al., 1992): G(a) F pqrst apaqarasat + F pqrstu apaqarasatau (15) The series expansion consists of one term of odd rand and one of even rank.It follows from eq. ( 9) that the sum of the squares of ap is always equal to 1.This property can be used to show that the 5th order term of eq. ( 15) implicitly includes terms of order and 3, and the 6th order term terms of order 0, 2 and 4.There is no "stress differential effect" when all coefficients Fpqrst are zero, i.e. the sign of the stress tensor is merely reversed if the sign of the strain rate tensor is reversed.The coefficients Fpqrst__ and Fpqrstu must be obtained by a least-squares fitting of G to the discrete values of M that have been found by the methods explained in the previous section.This method has-been used for many applications (Van Bael, 1994; Van Bael et al., 1991a,  1991b, 1994; Winters et al., 1994).It is similar to the method proposed by Arminjon  and Bacroix (1991), in which the F-coefficients of eq. ( 15) are calculated from the C-coefficients of the texture function f(g) by means of an analytical formula instead of by fitting G(a) to M. The method was originally limited to the 2nd order, but has been extended later to the 4th order (Arminjon et al., 1994).
Unfortunately, there may be problems due to loss of convexity when the order of the series expansion of G(a) is higher than 2. The problem has been reported by Van Bael (1994) as well as by Arminjon el al. (1994).Figure 2 shows a r-plane section of the yield locus of a polycrystalline material with a rather sharp texture.The n-plane is a plane which makes equal angles with the axes representing the three normal stresses.It contains all stress states which are deviatoric, i.e. for which o-l + o'22 + o'33 0.Moreover, all shear stress components are zero for the section of Figure 2. The shear strain components of the plastic strain rate tensor associated to these stress states are then zero as well, at least when the plastic properties of the material feature orthorhombic symmetry, which is the case for the example shown in Figure 2. The yield curve shown by this figure has been obtained by applying eq. ( 13) in combination with a 6-th order series expansion (eq.( 15)) for a series of strain modes, each represented by a unit vector a (eq.( 9)).Since the shear strain rates are zero, these vectors can also be represented in the section of stress-strain rate space shown by Figure 2. Since they are unit vectors, the end points of these vectors would be on a unit circle.For each point of this "strain mode" circle, eq. ( 13) gives a stress state.Travelling along the circle then generates all points of the yield curve.It is seen that the latter features some near-vertices, which are points with a very high curvature.Magnification of these vertices reveal a non-convex behaviour (Figure 2).Apparently, the stress sometimes "moves back" at one travels further along the strain mode circle.The resulting curve has been called a "fishtail".The risk of encountering such fishtails increases with the order of the series expansion and with the sharpness of the texture.A short comment on the fundamental reason for the formation of such "fishtails" will be given at the end of the next section.Figure 2 x-plane sections of the yield locus of a fcc polycrystal.The ODF is a Gaussian distribution centered around the (100)[011] orientation (0 16-5-) "Fishtails" can be seen at the vertices.
_Spline functions offer another possibility for interpolation within the discrete set of M(a) data.Strain rate space can be reduced to five dimensions for incompressible materials.A "unit vector" such as a can then be described by four independent angles (Van Houtte, 1994).So the.__spline functions_ have to be 4-dimensional in nature.A system was tested in which not M, but 1/M was represented by 4-dimensional parabolic spline functions which were especially designed for this problem.Let the values of these functions be R(a) (sort of "radius" of a surface of constant ).R(a) itself then is continuous, and so are its partial derivatives of first order; not so its partial derivatives of second or higher order.It was found that this method can more easily represent regions with high curvature than the series expansion method, but that it does not guarantee convexity either when stresses are calculated by means of eq. ( 13).Although these two interpolation schemes have been designed for use in strain rate space, they can also be used in deviatoric stress space.The tensor a then represents a direction in deviatoric stress space or a "stress mode".It will be called b in the next section.A function R(b) could then represent the angular dependency of the radius of a yield locus in a certain direction.
Transformation from strain rate space to deviatoric stress space A surface in strain rate space for which W is constant is called an equipotential surface (Van Houtte, 1994).Eq. ( 11) can be used to construct an expression in polar co-ordinates for it: OvM-M--(a) (16) since a represents the direction in strain rate space, and Dvr is proportional to{Dij Dij, the length of the radius vector (eq.( 9)).This makes it possible to construct the equipotential surface point by point.The question then arises, whether a similar formula can be found for deviatoric stress space, allowing to calculate the length of the radius vector for each direction in stress space.For strain rate insensitive materials which are plastically incompressible, the equipotential surface in stress space would be the yield locus, and it would be independent of the value of W. Unfortunately, it turns out that in general, it is not possible to analytically convert an equation known in strain rate space such as eq. ( 16 into an equation in stress space.As far as the authors know, it can only be done for the von Mises yield locus or for Hill's anisotropic yield locus (Van Houtte, 1994).The problem can however be solved in an algorithmic way.Such algorithm would guess the strain mode (characterised by a) that corresponds to the plastic strain which corresponds to a certain direction in stress space characterised by a unit vector b.A first possibility to construct such algorithm could be based on eq. ( 13) in combination with eq. ( 11).A more convenient (but otherwise equivalent) algorithm is formulated as a minimisation problem (Van Houtte et al., 1989, Van Houtte,  1994) The minimum must be found by varying the strain mode a.The direction in stress space is prescribed by a (deviatoric) "stress mode" bij.The plastic stress itself would be given by: {ij sbiy (18) So s is a scalar that characterises the length of the radius vector in stress space.The minimisation algorithm can even be used for a discrete set of M(a) data without interpolation.This would be called the "geometrical method" (Van Houtte et al., 1989).Figures 3a-b show an example for which interpolations based on spline functions have been used.A Dl-D22 section of an equipotential surface in strain rate space is given by Figure 3a.The corresponding cr-o'22 section in stress space can be seen in Figure 3b.The latter figure has been generated using the minimisation algorithm described above.A radius vector in Figure 3a corresponds to a normal to the yield locus in Figure 3b.As a result, "vertices" (example: point E) in Figure 3a correspond to flat portions of the yield curve in Figure 3b (example: region between A and B) and vice-versa.
The worst problems of loss of convexity occur when the interpolation scheme used in strain rate space is be unable to represent an almost linear portion of the equipotential surface (e.g. between F and G in Figure 3a) without inducing some oscillations around a straight line.These oscillations then cause a fish-tail at a vertex in stress space.(The vertex C in Figure 3b corresponds to the section FG in Figure 3a, but there is no fish-tail in this case). . (a) Figure 3 Equipotential surfaces in strain rate space and in stress space (yield locus) of an A13004 alloy which has been cold rolled.The stresses and strain rates are expressed with respect to coordinate axes which make angles of 45 and 135 to the rolling direction.
(a) D:D22 section in strain rate space.D33 is such that D+D22+D33 0. The shear strain rates are zero.
Figure 3 Equipotential surfaces in strain rate space and in stress space (yield locus) of an A13004 alloy which has been cold rolled.The stresses and strain rates are expressed with respect to coordinate axes which make angles of 45 and 135 to the rolling direction.
(b) (11-(22 section of the yield locus in non-deviatoric stress space.(33 0 for this section, as well as the shear stresses.

IMPLEMENTATION OF PLASTIC ANISOTROPY IN THE FINITE ELEMENT METHOD
Introduction Superscripts are used to distinguish between elastic and plastic strains.The total strain is meant when no superscript is used.The finite element method is used since many years for the calculation of stress distributions in bodies that undergo elastic deformation.Such bodies are subdivided in elements with simple geometrical shapes.The elements have nodal points at comer points and sometimes at intermediate points on the ribs of the elements.The principles of the FEM method are the following: l) The continuous deformation field of an element is assumed to be solely determined by the displacements of the nodes.Interpolation formulas are used to calculate the displacements or the strains at other points of the element.2) All forces acting on an element (except body forces) are assumed to act at the nodal points.
3) Using Hooke's law, a constitutive relationship is established between the displacements of, and the forces at the nodal points for each element.4) These constitutive relationships of the individual elements are combined into an overall equation for the entire structure: [u], [F] are column matrices that contain the components of all the modal displacements and of all the external forces acting on the nodal points.
[K] is a square matrix: the "stiffness matrix", based on the geometry and on Hooke's law.The matrix equation above is solved for the nodal displacements.Once these are known, the distribution of the strains can be calculated in each element, and, by using Hooke's law, the distribution of the stresses as well.
For more background reading, we recommend Zienkiewics (1977) and Rowe et al.
(1990) and the references listed there.
The elastic FE methods have served as a basis for developing methods that can be used for the simulation of metalworking processes.One distinguishes mainly between elastic-plastic, rigid-plastic and visco-plastic methods (Rowe et al., 1990).We have decided to use the elastic-plastic method.
The Elastic-Plastic Finite Element Method The total strain to be applied to the body will now be divided in small strain increments.In each step of the simulation, one first has to check whether the strain increment is totally or partially in the elastic range.This of course requires the knowledge of the yield locus.If it is in the plastic range, then one tries to solve for the increments of the nodal displacements using an elasto-plastic stiffness matrix in eq. ( 19).This matrix will usually be different for each step of the simulation.It is based on the stiffness matrices of the individual elements.For the elastic case, these were based on Hooke's law: l/j Ci/!k, (0) which is now replaced by: Aij Cijkl AEkl (21) in which the increment of total strain is the sum of the increments of elastic strain and plastic strain: Ae =/te + ze,' The Ciu-coefficients are sometimes called the tangent moduli.It is not clear beforehand that an expression such as eq.( 21) can be established at all, since it implies that the stress increment can be calculated from the total strain increment without previous knowledge of the decomposition of the latter in an elastic and a plastic part.Zienkiewics et al. (1969) and Zienkiewics (1977) have shown that it is possible indeed.Note that the Ciju are different for each element.Anisotropic Tangent Moduli Derived from the Yield Locus The formulas will now be given for the calculation of anisotropic tangent moduli from the texture-based yield loci described in the previous sections.For reasons of simplicity, the discussion will be restricted to infinitesimal strain increments.The modulus C#u is then def'med as follows: with d Cju de u (23) Cijkt must be estimated at the beginning of each new increment.At that moment, the deviatoric stress t:r', and the plastic strain rate tiy Dij are known from the previous step.It has been shown by Van Houtte et al. (1992) and by Van Bael (1994) Ciji-(rs dC Cel ) ars + atu rstu de The ai represent the strain mode at the end of the previous increment (eq. ( 9)).Strain hardening is introduced through the derivative of the critical resolved shear stress with respect to the (macroscopic) von Mises equivalent strain.The evolution of is often approximated as a function of the average total slip of the polycrystal (Gil  Sevillano et al., 1980; Tom6 et al., 1984).The following approximate expression for F is used: It is furthermore assumed that a function (') which would be valid for any strain mode can be obtained from tensile test data or from torsion test data (Gil Sevillano  et al., 1980; Tom6 et al., 1984).The derivative of to be used in eq. ( 26) can be obtained from such function in the following way: Finite Increments dz dz g(a) (28) de,,r d' The instantaneous elasto-plastic tangent modulus as given by eq. ( 26) only serves as the basis for an algorithm that must model the material behaviour in an element during a finite strain increment.The objectives of such algorithm could for example be: at the end of the strain increment, a stress state must be found which is exactly on the yield locus.It is necessary to take strain hardening into account.the modulus linking the (finite) stress increment to the (finite) strain increment must be the average of the elasto-plastic modulus (for infinitesimal strain increments) for the trajectory followed on the surface of the yield locus during the stain increment (or an estimate of it).
An example of such algorithm is the mean-normal method (Rice and Tracey, 1973;  Tracey, 1976; Rowe et al., 1990).It has first been used for isotropic materials with avon Mises yield locus.An anisotropic variant has been developed by Van Bael (1994), who used an anisotropic yield locus based on eq. ( 15) (series expansion in strain rate space) and on eq. ( 17).Other algorithms have been proposed by Bacroix and Gilormini  (1995) (for sheet metal forming only) and by Winters (1995).

Application
Several anisotropic elasto-plastic FE-codes have been developed.Van Bael et al.  (1991a-b) have developed an anisotropic version of the 3-dimensional elastic-plastic FEM code "epfep 3" of the University of Birmingham.Texture based anisotropy was brought in through a plastic potential in strain rate space, represented as a 6th-order series expansion (eq.( 15)).Eq. ( 17) was used to explore the yield locus in stress space.
The program was tested by means of the NAFEMS single-element benchmark for triaxial displacements (Van Bael et al., 1991a).In another study, a tensile test on an anisotropic material was simulated.The test material was a deep drawing steel with an average r-value of 1.8.The FEM program indeed predicted the correct ratio of width strain to thickness strain (Van Bael et al., 1991b).Picksley et al. (1994) have used the program to perform a simulation of an upsetting test of cylindrical titanium samples with a non-axisymmetric texture (Figures 4-5).Finally the program has also been used for a study of the strain distribution in standard tensile test samples made of anisotropic materials.It could be demonstrated that the strain distribution is not homogeneous, and that it depends on the degree of anisotropy.It was also found that the heterogeneity of the strain distribution would cause systematic errors on the measured values of Lankford factors (r-values), and that these errors depend on the position and the types of extensometers used (Van Bael et al., 1994).
This program suffers from the repeated use of eq. ( 17), which is calculation intensive.In later work, it was decided to develop a faster method.Winters (1995) also starts from eqs. (11-12) to generate the plastic potential from the texture.He however does not perform interpolations in strain rate space.Instead, he generates the entire yield locus in deviatoric stress space by means of the geometrical method using eq.( 17).This leads to a discrete representation of the yield locus in stress space, in the form of a set of radius vector lengths for a number of "stress modes" b.The same interpolation techniques can be used in deviatoric stress space as in strain rate space.A 6-th order series expansion was used, in combination with an algorithm for modelling the material behaviour during a finite strain increment (see previous section).All this was implemented in a 3-dimensional version of the elastic-plastic finite element program Figure 4 FE mesh used for the upsetting of cylindrical titanium specimens.
"Lagamine" of the University of Liege.Some results have already been published: a study on the strain distribution in anis.otropicplane strain tensile test samples (Winters  et al., 1995b), and the results of a simulation of a cup drawing test performed on anisotropic steel sheet (Figures 6-7).When compared to experimental observations of the earing profile, the results are very satisfactory, as an earing profile could be simulated that could never have been found by using a Hill-type yield locus.However, there are some small differences between predicted and calculated earing profiles, which are ascribed to the fact, that the texture evolution has not been taken into account during the simulation (Munhoven et al., 1995).

DISCUSSION
A serious problem is the risk of losing convexity which is inherent to the interpolation methods (series expansion, spline functions) used to obtain continuous representation of the plastic potentials in strain rate space.This is the method used by the anisotropic version of the "epfep 3" program.With our present knowledge, it can only be avoided totally by using Hill-type expressions or second order series expansions; but then, the accuracy of reproducing anisotropic yield loci, r-values and the like is insufficient.It is felt by the present authors that a 6-th order series expansion is the minimum requirement for a reproduction of the yield loci of engineering materials, especially in the case of aluminium alloys and titanium alloys.Loss of convexity would manifests itself by the appearance of "fishtails" at the vertices of the yield loci. Very small fishtails have been observed for some aluminium sheets with sharp textures.Other engineering materials with less sharp texture did not present the problem.
No problem related to loss of convexity was encountered in a new method (implemented in Lagamine), in which the plastic potential is first transferred as a whole form strain rate space to stress space by means of the "geometrical method", after which a series expansion was used in stress space.The method also has another advantage: it requires less calculation time than the previous method, which used a series expansion in strain rate space.It must be admitted that Bacroix and Gilormini (1995) developed an algorithm that integrates the constitutive law during a finite strain increment that can work with a plastic potential in strain rate space and which does not require a lot of calculation time (for a 2-dimensional application).
The two anisotropic codes developed so far do not yet feature a direct coupling between the evolving texture (the texture changes during the plastic strain) and the yield locus.Strategies will have to be developed to that purpose.It is very well possible, that this will be easier for methods which use a plastic potential in strain rate space than for those which use a yield locus described in stress space.

CONCLUSIONS
It is possible to derive a complete yield locus from the ODF that describes the crystallographic texture of a single phase metal.A representation of these yield loci which can be implemented in a code for elastic-plastic finite element simulations of metalworking processes has been developed.However, it is also necessary to adapt the code itself, because the formulas that calculate the tangent moduli must be sufficiently general.Two different anisotropic FE-codes have been developed so far.They are both 3-dimensional, elastic-plastic codes.The results obtained by using these codes are satisfactory and reach beyond that which can be achieved by using a Hill- type yield locus.However, it is still necessary to design ways to overcome the problem of loss of convexity (a problem only for sharp textures).
At present, the anisotropic FE codes do not use a permanent coupling between the yield locus and the texture of the material as it evolves during the simulation itself.This may cause minor discrepancies between the observed and calculated earing profiles in a cup drawing test.The problem could be worse for processes which involve liigher strains than cup drawing.It must be admitted that the method proposed by Kalidindi and Anand (1993) and by Dawson et al. (1994) automatically takes texture evolution into account, since it does not represent the material behaviour by a plastic potential, but directly uses polycrystal plasticity inside the FE program.But this method requires very powerful computing facilities.