- Email: [email protected]
An effective split of flow and die deformation calculations of aluminium extrusion
Journal of Materials Processing Technology 88 (1999) 67 – 76
An effective split of flow and die deformation calculations of aluminium extrusion H.G. Mooi a,*, P.T.G. Koenis b, J. Hue´tink a a
Department of Mechanical Engineering, Uni6ersity of Twente, PO Box 217, 7500 AE Enschede, The Netherlands b BOAL B.V., PO Box 75, 2678 ZH De Lier, The Netherlands Received 11 October 1997
Abstract The purpose of the research project of which this paper is a result is to gain insight in design and control parameters of aluminium extrusion by FE simulations. This is done by means of a thermal FE code based on an ALE algorithm (Dieka). A short overview is given of the material models as well as the FE methods that are used. It is shown that it is a very efficient way to separate the aluminium flow and die deformation calculations. Furthermore, it has made it plausible that the complex three-dimensional flow through the extrusion die can be studied very well by means of two-dimensional plane strain calculations. A comparison is made between two- and three-dimensional calculation results. Finally some results of die deformation calculations are shown that are based on different deformation mechanisms. © 1999 Elsevier Science S.A. All rights reserved. Keywords: Finite elements; Material models; Rate dependency
1. Introduction Aluminium extrusion is a forming process in which hot aluminium is forced through a die with an opening that closely resembles the profile that has to be formed. The profiles can have an enormous variety of cross sections. Recently, the complexity of the profiles has increased because of the higher requirements that are set onto the profiles. According to the experience of die designers and those who correct the dies, very small changes in the shape of the die have large effects on the profiles. The increasing complexity of the profiles together with the sensitivity of the profile to the exact shape of the die require better insight into the design of the dies. Therefore, the Eureka project ‘Die Design Control’ (EU 554) project was started. It was initiated by the extrusion company BOAL B.V. (in De Lier, The Notation: in this paper, (second and fourth order) tensors are represented by bold-face type characters while scalars are indicated by normal characters. * Corresponding author. Current address: TNO Road-Vehicles Research Institute, PO Box 6033, 2600 JA Delft, The Netherlands. Fax: + 31-15-262-4321. E-mail address: [email protected] (H.G. Mooi)
Netherlands) and the die producer Phoenix (Paderno, Italy). One of the aims of the project was the investigation of the thermomechanical phenomena that determine the behaviour of the aluminium and the die during production. By means of these insights the design of the extrusion dies was to be improved. In this paper some of the results of this project are described. A more profound elaboration can be found in [1]. Firstly, the material models used are briefly described in Section 2. Subsequently the Finite Element model as well as some numerical details are explained in Section 3. In Sections 4 and 5, respectively some results of the calculations of the aluminium flow, and the die deformations are shown. A specific aim of this paper was to show that the results of the two-dimensional calculations can be compared quite well with their three-dimensional counterparts. Finally some conclusions are drawn in Section 6.
2. Material models For the modelling of aluminium extrusion two material models were implemented into the Finite Element
0924-0136/99/$ - see front matter © 1999 Elsevier Science S.A. All rights reserved. PII: S 0 9 2 4 - 0 1 3 6 ( 9 8 ) 0 0 3 8 0 - X
68
H.G. Mooi et al. / Journal of Materials Processing Technology 88 (1999) 67–76
code, Dieka. Firstly, a thermal rate-dependent model was built in for the description of aluminium at elevated temperatures. Secondly, a creep model was implemented capable of modeling time-dependent deformation of the die. Both models can be incorporated into a viscoplastic or an elasto-viscoplastic formulation, treated in Sections 2.1 and 2.2, respectively. The viscoplastic and elasto-viscoplastic formulations are founded on the theory of thermodynamics, for more details see [1,2]. Large deformations and thermal influences are taken into account. A new Helmholz free energy function was defined from which the isotropic elastic behaviour can be derived. Contact behaviour is described by means of contact elements based on a penalty formulation. The Coulomb friction law was applied.
2.1. Viscoplastic models The viscoplastic models are based on the split of the Cauchy stress tensor, s, into an isotropic and a deviatoric part: s = −p1 + s
(1)
in which p is the hydrostatic pressure, 1 is the second order unity tensor and s is the deviatoric stress. In Dieka, the pressure is updated by means of the following differential equation: p; = − Cbtr(D)+ a6CbT:
(2)
where Cb is the bulk modulus, tr (…) represents the trace operator, D is the rate of deformation tensor, a6 is the coefficient of thermal expansion, T is the (absolute) temperature and T: represents its material time derivative. An elastic formulation for the isotropic part is obtained, since the pressure is updated by means of an incremental procedure such as: p = p0 +Dp
(3)
with p0 being the ‘old’ pressure at the beginning of a calculation step and Dp the linearized pressure increment following from Eq. (2). Note that this approach is reminiscent of the widely applied penalty formulation for the plastic incompressibility, but it is quite different. For the approach in this paper, the pressure was calculated by the physical expression for the pressure rate, Eq. (2). In the case of a penalty formulation, the pressure is calculated by the unphysical expression: p= − C*tr(D) b
(4)
in which C*b represents the elastic penalty and T0 the initial reference temperature of the problem. It is stressed that this expression is NOT used in the underlying formulation. The deviatoric part of the stress split in Eq. (1) is determined by means of a simple linear viscous relation with a non-constant viscosity:
(5)
s=2h(T, D, h)Dd
in which h is the viscosity function that may depend on the temperature, the rate of deformation and the deformation history, gathered in the parameter h. Note that only the deviatoric part of the rate of deformation appears in this part.
2.1.1. Rate dependent model For low strain rates the material behaviour of aluminium can be well described by the curve-fit: s¯ = C(e¯ , T)e¯; n(T)
(6)
where s¯ = 3/2 s:s and e¯; = 2/3 D:D are the equivalent stress and strain rate. The material parameters C and n are temperature dependent, as indicated. The history or path dependent parameters are now only represented by the equivalent strain. If the definitions of the equivalent stress is worked out by means of Eq. (5), the following relation is obtained: s¯ =
'
3 s:s= 2h 2
'
2 D:D=3h e¯; 3
(7)
This is equal to the requirement that the scalar product of the stress and deformation rate tensors are equal to the product of their equivalent counterparts s:D =s¯ e¯; . Comparing Eqs. (6) and (7), an equation for the viscosity h is derived: h(e¯ , e¯; , T)=
C(e¯ , T) ; n(T) − 1 e¯ 3
(8)
Special cases of this model are: (i) rigid plasticity where n(T)=0 and the constant C equals the flow stress C= sy ;and (ii) constant Newtonian viscosity for which n(T)= 1. For aluminium, rather elaborate flow stress measurements are documented in a material atlas by Akeret [3]. The results in this atlas suggest, however, a softening effect upon high deformation (equivalent strain up to 700%; see curve b in Fig. 1). Opposed to Akeret, in the opinion of the present authors, softening is merely due to the fact that the experiments are not isothermal (as suggested). The decrease of the flow stresses is, therefore, principally caused by heating of the specimen by dissipation. Simple calculations confirmed this assumption. For these calculations a small piece of material was stretched adiabatically. The flow stress was assumed to be a function only of temperature s¯ = s¯ (T): hardening and softening were not taken into account. Two extreme cases were verified, an isothermal case in which all the heat is immediately led away (line a in Fig. 1) and an adiabatic case in which all heat is used to heat the specimen (line c in Fig. 1). Both, of course, do not resemble the experiments completely, because no special insulation was applied in the experiments. Never-
H.G. Mooi et al. / Journal of Materials Processing Technology 88 (1999) 67–76
theless, the experiments were performed so quickly (up to deformation rates of 10 s − 1) that not all of the heat could be transported away completely. According to the adiabatic calculation (line c) the flow stress decreased about 25% upon 700% plastic strain, whereas the flow stress measured experimentally by Akeret decreased only about 15%, see Fig. 1. Therefore, it might be concluded that it is possible that the softening described by Akeret is solely due to the heating of the test specimen. This confirmed the assumption that the flow stress in the experiments probably only decreased because of the temperature increase, which was calculated to be on the order of 50°C. Therefore, in this paper the viscosity in Eq. (8) is not taken to be dependent on the equivalent strain e¯ . This means that this material law is not history dependent. However, other authors suggest that the softening effect is not only due to the possible heating of the specimen, see Verlinden [4]. Nevertheless, he states that the exact reason for the softening is not known. Therefore, in this study the isothermal softening is neglected all the same. It should be stressed that, if necessary, history effects such as hardening or softening can be taken into account in the formulation that was used.
2.1.2. Creep Creep can be described very efficiently by using, for example, the Baily–Norton creep law that connects the equivalent creep strain to the equivalent stress and the absolute time t: e¯ c = Ae − Qc/RTs¯ m t n
(9)
where A, m and n are material parameters, Qc is the activation energy for creep, and R is the universal gas constant. The creep models can be incorporated into the finite element model in an analogous way by circumscribing Eq. (9) to a rate formulation and eliminating the time by means of the thus-obtained equations. If this result is compared with the general deviatoric part of the constitutive equation Eq. (5) the viscosity that has to be used for creep is acquired:
− m) (1 − n) e¯ c 1 e¯; (n c h= An n 3
1 m
This model can be built very easily as a viscous law but has some serious disadvantages. Firstly, because no limit stress is applied under which no creep occurs, creep deformation always occurs and, therefore, overestimates the real creep. Secondly, forces can hardly be prescribed with this model since the iteration process turned out to converge very badly in that case. Of course, this is not tolerable for creep calculations because it is obligatory to prescribe forces for creep calculations. Finally, this model needs an excessive number of iterations. A solution for these problems is to make use of elasto-viscoplastic models. Because of the fact that they are combined with an elastic fraction they behave more stably than the pure viscous models. In this paper, the empirical creep data is obtained from Berns [5].
2.2. Elasto-6iscoplastic models The starting point of elasto-viscoplastic models is Hookean-like rate equation: r; $ s = s+ L:De + fT: r
(11)
$
in which s represents the Jaumann rate of the Cauchy stress, r is the mass density, L is a material tensor, De represents the elastic part of the rate of deformation tensor and K is a thermal material tensor. This equation follows from thermodynamic constitutive principles, see [1,2]. Plasticity and viscoplastic effects are added to the formulation by assuming that the rate of deformation tensor can be additively split into an elastic, a plastic and a viscoplastic part: D= De + Dp + Dc
(12)
If a.o. the normality rule for plasticity is applied, the constitutive equation for isotropic Von Mises material follows: r; $ s = s+ Y:D + fT: − fc r
(10)
Generally, n ranges from 0 to 1 and m from 4 to 8. The latter power is numerically rather tedious because of its enormous non-linearity. Furthermore, in the first step the creep strain is zero, so that the strain-hardening formulation gives no results. This means that special measures again have to be taken for very small creep strains as well as in the first step, when oc is zero. In that case an estimation for the viscosity based on a formulation in which the time is not eliminated might be used.
69
Fig. 1. Material testing: experiment versus calculations.
(13)
H.G. Mooi et al. / Journal of Materials Processing Technology 88 (1999) 67–76
70
where: C = L− (1− h)Y fc = f=
3j Y:s 2s¯ 2
1 (G s − a6Cb 1− (1 − h)Y u G (T
(14)
in which E is the fourth order elasticity tensor, h a hardening parameter, Y a fourth order plastic yield tensor, fc is the representation of the creep terms, j is a dissipation function for creep, G the shear modulus and Y u is the thermal yield tensor. For small volume changes (r; negligible with respect to other terms) and for the case of small elastic deformations, the L tensor reduces to the well-known elasticity tensor, whereas the other material parameters change to: Y=
3G s s s¯ 2
Yu= h=
1 s¯ (G (sy − 3G G (T (T
Et Et + 3G
(15)
in which sy is the yield stress and Et is the plastic tangent modulus, i.e. the slope of the equivalent stress strain curve. The creep dissipation function, fc, is defined by: j =s:Dc = s¯ e¯; c
(16)
The equivalent creep strain rate is obtained from the rate form of Eq. (9), out of which the time is eliminated using Eq. (9) itself. Analogously, a dissipation function for rate dependency may be defined so that the rate dependent effects can also be incorporated into the elasto-viscoplastic formulation. At this point the necessary material models to describe aluminium extrusion are described.
3. Finite element method The investigation of the thermomechanical phenomena occurring during aluminium extrusion was performed by means of an Arbitrary Lagrangian Eulerian (ALE) finite-element code Dieka. This code has been developed during the last 15 years by Hue´tink et al. [6] and [7]. Other examples can be found by Hora [8] or Chenot [9]. An ALE formulation is a method with properties of both Lagrangian as well as Eulerian codes. The mesh may be completely bound to the material (Lagrangian) or the material may completely flow through the mesh (Eulerian). An arbitrary mixture between both methods may also be employed. This can be utilized in order to anticipate too much mesh distor-
tion. The ALE formulation is very suitable for the modelling of aluminium extrusion. On one hand, the Eulerian part can be applied to the flow character of the deformation of the hot aluminium within the die. On the other hand, the material has to be followed after it leaves the die which requires an ALE description. Furthermore, the (smaller) deformations of the die can be characterized best by means of the Lagrangian formulation. Dieka is based on a transient Galerkin formulation. The problems are discretized in time by means of finite time steps. The ALE formulation is implemented as a decoupled forward ALE, see also Baaijens [10]. This means that firstly the Lagrangian part of the step is performed. In this part the material increments are calculated by solving the weak form of the thermomechanical problem. Subsequently, the state variables are numerically integrated. Incrementally objective integration schemes as well as an implicit scheme for the integration of the creep strain are employed for this purpose. After the Lagrangian part, a new grid may be defined. As opposed to remeshing, with the ALE method the mesh topology must remain unchanged. The state variables of the problem have to be redistributed over this new grid. Because the new grid has the same topology this redistribution can be performed by means of convection schemes. These schemes are based on finite volume techniques. The best scheme turned out to be a limited second order k= 1/3-scheme with second-order time stepping, see Koren [11] and Stoker [12]. This means that the solution is second-order accurate in smooth regions, whereas near sharp layers as little oscillations as possible are introduced. The redistribution procedure is called the Eulerian part of the step. Although a true stationary state is not reached during aluminium extrusion, the stationary temperature distribution gives a good upper limit of the possible temperatures. If transient calculations were used for the modelling of aluminium extrusion, far too many calculation steps would be necessary to obtain the steady temperature values. For this reason, a stationary temperature solver was implemented into Dieka. For this solver convection–diffusion equations have to be solved. In order to stabilize the oscillatory behaviour of the Galerkin solution for convection–diffusion equations upwind techniques (SUPG and GLS by Hughes, see [13]) were also added. The stationary methods are tested on a typical shock problem. The use of the stationary solvers enhanced the efficiency of the calculations by about 30 times. Since the main objective of this paper is not the thermal aspects of aluminium extrusion, the reader is referred to e.g. see [1,14] for more details or results on this method.
H.G. Mooi et al. / Journal of Materials Processing Technology 88 (1999) 67–76
71
4. Two- and three-dimensional flow calculations Finally, some results of calculations of the extrusion of an aluminium tube are presented. In order to avoid too large calculations, the flow through the die and the deformation of the die were calculated separately. Both two- and three-dimensional flow calculations were performed. In this section the results of both approaches are compared. Tubes are extruded by means of a hollow die: the inner surface is formed by a core which is suspended by three or more legs in the die. So, the aluminium flow is split by legs, and after the legs the flow has to weld together again (in so-called welding chambers) see Fig. 2. Two finite element meshes were used for the two-dimensional plane strain calculations: a coarse mesh with 104 elements and a fine mesh with 462 elements. The elements were four node linear elements, with reduced integration for isotropic deformation. The coarse mesh, with the boundary conditions, is shown in Fig. 3. Typical dimensions of the mesh are given in millimetres. The depicted mesh is rather coarse but the results turned out to be quite reasonable. The dotted lines depict contact elements, which are used to describe friction. Of course, only one half is modelled because of symmetry. The boundary conditions were prescribed for the aluminium, according to Fig. 3 and Table 1. The rate dependent viscoplastic model was used (see Section 2.1), with the parameters according to Table 2. The working temperature was assumed to be constant at 450°C. This is a rather rough assumption, but these calculations were merely meant to verify the fact that
Fig. 3. Course mesh for the two-dimensional calculations.
the flow could be calculated by means of two-dimensional calculations. In the table the friction parameters are also given. Note that an upper limit for the frictional shear stress (tmax) is used. This resembles the maximum shear stress of aluminium for shear, so that a very thin shear layer is simulated at the wall. A high value of Poissons ratio was taken in order to obtain the mechanically stationary state quickly. This is allowed as Poissons ratio is only used to calculate the elastic compression (isotropic) which does not influence the deviatoric viscous behaviour. The final pressure is determined by the boundary conditions only. Recall that elastic deformation is only isotropic for the viscous model for aluminium, see Section 2.1. Thus, the calculated elastic compression is much too small, which was verified not to influence the final results. In order to verify the two-dimensional calculations, three-dimensional flow calculations were also performed. Because only the flow characteristics needed to be verified, it sufficed to consider the isothermal case. Because of the symmetry of the die, only one eighth needed to be modelled. Again a coarse and a fine mesh were used with eight-node bulk and eight-node contact elements (1319 and 4995 elements, respectively). A commercial pre-processor (Patran) was used for the meshing. The coarse mesh is shown in Fig. 4. For the fine mesh, all elements of the coarse mesh were divided into eight elements.
Table 1 Boundary conditions for 2D aluminium flow calculations G1 G2 G3 G4 G5 Fig. 2. Porthole extrusion with a 2D calculation plane (section A – A).
Inflow Die wall Die outlet Outlet Leg
Fx =0; 6y =8 mm/s 6x =0; 6y =0 Contact elements Fx =0; Fy =0 6x =0; 6y =0
H.G. Mooi et al. / Journal of Materials Processing Technology 88 (1999) 67–76
72
Table 2 Material parameters isothermal extrusion Aluminium
E (MPa) n C (MPa s−n) n
Friction
m tmax (MPa)
70 000 0.4999 33.08 0.15 0.10 10.0
The boundary conditions were analogous to the boundary conditions for the two-dimensional calculations (see Table 1). The entrance velocity was prescribed in the z-direction and at the largest part of the wall sticking was assumed. At the planes of symmetry corresponding boundary conditions were prescribed. Near the outlet, contact elements were used with the same properties as given in Table 1. The leg was modelled with a simple rectangular shape. The influence of this assumption was verified to be only a few percent. Note that these kinds of simplifications are important for the efficiency of three-dimensional meshing. The material parameters were taken to be identical to the values for the two-dimensional calculations, see Table 2. A P-GMRES iterative solver was used for these calculations. This was necessary because the calculation times of a direct solver for these kinds of problems are much too large. Since iterative solvers are very sensitive to the condition number of the matrix, Poissons ratio had to be chosen e.g. 0.47 instead of 0.4999. Subsequently, some results of both the two- and three-dimensional calculation are treated. First, in Figs. 5 and 6 the two- and three-dimensional pressure distributions are shown. The pressure is the main deformation load of the core of a porthole die, since the core has to bear the whole pressure. As can be seen, the pressures are very high which means that the core and the suspending legs are loaded quite heavily. This can be seen later on in the deformation calculations of the die, see Section 5.
Fig. 4. Course 3D mesh for the flow calculation.
Fig. 5. Pressure distribution for the 2D calculation (MPa).
The pressure distributions are comparable. Especially the pressure decrease and the linear isobars are striking. An overview of some results is given in Table 3. It can be seen that the differences between using the coarse and the fine grid are almost negligible for both two- and three-dimensional calculations. The differences in the three-dimensional case are due to the contact definitions in the outflow, which behave in an unstable manner. This means on the whole that mesh refinement does not greatly influence the flow calculations. Furthermore, it can be seen that the two-dimensional results approximate the three-dimensional results rather well. The pressure results agree to a sufficient degree of accuracy. Since the pressure on the core causes the largest part of the deformation of the die, this is an important result. It means that the two-dimensional calculations can be used to evaluate the loads on a die instead of three-dimensional calculations up to some degree. Care should be taken in drawing this conclusion too harshly. If, for example, a closer look is taken at the shear stress distribution, it has to be concluded immediately that this stress distribution can be only equivalent in the plane of symmetry of the 3D calculation, see Figs. 7 and 8.
Fig. 6. Pressure distribution (MPa; max 674 MPa).
H.G. Mooi et al. / Journal of Materials Processing Technology 88 (1999) 67–76
73
Table 3 Comparison of 2D and 3D results (course and fine; all results in MPa)
Course (2D) Fine (2D) Course (3D) Fine (3D)
Max. Von M.
Min. txy or txza
Max. txy or txzb
Max. p
76.3 82.4 59.2 77.8
−35.9 −39.3 −21.5 −29.9
21.8 21.2 21.7 22.7
444 440 589 674
a
Shear stress at the surface opposite to the leg in the flow direction. b Shear stress at the leg in the flow direction.
These values are taken for the comparison in Table 3. Since, therefore, the shear stresses on other planes are neglected, it could be expected that the pressures are erroneously estimated by the 2D calculations. This can be concluded from Table 3 also. It was also observed that the calculations were highly sensitive to the definition of the contact in the outflow region. Since this influences the results quite significantly, at the moment contact algorithms are under constructions that may be less sensitive. For further results on the 2D calculations the reader is referred to [1,15]. Of course, from the 3D calculations results can be obtained that can not be achieved by means of the 2D calculations. An example is the nonuniform outflow velocity at the outlet of the die. This results in differences in outflow velocity of the profile that is formed, see Fig. 9. In this figure the difference of the displacement of points on the profile at the outlet are shown. The coordinate s is taken along the outer surface of the tube in the tangential direction, see Fig. 6. The different curves represent the displacements of points in the outlet at five different moments in time. It can be seen that points behind the leg in the outlet (arc length s = 32) have smaller displacements than points where the aluminium flows directly along the legs (arc length s = 0). If the profile were not symmetric, these length
difference could cause curved and/or twisted profiles. Therefore, these length differences are interesting for extrusion practice. Generally, it can be stated that the two-dimensional calculations can be used for the estimation of the loads that an extrusion die is exposed to. Furthermore, it is not necessary to use refined meshes because the results do not significantly improve upon this refinement. For three-dimensional calculations there is a difference in calculation times of a few hours or a few days. The contact conditions are an important factor for the results. Therefore, close attention should be paid, for example, to oscillatory behaviour from time step to time step because this indicates irregular behaviour of the contact elements. Since extrusion is essentially a thermal process, thermal degrees of freedom should also be incorporated. Because two-dimensional calculations are shown to give reasonable results, it suffices to
Fig. 7. Shear stress distribution for the 2D calculation.
Fig. 9. Length differences in outflow.
Fig. 8. Shear stress distribution in 3D calculation (MPa).
74
H.G. Mooi et al. / Journal of Materials Processing Technology 88 (1999) 67–76 Table 4 Material parameters for the die
Fig. 10. Mesh for the three-dimensional deformation calculation of the die (upper and side view).
expand only the two-dimensional calculations with thermal degrees of freedom if the order of magnitude of the difference between two- and three-dimensional calculations are known. An example of this is treated in [1] and [15].
5. Three-dimensional deformation calculations In order to complete the consideration on the extrusion of a tube, the three-dimensional deformation calculations on the extrusion die are briefly treated also. The results of the (thermal) flow calculations were used to determine the loads for these calculations. Different deformation mechanisms of the die were considered: i.e. elastic, creep and ratcheting deformation. For more details on the ratcheting calculations, the reader is referred to [1]. Because of symmetry, only one eighth of the porthole part of the die has to be modelled. The three-dimensional mesh with the boundary conditions is shown in Fig. 10. Some details are left out of the model, such as the exact bearing length, the exact shape of the legs and most of the roundings. These simplifications do not have a great influence on the results. The mesh contains 628 nodes and 353 eight node brick elements. A refined model was also used for the same type of calculations with 3404 nodes and 2535 elements. The model was simply supported at the outer two rows of nodes on the lower side. Since these boundary conditions (rigidly suppressed) were too stiff compared to the elastic support of the surrounding tooling, e.g. the press, a smaller part of the die was supported than in practice. At the symmetry planes, symmetry conditions were prescribed (upper view in Fig. 10). The pressure in front of the die was calculated by means of the two-dimensional calculations with temperature dependent material parameters like in Section 4. This pressure was prescribed as a distributed load on the core of the die (q1 in Fig. 10). The dis-
Elastic
E (MPa) n sy,0(MPa)
176 000 0.3 850
Creep
C (MPa−ms−n) m n
3.2×10−26 8.0 0.28
tributed load on the leg consisted of two parts. Firstly, the pressure difference over the leg contributed to this load. These pressures also follow from the two-dimensional flow calculations. Note that on the downstream part of the leg a pressure also works in the direction opposite to the extrusion direction. Secondly, the friction forces on the side of the leg were summed and modelled as a distributed load on the front of the leg. Both loads on the leg were combined in the distributed load q2. The distributed loads q1 and q2 were taken to be 230 and 350 MPa, respectively. The material models used, were purely elastic, elastoviscoplastic and viscous creep models. The latter two were treated in Section 2. The material parameters are taken according to Table 4. A first calculation is used to check whether the deformations are immediately plastic by means of the Von Mises stress distribution, see Figs. 11 and 12. For this calculation an elasto-plastic model only is used. The yield strength at 500°C for this material is 850 MPa, see Table 3. Therefore, no direct plastic strains are present for these loads. However, it is clear that if thermal loads are added or the loads are increased only slightly, plastic deformation might be possible, see e.g. [1,15]. The main conclusion of this first calculation is that the porthole part of the die does not directly deform plastically. In practice permanent deformation of the die is, in fact, observed. Therefore, it is also interesting to investigate the creep deformation of the die under the earlier mentioned loads. In order to do so, the same
Fig. 11. Von Mises stresses (in MPa; max 810 MPa).
H.G. Mooi et al. / Journal of Materials Processing Technology 88 (1999) 67–76
75
compared. Peaks in the extrusion load might result in even higher creep deformation. The viscous and elastoviscoplastic description of creep gave results in the same order of magnitude. The latter should be favoured because of its higher accuracy and the fact that loads can always be prescribed. The viscous model, though, can be used if a quick indication of the creep deformation is needed. In that case, the first step with (prescribed displacements) should be made sufficiently small. Both models tend to converge very badly, which might be cured by using multi-phase models as treated in [6]. Fig. 12. Deformed three-dimensional mesh (exaggerated).
problem was calculated with a material model for creep (see Section 2). The material parameters are given in Table 4. Both elasto-viscoplastic and viscous creep formulations were applied. The displacement in time of the point on the centre line on the upper surface of the die for the elasto-viscoplastic and the viscous models is given in Fig. 13. Note that the upper line also incorporates an elastic part. The viscous and the elasto-viscoplastic model give somewhat different results. This can mainly be explained by the fact that the creep strain for the viscous formulation is incorrect in the first step. This is due to the fact that displacements needed to be prescribed instead of loads for the sake of convergence. Note the steepness of the elasto-viscoplastic line due to the highly non-linear behaviour of the creep law; this requires very small time steps. The displacements by creep deformation are still about an order smaller than the elastic deformations that measured a few tenths of a millimetre, see also Fig. 13. It can thus be concluded that creep plays at least some role in the permanent deformation of the dies. In practice, permanent deformation of porthole dies is measured as a few millimetres, but the history of these dies over 180 h was unknown so that it could not be
6. Conclusions The ALE method, as developed by Hue´tink, is very applicable to the modelling of complex aluminium extrusion. It was shown that interesting insights could be gained in the aluminium flow and die deformation aspects of aluminium extrusion by separating both calculations. The two-dimensional flow calculations were shown to give very reasonable estimates of the real three-dimensional stress state and deformations. Therefore, it was tolerable to use the two-dimensional (thermal!) calculations to determine the pressures and shear forces that are exerted on the die during extrusion. The deformation of the die was shown to be mainly elastic; the permanent deformations could be ascribed to creep. Future research should focus on more complex geometries, for which three-dimensional possibilities of present FE codes should be enhanced anyway. Also, the friction conditions might need closer attention.
Acknowledgements The research project was carried out within the framework of the Eureka project no. EU 554, which was initiated by BOAL B.V. in De Lier, The Netherlands, and Phoenix s.p.a. in Paderno, Italy. The support is gratefully acknowledged.
References
Fig. 13. Displacement of the centre line of the die.
[1] H.G. Mooi, Finite Elements Simulations of aluminium extrusion, Ph.D. thesis, University of Twente, The Netherlands, 1996. [2] H.G. Mooi, J. Hue´tink, Thermodynamically based material models for large deformations with application to Aluminium Extrusion, J. Mech. Phys. Solids (submitted). [3] R. Akeret, H. Jung, G. Scharf, Atlas of hot working properties of nonferrous metals, vol. 1, Deutsche Gesellschaft fu¨r Metallkunde (DGM) 1978. [4] B. Verlinden, A. Suhadi, L. Delaey, A generalized constitutive equation for an AA6060 Aluminium Alloy, vol. 28, Scripta Metallurgica , 1993, pp. 1441 – 1446.
76
H.G. Mooi et al. / Journal of Materials Processing Technology 88 (1999) 67–76
[5] H. Berns, F. Pschenitzka, Das Kriechverhalten von Warmarbeitssta¨hlen mit 5% Chrom (in German), Zeitschrift Werkstofftechn., no. (11), Verlag Chemie, Germany, 1980, pp. 258 – 264. [6] J. Hue´tink, On the simulation of thermomechanical forming processes, Ph.D. thesis, University of Twente, The Netherlands, 1986. [7] J. Hue´tink, Analysis of metal forming proc. based on a comb. Euler-Lagrangian FE formulation, In: J.F.T. Pittman, R.D. Wood, J.M. Alexander, O.C. Zienkiewicz (Eds.), Num. meth. in indust. forming processes, Pineridge Press, Swansea, UK, 1982, pp. 501 – 509. [8] P. Hora, L. Tong, J. Reissner, The FE simulation of 3-D bulk forming processes using the ALE formulation, Proc. of Metal Forming Process Simulation in Industry, Kro¨plin et al. (Eds.), Baden-Baden, Germany, Werner Fastenrath Druckerei, Wuppertal, 1994, pp. 46 – 59. [9] J.-L. Chenot, M. Bellet, The ALE method for the numerical simulation of material forming processes, in: Simulation of Materials Processing: Theory, Methods and Applications, Proc. Numiform ‘95, Balkema, Rotterdam, 1994, pp. 39–48.
.
[10] F.P.T. Baaijens, An U-ALE formulation of 3-D unsteady viscoelastic flow, Int. J. Numer. Methods Eng. 36 (1993) 1115–1143. [11] B. Koren, A robust upwind discretization method for Advection, Diffusion and Source terms, in: Notes on Numerical Fluid Mechanics, 45 (1993) 117 – 138, Vieweg, Braunschweig. [12] H.C. Stoker, J. Hue´tink, The arbitrary Lagrangian-Eulerian method in large deformation problems, in: Proc. EccoMass ‘96, Wiley, Paris, 1996. [13] T.J.R. Hughes, L.P. Franca, G.M. Hulbert, A new finite element formulation for computational fluid dynamics: VIII. the Galerkin/ Least-Squares methods for advective-diffusion equations, Comp. Methods. Appl. Mech. Eng. 73 (1989) 173 – 189. [14] H.G. Mooi, J. Hue´tink, Simulation of complex aluminium extrusion using an arbitrary Eulerian Lagrangian formulation, in: Simulation of Materials Processing: Theory, Methods and Applications, Proc. Numiform ‘95, Balkema, Rotterdam, 1995, pp. 869 – 874. [15] H.G. Mooi, A.J. Den Bakker, K.E. Nilsen J. Hue´tink, Simulation of Aluminium Extrusion by means of the Finite Element Method (FEM), in: Proc. 6th Int. Al. Extr. Technol. Seminar, Chicago, 1996, pp. 67 – 73.
An effective split of flow and die deformation calculations of aluminium extrusion H.G. Mooi a,*, P.T.G. Koenis b, J. Hue´tink a a
Department of Mechanical Engineering, Uni6ersity of Twente, PO Box 217, 7500 AE Enschede, The Netherlands b BOAL B.V., PO Box 75, 2678 ZH De Lier, The Netherlands Received 11 October 1997
Abstract The purpose of the research project of which this paper is a result is to gain insight in design and control parameters of aluminium extrusion by FE simulations. This is done by means of a thermal FE code based on an ALE algorithm (Dieka). A short overview is given of the material models as well as the FE methods that are used. It is shown that it is a very efficient way to separate the aluminium flow and die deformation calculations. Furthermore, it has made it plausible that the complex three-dimensional flow through the extrusion die can be studied very well by means of two-dimensional plane strain calculations. A comparison is made between two- and three-dimensional calculation results. Finally some results of die deformation calculations are shown that are based on different deformation mechanisms. © 1999 Elsevier Science S.A. All rights reserved. Keywords: Finite elements; Material models; Rate dependency
1. Introduction Aluminium extrusion is a forming process in which hot aluminium is forced through a die with an opening that closely resembles the profile that has to be formed. The profiles can have an enormous variety of cross sections. Recently, the complexity of the profiles has increased because of the higher requirements that are set onto the profiles. According to the experience of die designers and those who correct the dies, very small changes in the shape of the die have large effects on the profiles. The increasing complexity of the profiles together with the sensitivity of the profile to the exact shape of the die require better insight into the design of the dies. Therefore, the Eureka project ‘Die Design Control’ (EU 554) project was started. It was initiated by the extrusion company BOAL B.V. (in De Lier, The Notation: in this paper, (second and fourth order) tensors are represented by bold-face type characters while scalars are indicated by normal characters. * Corresponding author. Current address: TNO Road-Vehicles Research Institute, PO Box 6033, 2600 JA Delft, The Netherlands. Fax: + 31-15-262-4321. E-mail address: [email protected] (H.G. Mooi)
Netherlands) and the die producer Phoenix (Paderno, Italy). One of the aims of the project was the investigation of the thermomechanical phenomena that determine the behaviour of the aluminium and the die during production. By means of these insights the design of the extrusion dies was to be improved. In this paper some of the results of this project are described. A more profound elaboration can be found in [1]. Firstly, the material models used are briefly described in Section 2. Subsequently the Finite Element model as well as some numerical details are explained in Section 3. In Sections 4 and 5, respectively some results of the calculations of the aluminium flow, and the die deformations are shown. A specific aim of this paper was to show that the results of the two-dimensional calculations can be compared quite well with their three-dimensional counterparts. Finally some conclusions are drawn in Section 6.
2. Material models For the modelling of aluminium extrusion two material models were implemented into the Finite Element
0924-0136/99/$ - see front matter © 1999 Elsevier Science S.A. All rights reserved. PII: S 0 9 2 4 - 0 1 3 6 ( 9 8 ) 0 0 3 8 0 - X
68
H.G. Mooi et al. / Journal of Materials Processing Technology 88 (1999) 67–76
code, Dieka. Firstly, a thermal rate-dependent model was built in for the description of aluminium at elevated temperatures. Secondly, a creep model was implemented capable of modeling time-dependent deformation of the die. Both models can be incorporated into a viscoplastic or an elasto-viscoplastic formulation, treated in Sections 2.1 and 2.2, respectively. The viscoplastic and elasto-viscoplastic formulations are founded on the theory of thermodynamics, for more details see [1,2]. Large deformations and thermal influences are taken into account. A new Helmholz free energy function was defined from which the isotropic elastic behaviour can be derived. Contact behaviour is described by means of contact elements based on a penalty formulation. The Coulomb friction law was applied.
2.1. Viscoplastic models The viscoplastic models are based on the split of the Cauchy stress tensor, s, into an isotropic and a deviatoric part: s = −p1 + s
(1)
in which p is the hydrostatic pressure, 1 is the second order unity tensor and s is the deviatoric stress. In Dieka, the pressure is updated by means of the following differential equation: p; = − Cbtr(D)+ a6CbT:
(2)
where Cb is the bulk modulus, tr (…) represents the trace operator, D is the rate of deformation tensor, a6 is the coefficient of thermal expansion, T is the (absolute) temperature and T: represents its material time derivative. An elastic formulation for the isotropic part is obtained, since the pressure is updated by means of an incremental procedure such as: p = p0 +Dp
(3)
with p0 being the ‘old’ pressure at the beginning of a calculation step and Dp the linearized pressure increment following from Eq. (2). Note that this approach is reminiscent of the widely applied penalty formulation for the plastic incompressibility, but it is quite different. For the approach in this paper, the pressure was calculated by the physical expression for the pressure rate, Eq. (2). In the case of a penalty formulation, the pressure is calculated by the unphysical expression: p= − C*tr(D) b
(4)
in which C*b represents the elastic penalty and T0 the initial reference temperature of the problem. It is stressed that this expression is NOT used in the underlying formulation. The deviatoric part of the stress split in Eq. (1) is determined by means of a simple linear viscous relation with a non-constant viscosity:
(5)
s=2h(T, D, h)Dd
in which h is the viscosity function that may depend on the temperature, the rate of deformation and the deformation history, gathered in the parameter h. Note that only the deviatoric part of the rate of deformation appears in this part.
2.1.1. Rate dependent model For low strain rates the material behaviour of aluminium can be well described by the curve-fit: s¯ = C(e¯ , T)e¯; n(T)
(6)
where s¯ = 3/2 s:s and e¯; = 2/3 D:D are the equivalent stress and strain rate. The material parameters C and n are temperature dependent, as indicated. The history or path dependent parameters are now only represented by the equivalent strain. If the definitions of the equivalent stress is worked out by means of Eq. (5), the following relation is obtained: s¯ =
'
3 s:s= 2h 2
'
2 D:D=3h e¯; 3
(7)
This is equal to the requirement that the scalar product of the stress and deformation rate tensors are equal to the product of their equivalent counterparts s:D =s¯ e¯; . Comparing Eqs. (6) and (7), an equation for the viscosity h is derived: h(e¯ , e¯; , T)=
C(e¯ , T) ; n(T) − 1 e¯ 3
(8)
Special cases of this model are: (i) rigid plasticity where n(T)=0 and the constant C equals the flow stress C= sy ;and (ii) constant Newtonian viscosity for which n(T)= 1. For aluminium, rather elaborate flow stress measurements are documented in a material atlas by Akeret [3]. The results in this atlas suggest, however, a softening effect upon high deformation (equivalent strain up to 700%; see curve b in Fig. 1). Opposed to Akeret, in the opinion of the present authors, softening is merely due to the fact that the experiments are not isothermal (as suggested). The decrease of the flow stresses is, therefore, principally caused by heating of the specimen by dissipation. Simple calculations confirmed this assumption. For these calculations a small piece of material was stretched adiabatically. The flow stress was assumed to be a function only of temperature s¯ = s¯ (T): hardening and softening were not taken into account. Two extreme cases were verified, an isothermal case in which all the heat is immediately led away (line a in Fig. 1) and an adiabatic case in which all heat is used to heat the specimen (line c in Fig. 1). Both, of course, do not resemble the experiments completely, because no special insulation was applied in the experiments. Never-
H.G. Mooi et al. / Journal of Materials Processing Technology 88 (1999) 67–76
theless, the experiments were performed so quickly (up to deformation rates of 10 s − 1) that not all of the heat could be transported away completely. According to the adiabatic calculation (line c) the flow stress decreased about 25% upon 700% plastic strain, whereas the flow stress measured experimentally by Akeret decreased only about 15%, see Fig. 1. Therefore, it might be concluded that it is possible that the softening described by Akeret is solely due to the heating of the test specimen. This confirmed the assumption that the flow stress in the experiments probably only decreased because of the temperature increase, which was calculated to be on the order of 50°C. Therefore, in this paper the viscosity in Eq. (8) is not taken to be dependent on the equivalent strain e¯ . This means that this material law is not history dependent. However, other authors suggest that the softening effect is not only due to the possible heating of the specimen, see Verlinden [4]. Nevertheless, he states that the exact reason for the softening is not known. Therefore, in this study the isothermal softening is neglected all the same. It should be stressed that, if necessary, history effects such as hardening or softening can be taken into account in the formulation that was used.
2.1.2. Creep Creep can be described very efficiently by using, for example, the Baily–Norton creep law that connects the equivalent creep strain to the equivalent stress and the absolute time t: e¯ c = Ae − Qc/RTs¯ m t n
(9)
where A, m and n are material parameters, Qc is the activation energy for creep, and R is the universal gas constant. The creep models can be incorporated into the finite element model in an analogous way by circumscribing Eq. (9) to a rate formulation and eliminating the time by means of the thus-obtained equations. If this result is compared with the general deviatoric part of the constitutive equation Eq. (5) the viscosity that has to be used for creep is acquired:
− m) (1 − n) e¯ c 1 e¯; (n c h= An n 3
1 m
This model can be built very easily as a viscous law but has some serious disadvantages. Firstly, because no limit stress is applied under which no creep occurs, creep deformation always occurs and, therefore, overestimates the real creep. Secondly, forces can hardly be prescribed with this model since the iteration process turned out to converge very badly in that case. Of course, this is not tolerable for creep calculations because it is obligatory to prescribe forces for creep calculations. Finally, this model needs an excessive number of iterations. A solution for these problems is to make use of elasto-viscoplastic models. Because of the fact that they are combined with an elastic fraction they behave more stably than the pure viscous models. In this paper, the empirical creep data is obtained from Berns [5].
2.2. Elasto-6iscoplastic models The starting point of elasto-viscoplastic models is Hookean-like rate equation: r; $ s = s+ L:De + fT: r
(11)
$
in which s represents the Jaumann rate of the Cauchy stress, r is the mass density, L is a material tensor, De represents the elastic part of the rate of deformation tensor and K is a thermal material tensor. This equation follows from thermodynamic constitutive principles, see [1,2]. Plasticity and viscoplastic effects are added to the formulation by assuming that the rate of deformation tensor can be additively split into an elastic, a plastic and a viscoplastic part: D= De + Dp + Dc
(12)
If a.o. the normality rule for plasticity is applied, the constitutive equation for isotropic Von Mises material follows: r; $ s = s+ Y:D + fT: − fc r
(10)
Generally, n ranges from 0 to 1 and m from 4 to 8. The latter power is numerically rather tedious because of its enormous non-linearity. Furthermore, in the first step the creep strain is zero, so that the strain-hardening formulation gives no results. This means that special measures again have to be taken for very small creep strains as well as in the first step, when oc is zero. In that case an estimation for the viscosity based on a formulation in which the time is not eliminated might be used.
69
Fig. 1. Material testing: experiment versus calculations.
(13)
H.G. Mooi et al. / Journal of Materials Processing Technology 88 (1999) 67–76
70
where: C = L− (1− h)Y fc = f=
3j Y:s 2s¯ 2
1 (G s − a6Cb 1− (1 − h)Y u G (T
(14)
in which E is the fourth order elasticity tensor, h a hardening parameter, Y a fourth order plastic yield tensor, fc is the representation of the creep terms, j is a dissipation function for creep, G the shear modulus and Y u is the thermal yield tensor. For small volume changes (r; negligible with respect to other terms) and for the case of small elastic deformations, the L tensor reduces to the well-known elasticity tensor, whereas the other material parameters change to: Y=
3G s s s¯ 2
Yu= h=
1 s¯ (G (sy − 3G G (T (T
Et Et + 3G
(15)
in which sy is the yield stress and Et is the plastic tangent modulus, i.e. the slope of the equivalent stress strain curve. The creep dissipation function, fc, is defined by: j =s:Dc = s¯ e¯; c
(16)
The equivalent creep strain rate is obtained from the rate form of Eq. (9), out of which the time is eliminated using Eq. (9) itself. Analogously, a dissipation function for rate dependency may be defined so that the rate dependent effects can also be incorporated into the elasto-viscoplastic formulation. At this point the necessary material models to describe aluminium extrusion are described.
3. Finite element method The investigation of the thermomechanical phenomena occurring during aluminium extrusion was performed by means of an Arbitrary Lagrangian Eulerian (ALE) finite-element code Dieka. This code has been developed during the last 15 years by Hue´tink et al. [6] and [7]. Other examples can be found by Hora [8] or Chenot [9]. An ALE formulation is a method with properties of both Lagrangian as well as Eulerian codes. The mesh may be completely bound to the material (Lagrangian) or the material may completely flow through the mesh (Eulerian). An arbitrary mixture between both methods may also be employed. This can be utilized in order to anticipate too much mesh distor-
tion. The ALE formulation is very suitable for the modelling of aluminium extrusion. On one hand, the Eulerian part can be applied to the flow character of the deformation of the hot aluminium within the die. On the other hand, the material has to be followed after it leaves the die which requires an ALE description. Furthermore, the (smaller) deformations of the die can be characterized best by means of the Lagrangian formulation. Dieka is based on a transient Galerkin formulation. The problems are discretized in time by means of finite time steps. The ALE formulation is implemented as a decoupled forward ALE, see also Baaijens [10]. This means that firstly the Lagrangian part of the step is performed. In this part the material increments are calculated by solving the weak form of the thermomechanical problem. Subsequently, the state variables are numerically integrated. Incrementally objective integration schemes as well as an implicit scheme for the integration of the creep strain are employed for this purpose. After the Lagrangian part, a new grid may be defined. As opposed to remeshing, with the ALE method the mesh topology must remain unchanged. The state variables of the problem have to be redistributed over this new grid. Because the new grid has the same topology this redistribution can be performed by means of convection schemes. These schemes are based on finite volume techniques. The best scheme turned out to be a limited second order k= 1/3-scheme with second-order time stepping, see Koren [11] and Stoker [12]. This means that the solution is second-order accurate in smooth regions, whereas near sharp layers as little oscillations as possible are introduced. The redistribution procedure is called the Eulerian part of the step. Although a true stationary state is not reached during aluminium extrusion, the stationary temperature distribution gives a good upper limit of the possible temperatures. If transient calculations were used for the modelling of aluminium extrusion, far too many calculation steps would be necessary to obtain the steady temperature values. For this reason, a stationary temperature solver was implemented into Dieka. For this solver convection–diffusion equations have to be solved. In order to stabilize the oscillatory behaviour of the Galerkin solution for convection–diffusion equations upwind techniques (SUPG and GLS by Hughes, see [13]) were also added. The stationary methods are tested on a typical shock problem. The use of the stationary solvers enhanced the efficiency of the calculations by about 30 times. Since the main objective of this paper is not the thermal aspects of aluminium extrusion, the reader is referred to e.g. see [1,14] for more details or results on this method.
H.G. Mooi et al. / Journal of Materials Processing Technology 88 (1999) 67–76
71
4. Two- and three-dimensional flow calculations Finally, some results of calculations of the extrusion of an aluminium tube are presented. In order to avoid too large calculations, the flow through the die and the deformation of the die were calculated separately. Both two- and three-dimensional flow calculations were performed. In this section the results of both approaches are compared. Tubes are extruded by means of a hollow die: the inner surface is formed by a core which is suspended by three or more legs in the die. So, the aluminium flow is split by legs, and after the legs the flow has to weld together again (in so-called welding chambers) see Fig. 2. Two finite element meshes were used for the two-dimensional plane strain calculations: a coarse mesh with 104 elements and a fine mesh with 462 elements. The elements were four node linear elements, with reduced integration for isotropic deformation. The coarse mesh, with the boundary conditions, is shown in Fig. 3. Typical dimensions of the mesh are given in millimetres. The depicted mesh is rather coarse but the results turned out to be quite reasonable. The dotted lines depict contact elements, which are used to describe friction. Of course, only one half is modelled because of symmetry. The boundary conditions were prescribed for the aluminium, according to Fig. 3 and Table 1. The rate dependent viscoplastic model was used (see Section 2.1), with the parameters according to Table 2. The working temperature was assumed to be constant at 450°C. This is a rather rough assumption, but these calculations were merely meant to verify the fact that
Fig. 3. Course mesh for the two-dimensional calculations.
the flow could be calculated by means of two-dimensional calculations. In the table the friction parameters are also given. Note that an upper limit for the frictional shear stress (tmax) is used. This resembles the maximum shear stress of aluminium for shear, so that a very thin shear layer is simulated at the wall. A high value of Poissons ratio was taken in order to obtain the mechanically stationary state quickly. This is allowed as Poissons ratio is only used to calculate the elastic compression (isotropic) which does not influence the deviatoric viscous behaviour. The final pressure is determined by the boundary conditions only. Recall that elastic deformation is only isotropic for the viscous model for aluminium, see Section 2.1. Thus, the calculated elastic compression is much too small, which was verified not to influence the final results. In order to verify the two-dimensional calculations, three-dimensional flow calculations were also performed. Because only the flow characteristics needed to be verified, it sufficed to consider the isothermal case. Because of the symmetry of the die, only one eighth needed to be modelled. Again a coarse and a fine mesh were used with eight-node bulk and eight-node contact elements (1319 and 4995 elements, respectively). A commercial pre-processor (Patran) was used for the meshing. The coarse mesh is shown in Fig. 4. For the fine mesh, all elements of the coarse mesh were divided into eight elements.
Table 1 Boundary conditions for 2D aluminium flow calculations G1 G2 G3 G4 G5 Fig. 2. Porthole extrusion with a 2D calculation plane (section A – A).
Inflow Die wall Die outlet Outlet Leg
Fx =0; 6y =8 mm/s 6x =0; 6y =0 Contact elements Fx =0; Fy =0 6x =0; 6y =0
H.G. Mooi et al. / Journal of Materials Processing Technology 88 (1999) 67–76
72
Table 2 Material parameters isothermal extrusion Aluminium
E (MPa) n C (MPa s−n) n
Friction
m tmax (MPa)
70 000 0.4999 33.08 0.15 0.10 10.0
The boundary conditions were analogous to the boundary conditions for the two-dimensional calculations (see Table 1). The entrance velocity was prescribed in the z-direction and at the largest part of the wall sticking was assumed. At the planes of symmetry corresponding boundary conditions were prescribed. Near the outlet, contact elements were used with the same properties as given in Table 1. The leg was modelled with a simple rectangular shape. The influence of this assumption was verified to be only a few percent. Note that these kinds of simplifications are important for the efficiency of three-dimensional meshing. The material parameters were taken to be identical to the values for the two-dimensional calculations, see Table 2. A P-GMRES iterative solver was used for these calculations. This was necessary because the calculation times of a direct solver for these kinds of problems are much too large. Since iterative solvers are very sensitive to the condition number of the matrix, Poissons ratio had to be chosen e.g. 0.47 instead of 0.4999. Subsequently, some results of both the two- and three-dimensional calculation are treated. First, in Figs. 5 and 6 the two- and three-dimensional pressure distributions are shown. The pressure is the main deformation load of the core of a porthole die, since the core has to bear the whole pressure. As can be seen, the pressures are very high which means that the core and the suspending legs are loaded quite heavily. This can be seen later on in the deformation calculations of the die, see Section 5.
Fig. 4. Course 3D mesh for the flow calculation.
Fig. 5. Pressure distribution for the 2D calculation (MPa).
The pressure distributions are comparable. Especially the pressure decrease and the linear isobars are striking. An overview of some results is given in Table 3. It can be seen that the differences between using the coarse and the fine grid are almost negligible for both two- and three-dimensional calculations. The differences in the three-dimensional case are due to the contact definitions in the outflow, which behave in an unstable manner. This means on the whole that mesh refinement does not greatly influence the flow calculations. Furthermore, it can be seen that the two-dimensional results approximate the three-dimensional results rather well. The pressure results agree to a sufficient degree of accuracy. Since the pressure on the core causes the largest part of the deformation of the die, this is an important result. It means that the two-dimensional calculations can be used to evaluate the loads on a die instead of three-dimensional calculations up to some degree. Care should be taken in drawing this conclusion too harshly. If, for example, a closer look is taken at the shear stress distribution, it has to be concluded immediately that this stress distribution can be only equivalent in the plane of symmetry of the 3D calculation, see Figs. 7 and 8.
Fig. 6. Pressure distribution (MPa; max 674 MPa).
H.G. Mooi et al. / Journal of Materials Processing Technology 88 (1999) 67–76
73
Table 3 Comparison of 2D and 3D results (course and fine; all results in MPa)
Course (2D) Fine (2D) Course (3D) Fine (3D)
Max. Von M.
Min. txy or txza
Max. txy or txzb
Max. p
76.3 82.4 59.2 77.8
−35.9 −39.3 −21.5 −29.9
21.8 21.2 21.7 22.7
444 440 589 674
a
Shear stress at the surface opposite to the leg in the flow direction. b Shear stress at the leg in the flow direction.
These values are taken for the comparison in Table 3. Since, therefore, the shear stresses on other planes are neglected, it could be expected that the pressures are erroneously estimated by the 2D calculations. This can be concluded from Table 3 also. It was also observed that the calculations were highly sensitive to the definition of the contact in the outflow region. Since this influences the results quite significantly, at the moment contact algorithms are under constructions that may be less sensitive. For further results on the 2D calculations the reader is referred to [1,15]. Of course, from the 3D calculations results can be obtained that can not be achieved by means of the 2D calculations. An example is the nonuniform outflow velocity at the outlet of the die. This results in differences in outflow velocity of the profile that is formed, see Fig. 9. In this figure the difference of the displacement of points on the profile at the outlet are shown. The coordinate s is taken along the outer surface of the tube in the tangential direction, see Fig. 6. The different curves represent the displacements of points in the outlet at five different moments in time. It can be seen that points behind the leg in the outlet (arc length s = 32) have smaller displacements than points where the aluminium flows directly along the legs (arc length s = 0). If the profile were not symmetric, these length
difference could cause curved and/or twisted profiles. Therefore, these length differences are interesting for extrusion practice. Generally, it can be stated that the two-dimensional calculations can be used for the estimation of the loads that an extrusion die is exposed to. Furthermore, it is not necessary to use refined meshes because the results do not significantly improve upon this refinement. For three-dimensional calculations there is a difference in calculation times of a few hours or a few days. The contact conditions are an important factor for the results. Therefore, close attention should be paid, for example, to oscillatory behaviour from time step to time step because this indicates irregular behaviour of the contact elements. Since extrusion is essentially a thermal process, thermal degrees of freedom should also be incorporated. Because two-dimensional calculations are shown to give reasonable results, it suffices to
Fig. 7. Shear stress distribution for the 2D calculation.
Fig. 9. Length differences in outflow.
Fig. 8. Shear stress distribution in 3D calculation (MPa).
74
H.G. Mooi et al. / Journal of Materials Processing Technology 88 (1999) 67–76 Table 4 Material parameters for the die
Fig. 10. Mesh for the three-dimensional deformation calculation of the die (upper and side view).
expand only the two-dimensional calculations with thermal degrees of freedom if the order of magnitude of the difference between two- and three-dimensional calculations are known. An example of this is treated in [1] and [15].
5. Three-dimensional deformation calculations In order to complete the consideration on the extrusion of a tube, the three-dimensional deformation calculations on the extrusion die are briefly treated also. The results of the (thermal) flow calculations were used to determine the loads for these calculations. Different deformation mechanisms of the die were considered: i.e. elastic, creep and ratcheting deformation. For more details on the ratcheting calculations, the reader is referred to [1]. Because of symmetry, only one eighth of the porthole part of the die has to be modelled. The three-dimensional mesh with the boundary conditions is shown in Fig. 10. Some details are left out of the model, such as the exact bearing length, the exact shape of the legs and most of the roundings. These simplifications do not have a great influence on the results. The mesh contains 628 nodes and 353 eight node brick elements. A refined model was also used for the same type of calculations with 3404 nodes and 2535 elements. The model was simply supported at the outer two rows of nodes on the lower side. Since these boundary conditions (rigidly suppressed) were too stiff compared to the elastic support of the surrounding tooling, e.g. the press, a smaller part of the die was supported than in practice. At the symmetry planes, symmetry conditions were prescribed (upper view in Fig. 10). The pressure in front of the die was calculated by means of the two-dimensional calculations with temperature dependent material parameters like in Section 4. This pressure was prescribed as a distributed load on the core of the die (q1 in Fig. 10). The dis-
Elastic
E (MPa) n sy,0(MPa)
176 000 0.3 850
Creep
C (MPa−ms−n) m n
3.2×10−26 8.0 0.28
tributed load on the leg consisted of two parts. Firstly, the pressure difference over the leg contributed to this load. These pressures also follow from the two-dimensional flow calculations. Note that on the downstream part of the leg a pressure also works in the direction opposite to the extrusion direction. Secondly, the friction forces on the side of the leg were summed and modelled as a distributed load on the front of the leg. Both loads on the leg were combined in the distributed load q2. The distributed loads q1 and q2 were taken to be 230 and 350 MPa, respectively. The material models used, were purely elastic, elastoviscoplastic and viscous creep models. The latter two were treated in Section 2. The material parameters are taken according to Table 4. A first calculation is used to check whether the deformations are immediately plastic by means of the Von Mises stress distribution, see Figs. 11 and 12. For this calculation an elasto-plastic model only is used. The yield strength at 500°C for this material is 850 MPa, see Table 3. Therefore, no direct plastic strains are present for these loads. However, it is clear that if thermal loads are added or the loads are increased only slightly, plastic deformation might be possible, see e.g. [1,15]. The main conclusion of this first calculation is that the porthole part of the die does not directly deform plastically. In practice permanent deformation of the die is, in fact, observed. Therefore, it is also interesting to investigate the creep deformation of the die under the earlier mentioned loads. In order to do so, the same
Fig. 11. Von Mises stresses (in MPa; max 810 MPa).
H.G. Mooi et al. / Journal of Materials Processing Technology 88 (1999) 67–76
75
compared. Peaks in the extrusion load might result in even higher creep deformation. The viscous and elastoviscoplastic description of creep gave results in the same order of magnitude. The latter should be favoured because of its higher accuracy and the fact that loads can always be prescribed. The viscous model, though, can be used if a quick indication of the creep deformation is needed. In that case, the first step with (prescribed displacements) should be made sufficiently small. Both models tend to converge very badly, which might be cured by using multi-phase models as treated in [6]. Fig. 12. Deformed three-dimensional mesh (exaggerated).
problem was calculated with a material model for creep (see Section 2). The material parameters are given in Table 4. Both elasto-viscoplastic and viscous creep formulations were applied. The displacement in time of the point on the centre line on the upper surface of the die for the elasto-viscoplastic and the viscous models is given in Fig. 13. Note that the upper line also incorporates an elastic part. The viscous and the elasto-viscoplastic model give somewhat different results. This can mainly be explained by the fact that the creep strain for the viscous formulation is incorrect in the first step. This is due to the fact that displacements needed to be prescribed instead of loads for the sake of convergence. Note the steepness of the elasto-viscoplastic line due to the highly non-linear behaviour of the creep law; this requires very small time steps. The displacements by creep deformation are still about an order smaller than the elastic deformations that measured a few tenths of a millimetre, see also Fig. 13. It can thus be concluded that creep plays at least some role in the permanent deformation of the dies. In practice, permanent deformation of porthole dies is measured as a few millimetres, but the history of these dies over 180 h was unknown so that it could not be
6. Conclusions The ALE method, as developed by Hue´tink, is very applicable to the modelling of complex aluminium extrusion. It was shown that interesting insights could be gained in the aluminium flow and die deformation aspects of aluminium extrusion by separating both calculations. The two-dimensional flow calculations were shown to give very reasonable estimates of the real three-dimensional stress state and deformations. Therefore, it was tolerable to use the two-dimensional (thermal!) calculations to determine the pressures and shear forces that are exerted on the die during extrusion. The deformation of the die was shown to be mainly elastic; the permanent deformations could be ascribed to creep. Future research should focus on more complex geometries, for which three-dimensional possibilities of present FE codes should be enhanced anyway. Also, the friction conditions might need closer attention.
Acknowledgements The research project was carried out within the framework of the Eureka project no. EU 554, which was initiated by BOAL B.V. in De Lier, The Netherlands, and Phoenix s.p.a. in Paderno, Italy. The support is gratefully acknowledged.
References
Fig. 13. Displacement of the centre line of the die.
[1] H.G. Mooi, Finite Elements Simulations of aluminium extrusion, Ph.D. thesis, University of Twente, The Netherlands, 1996. [2] H.G. Mooi, J. Hue´tink, Thermodynamically based material models for large deformations with application to Aluminium Extrusion, J. Mech. Phys. Solids (submitted). [3] R. Akeret, H. Jung, G. Scharf, Atlas of hot working properties of nonferrous metals, vol. 1, Deutsche Gesellschaft fu¨r Metallkunde (DGM) 1978. [4] B. Verlinden, A. Suhadi, L. Delaey, A generalized constitutive equation for an AA6060 Aluminium Alloy, vol. 28, Scripta Metallurgica , 1993, pp. 1441 – 1446.
76
H.G. Mooi et al. / Journal of Materials Processing Technology 88 (1999) 67–76
[5] H. Berns, F. Pschenitzka, Das Kriechverhalten von Warmarbeitssta¨hlen mit 5% Chrom (in German), Zeitschrift Werkstofftechn., no. (11), Verlag Chemie, Germany, 1980, pp. 258 – 264. [6] J. Hue´tink, On the simulation of thermomechanical forming processes, Ph.D. thesis, University of Twente, The Netherlands, 1986. [7] J. Hue´tink, Analysis of metal forming proc. based on a comb. Euler-Lagrangian FE formulation, In: J.F.T. Pittman, R.D. Wood, J.M. Alexander, O.C. Zienkiewicz (Eds.), Num. meth. in indust. forming processes, Pineridge Press, Swansea, UK, 1982, pp. 501 – 509. [8] P. Hora, L. Tong, J. Reissner, The FE simulation of 3-D bulk forming processes using the ALE formulation, Proc. of Metal Forming Process Simulation in Industry, Kro¨plin et al. (Eds.), Baden-Baden, Germany, Werner Fastenrath Druckerei, Wuppertal, 1994, pp. 46 – 59. [9] J.-L. Chenot, M. Bellet, The ALE method for the numerical simulation of material forming processes, in: Simulation of Materials Processing: Theory, Methods and Applications, Proc. Numiform ‘95, Balkema, Rotterdam, 1994, pp. 39–48.
.
[10] F.P.T. Baaijens, An U-ALE formulation of 3-D unsteady viscoelastic flow, Int. J. Numer. Methods Eng. 36 (1993) 1115–1143. [11] B. Koren, A robust upwind discretization method for Advection, Diffusion and Source terms, in: Notes on Numerical Fluid Mechanics, 45 (1993) 117 – 138, Vieweg, Braunschweig. [12] H.C. Stoker, J. Hue´tink, The arbitrary Lagrangian-Eulerian method in large deformation problems, in: Proc. EccoMass ‘96, Wiley, Paris, 1996. [13] T.J.R. Hughes, L.P. Franca, G.M. Hulbert, A new finite element formulation for computational fluid dynamics: VIII. the Galerkin/ Least-Squares methods for advective-diffusion equations, Comp. Methods. Appl. Mech. Eng. 73 (1989) 173 – 189. [14] H.G. Mooi, J. Hue´tink, Simulation of complex aluminium extrusion using an arbitrary Eulerian Lagrangian formulation, in: Simulation of Materials Processing: Theory, Methods and Applications, Proc. Numiform ‘95, Balkema, Rotterdam, 1995, pp. 869 – 874. [15] H.G. Mooi, A.J. Den Bakker, K.E. Nilsen J. Hue´tink, Simulation of Aluminium Extrusion by means of the Finite Element Method (FEM), in: Proc. 6th Int. Al. Extr. Technol. Seminar, Chicago, 1996, pp. 67 – 73.