- Email: [email protected]
The practicalities of ring rolling simulation for profiled rings
Journal of Materials Processing Technology 125–126 (2002) 619–625
The practicalities of ring rolling simulation for profiled rings K. Davey*, M.J. Ward1 Department of Mechanical, Aerospace and Manufacturing Engineering, UMIST, P.O. Box 88, Sackville Street, M60 1QD Manchester, UK Received 1 December 2001; accepted 6 January 2002
Abstract Ring rolling finite element simulation for profiled rings is currently not practicable with existing Lagrangian codes and on any realistic computational platform. The main reason for this is that a large number of incremental stages are typically required for completion of a full simulation and the computational cost per increment is high. The nature of ring rolling means that the amount of deformation that takes place in a given increment is relatively small when compared with typical metal forming processes. Small errors generated over each increment can grow to give highly inaccurate predictions after not so many ring revolutions. Inaccuracies manifest themselves in volume growth, incomplete profile filling, overall dimensional inaccuracy and inaccurate force and stress level predictions. This paper describes measures that make the analysis of profiled ring rolling a more practicable proposition. A model has been developed comprising the finite element flow formulation, an arbitrary Lagrangian–Eulerian update strategy, and a novel iterative solution scheme called the successive preconditioned conjugate gradient method. The focus of the paper is on methods developed to produce accurate predictions and that ensure numerical stability for an arbitrary high number of ring revolutions. The approach adopted includes methods for volume control, circular movement stability, circular interpolation techniques for accurate transfer of state variables, contact evolution, etc. The methods are validated by comparison with earlier experimental work and previously developed models for both pure radial, and radial– axial ring rolling. In addition, some results are presented for the analysis of commercially produced railway wheels and rings. # 2002 Elsevier Science B.V. All rights reserved. Keywords: Ring rolling; Finite elements; Arbitrary Lagrangian–Eulerian
1. Introduction The ring rolling process is widely used in the production of large diameter seamless rings, with applications in the aerospace and automotive industries, as well as for producing railway wheels and tyres. Its simplest form, pure radial ring rolling involves compression of a small diameter ringshaped workpiece between a mandrel and drive roll. This causes the ring to grow whilst allowing a profile to be developed on its surface. In the more widely used radial– axial process compression also takes place in the axial direction between conical–axial rolls. The radial–axial process is illustrated in Fig. 1. The nature of ring rolling makes many of the simplifying assumptions typically employed in the analysis of other rolling processes inapplicable. In particular, the process is non-steady state throughout and involves radial deformation between rolls of different diameter, only one of which is typically powered. *
Corresponding author. Tel.: þ161-236-3311; fax: þ161-228-7040. E-mail address: [email protected] (K. Davey). 1 Present address: Rolls Royce plc, Derby, UK.
Finite element (FE) simulation of this process is particularly difficult, suffering from excessive computational requirements and stability difficulties. The first attempts at solving this problem focused on 2D simplifications [1,2] or models in which only the roll gap region was considered [3,4]. Modelling the entire process using conventional Lagrangian FE codes is inefficient, involving highdensity meshes, while purely Eulerian methods can fail to accurately capture evolving geometry. As a result of this, most recent workers have favoured an arbitrary Lagrangian– Eulerian (ALE) based approach [5–11].
2. ALE flow formulation The ALE description is an attempt to combine the advantages of the Eulerian and Lagrangian formulation methods. In the ALE formulation, a reference system is defined such that co-ordinates are neither fixed in space nor attached to the material. This computational reference system (CRS) is used in conjunction with a second reference system, the material reference system (MRS) which behaves as a Lagrangian mesh. The ALE formulation permits inde-
0924-0136/02/$ – see front matter # 2002 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 4 - 0 1 3 6 ( 0 2 ) 0 0 3 5 6 - 4
620
K. Davey, M.J. Ward / Journal of Materials Processing Technology 125–126 (2002) 619–625
Fig. 1. Schematic illustration of the ring rolling process.
pendent mesh movement in the CRS, which in turn allows the use of a non-uniform mesh with element concentration in roll gap regions. In this paper, the operator split ALE formulation has been used in conjunction with the flow formulation. The flow formulation has been selected for this model as it is believed [9] to deliver good accuracy with a high level of computational efficiency. Under the ALE formulation, independent mesh movement in the CRS is permitted. The resulting coupled ALE method of Hu and Liu [6] involves the solution of this entire system, it is, however, convenient to decouple the system by adopting an operator split method first described by Benson [9]. This involves an initial Lagrangian step followed by a purely convective one, in one deformation increment. Details of this approach are described by Davey and Ward [10]. The process of uncoupling of the MRS and CRS requires that, in addition to the incremental updating procedure, convection must be taken into account in order to update state variables at CRS nodes. The convection of state variables in this way allows the CRS to be positioned conveniently throughout the simulation. In ring rolling a dense mesh is desired around the roll gap region.
3. Conditions of analysis The specific nature of the ring rolling process means that a number of issues, which do not arise in the simulation of other metal forming processes become important. These issues together with the approaches taken towards handling them in the model are described in this section. 3.1. The ALE update procedure Under the operator split ALE formulation, each increment proceeds exactly as for the pure Lagrangian description using the CRS until the end of the solution phase. The velocity field determined by this solution stage, v, relates to the movement of material, i.e. the MRS, as the Lagrangian description assumes that the mesh moves with the material and that the convective velocity, c is zero. Updating the CRS based on v would result in a Lagrangian approach where the dense mesh region would gradually move out of the roll
gap. It is desirable in ring rolling simulation that the CRS does not move in the circumferential direction. It is also required that the boundary of the CRS must coincide with the boundary of the material. Thus a convective velocity field, must be defined such that the CRS velocity field exhibits both of these properties, i.e. the resulting CRS should be stationary in the circumferential direction and have a boundary that is coincident with that of the material. The process of achieving this involves the definition of a cylindrical co-ordinate system for both MRS and CRS updates. This is described in detail by Davey and Ward [10] and Ward [11]. 3.2. Mandrel velocity In ring rolling the mandrel is idle and therefore its angular velocity is unknown. Yang and Kim [1] used a method to generate the rotation speed of the mandrel where the tangential mandrel speed is treated as an additional R unknown. Minimisation of the velocity functional F ¼ O Eð_eij Þ dV R R 2 ~ ~ v dS þ k ð_ e Þ dV, is required, where Eð_eij Þ is a t f r v ii Sf O power density function such that sij ¼ @E=@ e_ ij , kv is a large positive constant termed the penalty constant, e_ ii is the volumetric strain rate, Sf the surface in contact with the rolls, ~ vr the relative velocity vector at Sf and ~ tf the frictional stress. The functional F must be regarded as a function of both the vector of nodal velocities, V, and the magnitude of mandrel tangential velocity vmt. In this case the stationary condition is @F ¼ 0; @V
@F ¼0 @vmt
(1)
This equation is only applicable to a non-profiled cylindrical mandrel, since variations in mandrel radius would mean that a single tangential velocity was impossible. The method is, however, easily adapted to profiled rollers. The tangential velocity of any point on the mandrel is given by aI or2 tI, where aI is the perpendicular distance from point I to the mandrel axis, or2 the angular velocity of the mandrel and tI is a unit vector in the tangential direction of mandrel motion at point I. Given that in this analysis, rigid rollers are employed the mandrel radius is constant for a particular point on its surface. The unit tangent vector tI can also be determined from knowledge of the angular position of the point. Therefore, Eq. (1) can be rewritten in terms of angular velocity @F ¼ 0; @V
@F ¼0 @om
(2)
This represents a non-linear system of n þ 1 equations and n þ 1 unknowns, where n is equal to the number of degrees of freedom multiplied by the number of CRS nodes. Whilst the basis of this approach was described in 1988 [1], it is not the approach that has been adopted by all other workers in the field. As described by the authors [10], three clear approaches therefore present themselves, zero friction stationary
K. Davey, M.J. Ward / Journal of Materials Processing Technology 125–126 (2002) 619–625
interface, stationary interface with zero friction in the rotation direction and friction in the axial direction, and application of friction in all directions based on a computed mandrel speed. For clarity these approaches will be referred to as zero friction, axial friction and full friction, respectively, for the rest of this paper. In order to compare the three methods of handling the mandrel interface, examples of the pure radial ring rolling tests were repeated for a limited number of revolutions. The tests chosen involved rolling with a grooved mandrel fed at a rate of 0.175 mm/rev. The mandrel diameter in this test was relatively close to the inside diameter of the workpiece at startup, the ratio of workpiece inner diameter to mandrel diameter was 1.1 at the start of the test. The test involved simulation under the three friction conditions. Predicted geometry was then compared. The predicted profile of the axial edge of the workpiece at increasing thickness reductions are shown in Fig. 2 under each of the friction conditions. It can be seen that all three friction models give visually similar results, with the greatest difference in predictions between the zero and full friction approaches. It is easy to see why other authors have come to the conclusion that the zero friction model is an acceptable simplification. However, some accuracy loss does seem likely. Comparison of displacements indicates that in both simplified models, deviation from the predictions of the full friction approach is seen. 3.3. Velocity boundary conditions The rotation speed of a ring under rolling conditions depends on the speed of the drive roll. In order for ring and drive roll velocities to be matched the determination of a neutral curve on the drive roll is required. The existence of a neutral curve depends on torque generated by the drive roll. In practice, the neutral curve is determined at each increment and a match between roll and workpiece velocities is maintained by the application of boundary conditions at nodes remote from the roll gap(s) in the angular direction. The mechanism for this comparison and the limitations of this approach are described in detail by Ward [11].
Fig. 2. Axial profile for deformed radial section under different friction models after 9% thickness reduction.
621
3.4. Contact Rollers are defined in a piecewise manner using quadrilateral contact segments. Segments of this type can be easily generated from a 2D CAD representation, by first dividing the roller profile into linear segments and then sweeping this profile cylindrically around the roll axis. Three stages of contact modelling can be identified: Detection of contact. Imposition of contact boundary conditions. Control of contact after time integration. 3.4.1. Detection of contact The initial detection of contact between workpiece and die involves comparing the separation of every material node on the surface with each roller segment. If the separation is less than a pre-specified contact tolerance, then die is assumed to be in contact with the segment. The method used here for detecting contact has been specifically developed to take advantage of the simplified case where dies are known to be annular rollers. Essentially, three quantities are determined: the distance of the node away from the central axis of the roller, the angular position of the node relative to a reference plane, and the lateral position of the node along the central axis. These three quantities allow the existence of a state of contact between any node and roller segment to be detected. The method is described in detail by Ward [11]. 3.4.2. Imposition of contact boundary conditions Friction between material and roller surfaces is modelled in the way described by various authors such as Kobayashi et al. [12]. This requires a mixed boundary condition, in which displacement or velocity is defined normal to the roller surface via a Dirichlet condition, and stress is applied tangentially to the surface with a Neumann condition. Application of these conditions requires that the co-ordinate system for the contacting node is transformed into a local one based on tangential and normal directions. If, for a given contacting node I, the velocity vector in the global coordinate system is given by vI, then the corresponding vector in the transformed system, v0I is given by v0I ¼ TI vI in which TI is the transformation matrix for node I. The matrix TI is determined based on the roller segment that is associated with node I at the contact search stage. The matrix is assembled as TI ¼ ½ t1S t2S nS T where t1S ; t2S and nS are base-tangent vectors and a normal vector for the contacting die segment. These vectors can be easily generated based on the co-ordinates of the points that define the roller segment. The Dirichlet boundary condition for velocity at the contact interface is given by vnI ¼ nS vD in which vnI is the normal component of velocity for node I and vD the feed velocity of the roller. This ensures that the contacting node moves with the roller in its feed direction. The Neumann boundary condition is derived from discretisation of the surface integral of tangential frictional
622
K. Davey, M.J. Ward / Journal of Materials Processing Technology 125–126 (2002) 619–625
stress. For an infinitesimal area, dS over which frictional shear stress may be assumed uniform, the resultant force dfdS is given as mf s 2 jvr j arctan dfdS ¼ pffiffiffi dS (3) v0 3 p where s is the effective stress and mf is a constant friction factor. Force components from this friction model are in the opposite direction to those of the corresponding components of relative sliding velocity. The component of i i dfdS in direction i, dfdS is therefore given by dfdS ¼ pffiffiffi i ð2mf s= 3pÞ½arctanðjvr j=v0 Þ ðvr =jvr jÞ dS. Integration over the surface SC allows the traction force in direction i, fi as i Z 2mf s jvr j vr i pffiffiffi arctan f ¼ dS (4) v0 jvr j 3p SC A potential problem in determining this integral is that the true contacting area is unknown. The approach of Oh [13] was taken here, in which integration is performed over the surface of the contacting element where the relative sliding velocity field is approximated from the nodal values by the shape functions, i.e. vr ¼
ns X NI vrI
(5)
I¼1
where ns is the number of surface nodes for the element. The components of relative sliding velocity are given by the tangential relative sliding velocity in the global system. Initially, the relative velocities in the local system, v0r1 1 and v0r1 2 must be determined from the components of the material and die velocities, v0I 1 and v0I 2 , and v0D 1 and v0D 2 , respectively, as follows: 0 1 01 vrl vI v0D 1 ¼ (6) v0rl 2 v0I 1 v0D 2 the components can then be transformed into the global system as 0 11 0 1 vrI vrI @ v2 A ¼ ½ t 1 t 2
(7) rI I I v0rI 2 3 vrI This approach works well for forging problems where each die point is moving with the same velocity, and the contacting surface is made up of many elements. In this case the rotation of the rollers means that each point on the roller has its own velocity. This is relatively easy to incorporate once the rotational and feed components of velocity are identified. The rotational component can be readily determined from the angular position of node I relative to a vertical plane through the roller’s central axis, its separation from the central axis, and the angular velocity of the roller. If these terms are labelled v0Drot 1 and v0Drot 2 in each of the tangent
directions, then Eq. (6) can be adapted to include roll rotation as 0 1 01 vI ðv0Dlin 1 þ v0Dort 1 Þ vrl ¼ (8) v0I 2 ðv0Dlin 2 þ v0Dort 2 Þ v0rl 2 where v0Dlin 1 and v0Dlin 2 are the components of die velocity that result from linear feed of the dies in each of the tangential directions. The other key difference with rolling contact is that the contacting surface represents a smaller proportion of the surface area of the workpiece, particularly early in the process. The method described above will only determine vr for those elements which have all surface nodes in contact with the die. This represents a problem in rolling, as a large portion of the true contact area would be ignored due to elements only being in partial contact. As well as giving an under estimate of contact area and therefore frictional effects, this also brings about significant convergence problems in the nonlinear solution scheme when new elements move into complete contact. The situation has been improved here by a simple adaptation in which Eq. (5) becomes the sum over the number of contacting nodes, rather than all surface nodes. A force can then be determined for a partially contacting element by Eq. (4). Eq. (4) can be evaluated to good accuracy by integrating over the contact surface numerically using 16-point Gaussian quadrature for the element faces. 3.4.3. Control of contact after temporal integration Following the ALE update phase a contact search must be performed to establish the new state of contact between each surface node and each roller to determine the new state of contact for all surface nodes and to associate contacting nodes with the correct roller segment. If this search identifies all surface nodes as either being in acceptable contact with a roller segment, or free from contact, then the boundary conditions can be determined for the next increment. The updating process can however result in nodes moving into a position where they are penetrating a roller. In this case one of two alternative methods can be adopted to correct this situation. These are termed sub-incrementation and node projection. Sub-incrementation involves reversing the update stage, by determining a new time step which, when the explicit integration scheme is applied places the penetrating node on the roller surface. Thus, an Euler update with a negative time step based on the velocity field determined for the current increment is performed followed by a final correction subincrement, required so that the remaining time for the original increment size is accounted for. As with a full increment, the sub-increment stage requires a non-linear solution phase. In practice, several sub-increments can be required before a contact problem of this nature is resolved. This can have a significant influence on run time. The nodal projection method involves applying a corrective displacement to the penetrating node so that it is placed
K. Davey, M.J. Ward / Journal of Materials Processing Technology 125–126 (2002) 619–625
on the roller surface. Methods of this type have been employed by Chaudhary and Bathe [14], and Chenot and Bellet [15]. In both cases, different methods of determining the projection direction are employed. In this case the annular nature of the rollers is exploited, with the projection being applied in a radial direction from the roll axis. The projection method potentially suffers from volume loss problems, as the projection will, in all cases, move the penetrating node in such a way as to reduce workpiece thickness. However, corrections of this type are infrequent, particularly under the ALE formulation, as the contacting region changes only very gradually throughout the process. In addition, the corrections are generally very small. In view of this situation the run time benefits from this approach over the sub-incrementation approach have been exploited throughout the trial runs carried out here.
623
Fig. 3. Comparisons between current model and the experiments of Mamalis et al. [17] for axial spread: CL, centre line; ID, inside diameter; OD, outside diameter.
5.1. Tyre rolling trial 4. Solution method The successive preconditioned conjugate gradient method (SPCGM) is a technique developed by the authors [16] and later applied for ring rolling [8]. It essentially combines the preconditioned conjugate gradient method with the quasi-Newton (QN) method. The linear SPCGM algorithm is designed to construct the search vectors based on preconditioners for each iteration which are determined by a QN update step. The initial preconditioner, B0 is chosen in this case as the fully LDLT-factorised stiffness matrix. QN updating is only carried out under certain circumstances. These have been examined in detail in [16]. The objective is to minimise the number of updates, and therefore the cost, while ensuring that the condition number of the preconditioned problem remains at an appropriate level.
5. Model testing This model has been successfully applied to the experimental pure radial ring rolling and radial–axial ring rolling of tellurium–lead specimens by Mamalis et al. [17] and Salimi [18], respectively, in the previous work by the authors [10]. In particular, comparison between the results of this model and the work of Mamalis et al. [17] has demonstrated that the model is capable of predicting axial flow and ring growth for increasing radial compression. Typical comparisons taken form Davey and Ward [10] for the case of low feed rate pure axial rolling, with a grooved mandrel, can be seen for axial flow in Fig. 3. In the following sections, the model is applied to the industrial processes of railway tyre and wheel manufacture to demonstrate how increasingly complex geometry can be handled. The trials were performed on a 400 MHz Pentium II PC with 128 MB of RAM.
The rolling of railway tyres involves a non-symmetrical starting workpiece rolled between radial and axial rolls, and with both mandrel and drive rolls profiled. The lack of symmetry presents a problem in that, for rapid convergence, nodes need to be selected where velocity boundary conditions in the axial direction can be applied. For such problems, boundary conditions are imposed in an angular direction in order to maintain the position of the neutral curve in the basic model. It was therefore decided that this plane of nodes, removed as it was from the main deforming region, could also have axial direction boundary conditions imposed with minimal effect on the predicted geometry of the workpiece. In addition to this difference in the imposition of boundary conditions, the asymmetry of the tyre rolling process required both upper and lower axial rolls to be modelled. No modification to the basic model was required for this purpose. The mesh configuration of material and rollers used in the tyre simulation together with the plane of boundary condition nodes is illustrated in Fig. 4.
Fig. 4. Mesh configuration used in tyre simulation.
624
K. Davey, M.J. Ward / Journal of Materials Processing Technology 125–126 (2002) 619–625
Fig. 5. Effective strain rate distribution during tyre rolling. Fig. 7. Mesh configuration for wheel rolling simulation.
Two 308 dense mesh regions were employed for this simulation. These are located in the two roll gap regions so as to adequately represent the deformation in those most rapidly deforming areas. The effectiveness of the dense regions in covering the areas of greatest strain rate can be seen in Fig. 5, in which the effective strain rate distribution during the process is shown. The resulting mesh was made up of 1620 nodes. The process took place over 40 revolutions at a variable feed rate, this required 4800 deformation increments and took approximately 3 weeks runtime. 5.2. Wheel rolling trial The wheel rolling process differs from conventional ring rolling due to the lack of a central mandrel, the use of profiled axial rolls, and the presence of idle pressure rolls. The configuration is illustrated in Fig. 6. These differences meant that a modified version of the rolling simulation program was required. In particular, the need for mandrel speed calculation was replaced by the need to handle the interfaces between the pressure rolls and the workpiece. Given the results shown in Section 3.2, it was considered best to determine pressure roll velocities in the same way as
Fig. 6. Schematic illustration of wheel rolling.
was used for mandrel speed, and apply a full friction model at each interface. In addition, the presence of the pressure rolls meant that the dense mesh region on that side of the wheel needed to be extended to cover the entire region in close proximity to the pressure and edge rollers. This meant that a 908 dense region was used on that side of the wheel, and a 308 dense region used in the proximity of the drive roll as illustrated in Fig. 7. This clearly has resulted in an increased problem size, a situation that was partially reduced by excluding the boss section of the wheel from the model. Historical observations of the process has lead to the general assumption that the boss section is unaffected by the rolling process. This assumption suggests that the inner radius of the meshed region, the boundary with the boss can be assumed to be fixed during the process. Hence velocity boundary conditions have been applied in all directions to the inner ring of nodes in the mesh. The importance of extending the dense mesh region can be seen in Fig. 8, which shows the effective strain rate in the workpiece during the process. Strain rate in the vicinity of
Fig. 8. Effective strain rate distribution during wheel rolling.
K. Davey, M.J. Ward / Journal of Materials Processing Technology 125–126 (2002) 619–625
the drive roll is much smaller, suggesting that the dense region close to that roll could be made smaller.
6. Conclusions 1. The ALE update strategy allows the use of a nonuniform mesh and means that the dense mesh region does not follow the material out of the roll gap, as would be the case for a standard updated Lagrangian approach. 2. The mandrel rotation speed is determined as an additional unknown in the problem in an extension of the method of Yang and Kim [1]. Simplified methods of handling contact at this interface result in noticeably different predicted ring profiles. 3. The resulting model gives good predictions of both radial and axial flows, and ring growth during simulations of quite complex industrial problems. 4. While run time remains the limiting factor in the industrial application of this process, the combination of an operator split ALE formulation, with a number of simplifying assumptions leads to a highly convergent system of equations. When such a system is presented to the SPCGM solver, considerable run time benefits have been achieved.
References [1] D.Y. Yang, K.H. Kim, Rigid plastic finite element analysis of plane strain ring rolling, Int. J. Mech. Sci. 30 (1988) 541–550. [2] T.C. Tszeng, T. Altan, Investigation of ring rolling by pseudo-planestrain FEM analysis, J. Mater. Process. Technol. 27 (1991) 151–161; S.G. Xu, J.C. Lian, J.B. Hawkyard, Simulation of ring rolling using a rigid–plastic finite element model, Int. J. Mech. Sci. 33 (5) (1991) 393–401.
625
[3] S.G. Xu, J.C. Lian, J.B. Hawkyard, Simulation of ring rolling using a rigid–plastic finite element model, Int. J. Mech. Sci. 33 (5) (1991) 393–401. [4] D.Y. Yang, K.H. Kim, J.B. Hawkyard, Simulation of T-section profile ring rolling by the 3D rigid–plastic finite element method, Int. J. Mech. Sci. 33 (7) (1991) 541–550. [5] N. Kim, S. Machida, S. Kobayashi, Ring-rolling process simulation by three-dimensional finite element method, Int. J. Mach. Tools Manuf. 30 (4) (1990) 569–577. [6] Y.K. Hu, W.K. Liu, ALE finite element formulation for ring rolling analysis, Int. J. Numer. Meth. Eng. 33 (1992) 1217–1236. [7] T. Lim, I. Pillinger, P. Hartley, A finite-element simulation of profile ring rolling using a hybrid mesh model, J. Mater. Process. Technol. 80–81 (1998) 199–205. [8] K. Davey, M.J. Ward, An efficient solution method for finite element ring rolling simulation, Int. J. Numer. Meth. Eng. 47 (2000) 1997–2018. [9] D.J. Benson, An efficient, accurate, simple ALE method for nonlinear finite element programs, Comput. Meth. Appl. Mech. Eng. 72 (1989) 305–350. [10] K. Davey, M.J. Ward, A practical method for finite element ring rolling simulation using the ALE flow formulation, Int. J. Mech. Sci. 44 (2002) 165–190. [11] M.J. Ward, Practical models for ring rolling of railway wheels and tyres, Eng.D. Thesis, UMIST, 1999. [12] S. Kobayashi, S. Oh, T. Altan, Metal Forming and the Finite Element Method, Oxford University Press, Oxford, 1989. [13] S.I. Oh, Finite element analysis of metal forming processes with arbitrarily shaped dies, Int. J. Mech. Sci. 26 (1984) 165–176. [14] A.B. Chaudhary, K.J. Bathe, A solution method for static and dynamic analysis of three-dimensional contact problems with friction, Comput. Struct. 24 (1986) 855–873. [15] J.-L. Chenot, M. Bellet, The viscoplastic approach for the finite element modelling of metal-forming processes, in: P. Hartley, I. Pillinger, C. Sturgess (Eds.), Numerical Modelling of Material Deformation Processes—Research, Development and Applications, Springer, Berlin, 1992. [16] K. Davey, M.J. Ward, A successive preconditioned conjugate gradient method for the minimisation of quadratic and non-linear functions, IMACS J. Appl. Numer. Math. 35 (2000) 129–156. [17] A.G. Mamalis, J.B. Hawkyard, Johnson, Spread and flow patterns in ring-rolling, Int. J. Mech. Sci. 18 (1976) 11–16. [18] M. Salimi, Ph.D. Thesis, UMIST, 1988.
The practicalities of ring rolling simulation for profiled rings K. Davey*, M.J. Ward1 Department of Mechanical, Aerospace and Manufacturing Engineering, UMIST, P.O. Box 88, Sackville Street, M60 1QD Manchester, UK Received 1 December 2001; accepted 6 January 2002
Abstract Ring rolling finite element simulation for profiled rings is currently not practicable with existing Lagrangian codes and on any realistic computational platform. The main reason for this is that a large number of incremental stages are typically required for completion of a full simulation and the computational cost per increment is high. The nature of ring rolling means that the amount of deformation that takes place in a given increment is relatively small when compared with typical metal forming processes. Small errors generated over each increment can grow to give highly inaccurate predictions after not so many ring revolutions. Inaccuracies manifest themselves in volume growth, incomplete profile filling, overall dimensional inaccuracy and inaccurate force and stress level predictions. This paper describes measures that make the analysis of profiled ring rolling a more practicable proposition. A model has been developed comprising the finite element flow formulation, an arbitrary Lagrangian–Eulerian update strategy, and a novel iterative solution scheme called the successive preconditioned conjugate gradient method. The focus of the paper is on methods developed to produce accurate predictions and that ensure numerical stability for an arbitrary high number of ring revolutions. The approach adopted includes methods for volume control, circular movement stability, circular interpolation techniques for accurate transfer of state variables, contact evolution, etc. The methods are validated by comparison with earlier experimental work and previously developed models for both pure radial, and radial– axial ring rolling. In addition, some results are presented for the analysis of commercially produced railway wheels and rings. # 2002 Elsevier Science B.V. All rights reserved. Keywords: Ring rolling; Finite elements; Arbitrary Lagrangian–Eulerian
1. Introduction The ring rolling process is widely used in the production of large diameter seamless rings, with applications in the aerospace and automotive industries, as well as for producing railway wheels and tyres. Its simplest form, pure radial ring rolling involves compression of a small diameter ringshaped workpiece between a mandrel and drive roll. This causes the ring to grow whilst allowing a profile to be developed on its surface. In the more widely used radial– axial process compression also takes place in the axial direction between conical–axial rolls. The radial–axial process is illustrated in Fig. 1. The nature of ring rolling makes many of the simplifying assumptions typically employed in the analysis of other rolling processes inapplicable. In particular, the process is non-steady state throughout and involves radial deformation between rolls of different diameter, only one of which is typically powered. *
Corresponding author. Tel.: þ161-236-3311; fax: þ161-228-7040. E-mail address: [email protected] (K. Davey). 1 Present address: Rolls Royce plc, Derby, UK.
Finite element (FE) simulation of this process is particularly difficult, suffering from excessive computational requirements and stability difficulties. The first attempts at solving this problem focused on 2D simplifications [1,2] or models in which only the roll gap region was considered [3,4]. Modelling the entire process using conventional Lagrangian FE codes is inefficient, involving highdensity meshes, while purely Eulerian methods can fail to accurately capture evolving geometry. As a result of this, most recent workers have favoured an arbitrary Lagrangian– Eulerian (ALE) based approach [5–11].
2. ALE flow formulation The ALE description is an attempt to combine the advantages of the Eulerian and Lagrangian formulation methods. In the ALE formulation, a reference system is defined such that co-ordinates are neither fixed in space nor attached to the material. This computational reference system (CRS) is used in conjunction with a second reference system, the material reference system (MRS) which behaves as a Lagrangian mesh. The ALE formulation permits inde-
0924-0136/02/$ – see front matter # 2002 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 4 - 0 1 3 6 ( 0 2 ) 0 0 3 5 6 - 4
620
K. Davey, M.J. Ward / Journal of Materials Processing Technology 125–126 (2002) 619–625
Fig. 1. Schematic illustration of the ring rolling process.
pendent mesh movement in the CRS, which in turn allows the use of a non-uniform mesh with element concentration in roll gap regions. In this paper, the operator split ALE formulation has been used in conjunction with the flow formulation. The flow formulation has been selected for this model as it is believed [9] to deliver good accuracy with a high level of computational efficiency. Under the ALE formulation, independent mesh movement in the CRS is permitted. The resulting coupled ALE method of Hu and Liu [6] involves the solution of this entire system, it is, however, convenient to decouple the system by adopting an operator split method first described by Benson [9]. This involves an initial Lagrangian step followed by a purely convective one, in one deformation increment. Details of this approach are described by Davey and Ward [10]. The process of uncoupling of the MRS and CRS requires that, in addition to the incremental updating procedure, convection must be taken into account in order to update state variables at CRS nodes. The convection of state variables in this way allows the CRS to be positioned conveniently throughout the simulation. In ring rolling a dense mesh is desired around the roll gap region.
3. Conditions of analysis The specific nature of the ring rolling process means that a number of issues, which do not arise in the simulation of other metal forming processes become important. These issues together with the approaches taken towards handling them in the model are described in this section. 3.1. The ALE update procedure Under the operator split ALE formulation, each increment proceeds exactly as for the pure Lagrangian description using the CRS until the end of the solution phase. The velocity field determined by this solution stage, v, relates to the movement of material, i.e. the MRS, as the Lagrangian description assumes that the mesh moves with the material and that the convective velocity, c is zero. Updating the CRS based on v would result in a Lagrangian approach where the dense mesh region would gradually move out of the roll
gap. It is desirable in ring rolling simulation that the CRS does not move in the circumferential direction. It is also required that the boundary of the CRS must coincide with the boundary of the material. Thus a convective velocity field, must be defined such that the CRS velocity field exhibits both of these properties, i.e. the resulting CRS should be stationary in the circumferential direction and have a boundary that is coincident with that of the material. The process of achieving this involves the definition of a cylindrical co-ordinate system for both MRS and CRS updates. This is described in detail by Davey and Ward [10] and Ward [11]. 3.2. Mandrel velocity In ring rolling the mandrel is idle and therefore its angular velocity is unknown. Yang and Kim [1] used a method to generate the rotation speed of the mandrel where the tangential mandrel speed is treated as an additional R unknown. Minimisation of the velocity functional F ¼ O Eð_eij Þ dV R R 2 ~ ~ v dS þ k ð_ e Þ dV, is required, where Eð_eij Þ is a t f r v ii Sf O power density function such that sij ¼ @E=@ e_ ij , kv is a large positive constant termed the penalty constant, e_ ii is the volumetric strain rate, Sf the surface in contact with the rolls, ~ vr the relative velocity vector at Sf and ~ tf the frictional stress. The functional F must be regarded as a function of both the vector of nodal velocities, V, and the magnitude of mandrel tangential velocity vmt. In this case the stationary condition is @F ¼ 0; @V
@F ¼0 @vmt
(1)
This equation is only applicable to a non-profiled cylindrical mandrel, since variations in mandrel radius would mean that a single tangential velocity was impossible. The method is, however, easily adapted to profiled rollers. The tangential velocity of any point on the mandrel is given by aI or2 tI, where aI is the perpendicular distance from point I to the mandrel axis, or2 the angular velocity of the mandrel and tI is a unit vector in the tangential direction of mandrel motion at point I. Given that in this analysis, rigid rollers are employed the mandrel radius is constant for a particular point on its surface. The unit tangent vector tI can also be determined from knowledge of the angular position of the point. Therefore, Eq. (1) can be rewritten in terms of angular velocity @F ¼ 0; @V
@F ¼0 @om
(2)
This represents a non-linear system of n þ 1 equations and n þ 1 unknowns, where n is equal to the number of degrees of freedom multiplied by the number of CRS nodes. Whilst the basis of this approach was described in 1988 [1], it is not the approach that has been adopted by all other workers in the field. As described by the authors [10], three clear approaches therefore present themselves, zero friction stationary
K. Davey, M.J. Ward / Journal of Materials Processing Technology 125–126 (2002) 619–625
interface, stationary interface with zero friction in the rotation direction and friction in the axial direction, and application of friction in all directions based on a computed mandrel speed. For clarity these approaches will be referred to as zero friction, axial friction and full friction, respectively, for the rest of this paper. In order to compare the three methods of handling the mandrel interface, examples of the pure radial ring rolling tests were repeated for a limited number of revolutions. The tests chosen involved rolling with a grooved mandrel fed at a rate of 0.175 mm/rev. The mandrel diameter in this test was relatively close to the inside diameter of the workpiece at startup, the ratio of workpiece inner diameter to mandrel diameter was 1.1 at the start of the test. The test involved simulation under the three friction conditions. Predicted geometry was then compared. The predicted profile of the axial edge of the workpiece at increasing thickness reductions are shown in Fig. 2 under each of the friction conditions. It can be seen that all three friction models give visually similar results, with the greatest difference in predictions between the zero and full friction approaches. It is easy to see why other authors have come to the conclusion that the zero friction model is an acceptable simplification. However, some accuracy loss does seem likely. Comparison of displacements indicates that in both simplified models, deviation from the predictions of the full friction approach is seen. 3.3. Velocity boundary conditions The rotation speed of a ring under rolling conditions depends on the speed of the drive roll. In order for ring and drive roll velocities to be matched the determination of a neutral curve on the drive roll is required. The existence of a neutral curve depends on torque generated by the drive roll. In practice, the neutral curve is determined at each increment and a match between roll and workpiece velocities is maintained by the application of boundary conditions at nodes remote from the roll gap(s) in the angular direction. The mechanism for this comparison and the limitations of this approach are described in detail by Ward [11].
Fig. 2. Axial profile for deformed radial section under different friction models after 9% thickness reduction.
621
3.4. Contact Rollers are defined in a piecewise manner using quadrilateral contact segments. Segments of this type can be easily generated from a 2D CAD representation, by first dividing the roller profile into linear segments and then sweeping this profile cylindrically around the roll axis. Three stages of contact modelling can be identified: Detection of contact. Imposition of contact boundary conditions. Control of contact after time integration. 3.4.1. Detection of contact The initial detection of contact between workpiece and die involves comparing the separation of every material node on the surface with each roller segment. If the separation is less than a pre-specified contact tolerance, then die is assumed to be in contact with the segment. The method used here for detecting contact has been specifically developed to take advantage of the simplified case where dies are known to be annular rollers. Essentially, three quantities are determined: the distance of the node away from the central axis of the roller, the angular position of the node relative to a reference plane, and the lateral position of the node along the central axis. These three quantities allow the existence of a state of contact between any node and roller segment to be detected. The method is described in detail by Ward [11]. 3.4.2. Imposition of contact boundary conditions Friction between material and roller surfaces is modelled in the way described by various authors such as Kobayashi et al. [12]. This requires a mixed boundary condition, in which displacement or velocity is defined normal to the roller surface via a Dirichlet condition, and stress is applied tangentially to the surface with a Neumann condition. Application of these conditions requires that the co-ordinate system for the contacting node is transformed into a local one based on tangential and normal directions. If, for a given contacting node I, the velocity vector in the global coordinate system is given by vI, then the corresponding vector in the transformed system, v0I is given by v0I ¼ TI vI in which TI is the transformation matrix for node I. The matrix TI is determined based on the roller segment that is associated with node I at the contact search stage. The matrix is assembled as TI ¼ ½ t1S t2S nS T where t1S ; t2S and nS are base-tangent vectors and a normal vector for the contacting die segment. These vectors can be easily generated based on the co-ordinates of the points that define the roller segment. The Dirichlet boundary condition for velocity at the contact interface is given by vnI ¼ nS vD in which vnI is the normal component of velocity for node I and vD the feed velocity of the roller. This ensures that the contacting node moves with the roller in its feed direction. The Neumann boundary condition is derived from discretisation of the surface integral of tangential frictional
622
K. Davey, M.J. Ward / Journal of Materials Processing Technology 125–126 (2002) 619–625
stress. For an infinitesimal area, dS over which frictional shear stress may be assumed uniform, the resultant force dfdS is given as mf s 2 jvr j arctan dfdS ¼ pffiffiffi dS (3) v0 3 p where s is the effective stress and mf is a constant friction factor. Force components from this friction model are in the opposite direction to those of the corresponding components of relative sliding velocity. The component of i i dfdS in direction i, dfdS is therefore given by dfdS ¼ pffiffiffi i ð2mf s= 3pÞ½arctanðjvr j=v0 Þ ðvr =jvr jÞ dS. Integration over the surface SC allows the traction force in direction i, fi as i Z 2mf s jvr j vr i pffiffiffi arctan f ¼ dS (4) v0 jvr j 3p SC A potential problem in determining this integral is that the true contacting area is unknown. The approach of Oh [13] was taken here, in which integration is performed over the surface of the contacting element where the relative sliding velocity field is approximated from the nodal values by the shape functions, i.e. vr ¼
ns X NI vrI
(5)
I¼1
where ns is the number of surface nodes for the element. The components of relative sliding velocity are given by the tangential relative sliding velocity in the global system. Initially, the relative velocities in the local system, v0r1 1 and v0r1 2 must be determined from the components of the material and die velocities, v0I 1 and v0I 2 , and v0D 1 and v0D 2 , respectively, as follows: 0 1 01 vrl vI v0D 1 ¼ (6) v0rl 2 v0I 1 v0D 2 the components can then be transformed into the global system as 0 11 0 1 vrI vrI @ v2 A ¼ ½ t 1 t 2
(7) rI I I v0rI 2 3 vrI This approach works well for forging problems where each die point is moving with the same velocity, and the contacting surface is made up of many elements. In this case the rotation of the rollers means that each point on the roller has its own velocity. This is relatively easy to incorporate once the rotational and feed components of velocity are identified. The rotational component can be readily determined from the angular position of node I relative to a vertical plane through the roller’s central axis, its separation from the central axis, and the angular velocity of the roller. If these terms are labelled v0Drot 1 and v0Drot 2 in each of the tangent
directions, then Eq. (6) can be adapted to include roll rotation as 0 1 01 vI ðv0Dlin 1 þ v0Dort 1 Þ vrl ¼ (8) v0I 2 ðv0Dlin 2 þ v0Dort 2 Þ v0rl 2 where v0Dlin 1 and v0Dlin 2 are the components of die velocity that result from linear feed of the dies in each of the tangential directions. The other key difference with rolling contact is that the contacting surface represents a smaller proportion of the surface area of the workpiece, particularly early in the process. The method described above will only determine vr for those elements which have all surface nodes in contact with the die. This represents a problem in rolling, as a large portion of the true contact area would be ignored due to elements only being in partial contact. As well as giving an under estimate of contact area and therefore frictional effects, this also brings about significant convergence problems in the nonlinear solution scheme when new elements move into complete contact. The situation has been improved here by a simple adaptation in which Eq. (5) becomes the sum over the number of contacting nodes, rather than all surface nodes. A force can then be determined for a partially contacting element by Eq. (4). Eq. (4) can be evaluated to good accuracy by integrating over the contact surface numerically using 16-point Gaussian quadrature for the element faces. 3.4.3. Control of contact after temporal integration Following the ALE update phase a contact search must be performed to establish the new state of contact between each surface node and each roller to determine the new state of contact for all surface nodes and to associate contacting nodes with the correct roller segment. If this search identifies all surface nodes as either being in acceptable contact with a roller segment, or free from contact, then the boundary conditions can be determined for the next increment. The updating process can however result in nodes moving into a position where they are penetrating a roller. In this case one of two alternative methods can be adopted to correct this situation. These are termed sub-incrementation and node projection. Sub-incrementation involves reversing the update stage, by determining a new time step which, when the explicit integration scheme is applied places the penetrating node on the roller surface. Thus, an Euler update with a negative time step based on the velocity field determined for the current increment is performed followed by a final correction subincrement, required so that the remaining time for the original increment size is accounted for. As with a full increment, the sub-increment stage requires a non-linear solution phase. In practice, several sub-increments can be required before a contact problem of this nature is resolved. This can have a significant influence on run time. The nodal projection method involves applying a corrective displacement to the penetrating node so that it is placed
K. Davey, M.J. Ward / Journal of Materials Processing Technology 125–126 (2002) 619–625
on the roller surface. Methods of this type have been employed by Chaudhary and Bathe [14], and Chenot and Bellet [15]. In both cases, different methods of determining the projection direction are employed. In this case the annular nature of the rollers is exploited, with the projection being applied in a radial direction from the roll axis. The projection method potentially suffers from volume loss problems, as the projection will, in all cases, move the penetrating node in such a way as to reduce workpiece thickness. However, corrections of this type are infrequent, particularly under the ALE formulation, as the contacting region changes only very gradually throughout the process. In addition, the corrections are generally very small. In view of this situation the run time benefits from this approach over the sub-incrementation approach have been exploited throughout the trial runs carried out here.
623
Fig. 3. Comparisons between current model and the experiments of Mamalis et al. [17] for axial spread: CL, centre line; ID, inside diameter; OD, outside diameter.
5.1. Tyre rolling trial 4. Solution method The successive preconditioned conjugate gradient method (SPCGM) is a technique developed by the authors [16] and later applied for ring rolling [8]. It essentially combines the preconditioned conjugate gradient method with the quasi-Newton (QN) method. The linear SPCGM algorithm is designed to construct the search vectors based on preconditioners for each iteration which are determined by a QN update step. The initial preconditioner, B0 is chosen in this case as the fully LDLT-factorised stiffness matrix. QN updating is only carried out under certain circumstances. These have been examined in detail in [16]. The objective is to minimise the number of updates, and therefore the cost, while ensuring that the condition number of the preconditioned problem remains at an appropriate level.
5. Model testing This model has been successfully applied to the experimental pure radial ring rolling and radial–axial ring rolling of tellurium–lead specimens by Mamalis et al. [17] and Salimi [18], respectively, in the previous work by the authors [10]. In particular, comparison between the results of this model and the work of Mamalis et al. [17] has demonstrated that the model is capable of predicting axial flow and ring growth for increasing radial compression. Typical comparisons taken form Davey and Ward [10] for the case of low feed rate pure axial rolling, with a grooved mandrel, can be seen for axial flow in Fig. 3. In the following sections, the model is applied to the industrial processes of railway tyre and wheel manufacture to demonstrate how increasingly complex geometry can be handled. The trials were performed on a 400 MHz Pentium II PC with 128 MB of RAM.
The rolling of railway tyres involves a non-symmetrical starting workpiece rolled between radial and axial rolls, and with both mandrel and drive rolls profiled. The lack of symmetry presents a problem in that, for rapid convergence, nodes need to be selected where velocity boundary conditions in the axial direction can be applied. For such problems, boundary conditions are imposed in an angular direction in order to maintain the position of the neutral curve in the basic model. It was therefore decided that this plane of nodes, removed as it was from the main deforming region, could also have axial direction boundary conditions imposed with minimal effect on the predicted geometry of the workpiece. In addition to this difference in the imposition of boundary conditions, the asymmetry of the tyre rolling process required both upper and lower axial rolls to be modelled. No modification to the basic model was required for this purpose. The mesh configuration of material and rollers used in the tyre simulation together with the plane of boundary condition nodes is illustrated in Fig. 4.
Fig. 4. Mesh configuration used in tyre simulation.
624
K. Davey, M.J. Ward / Journal of Materials Processing Technology 125–126 (2002) 619–625
Fig. 5. Effective strain rate distribution during tyre rolling. Fig. 7. Mesh configuration for wheel rolling simulation.
Two 308 dense mesh regions were employed for this simulation. These are located in the two roll gap regions so as to adequately represent the deformation in those most rapidly deforming areas. The effectiveness of the dense regions in covering the areas of greatest strain rate can be seen in Fig. 5, in which the effective strain rate distribution during the process is shown. The resulting mesh was made up of 1620 nodes. The process took place over 40 revolutions at a variable feed rate, this required 4800 deformation increments and took approximately 3 weeks runtime. 5.2. Wheel rolling trial The wheel rolling process differs from conventional ring rolling due to the lack of a central mandrel, the use of profiled axial rolls, and the presence of idle pressure rolls. The configuration is illustrated in Fig. 6. These differences meant that a modified version of the rolling simulation program was required. In particular, the need for mandrel speed calculation was replaced by the need to handle the interfaces between the pressure rolls and the workpiece. Given the results shown in Section 3.2, it was considered best to determine pressure roll velocities in the same way as
Fig. 6. Schematic illustration of wheel rolling.
was used for mandrel speed, and apply a full friction model at each interface. In addition, the presence of the pressure rolls meant that the dense mesh region on that side of the wheel needed to be extended to cover the entire region in close proximity to the pressure and edge rollers. This meant that a 908 dense region was used on that side of the wheel, and a 308 dense region used in the proximity of the drive roll as illustrated in Fig. 7. This clearly has resulted in an increased problem size, a situation that was partially reduced by excluding the boss section of the wheel from the model. Historical observations of the process has lead to the general assumption that the boss section is unaffected by the rolling process. This assumption suggests that the inner radius of the meshed region, the boundary with the boss can be assumed to be fixed during the process. Hence velocity boundary conditions have been applied in all directions to the inner ring of nodes in the mesh. The importance of extending the dense mesh region can be seen in Fig. 8, which shows the effective strain rate in the workpiece during the process. Strain rate in the vicinity of
Fig. 8. Effective strain rate distribution during wheel rolling.
K. Davey, M.J. Ward / Journal of Materials Processing Technology 125–126 (2002) 619–625
the drive roll is much smaller, suggesting that the dense region close to that roll could be made smaller.
6. Conclusions 1. The ALE update strategy allows the use of a nonuniform mesh and means that the dense mesh region does not follow the material out of the roll gap, as would be the case for a standard updated Lagrangian approach. 2. The mandrel rotation speed is determined as an additional unknown in the problem in an extension of the method of Yang and Kim [1]. Simplified methods of handling contact at this interface result in noticeably different predicted ring profiles. 3. The resulting model gives good predictions of both radial and axial flows, and ring growth during simulations of quite complex industrial problems. 4. While run time remains the limiting factor in the industrial application of this process, the combination of an operator split ALE formulation, with a number of simplifying assumptions leads to a highly convergent system of equations. When such a system is presented to the SPCGM solver, considerable run time benefits have been achieved.
References [1] D.Y. Yang, K.H. Kim, Rigid plastic finite element analysis of plane strain ring rolling, Int. J. Mech. Sci. 30 (1988) 541–550. [2] T.C. Tszeng, T. Altan, Investigation of ring rolling by pseudo-planestrain FEM analysis, J. Mater. Process. Technol. 27 (1991) 151–161; S.G. Xu, J.C. Lian, J.B. Hawkyard, Simulation of ring rolling using a rigid–plastic finite element model, Int. J. Mech. Sci. 33 (5) (1991) 393–401.
625
[3] S.G. Xu, J.C. Lian, J.B. Hawkyard, Simulation of ring rolling using a rigid–plastic finite element model, Int. J. Mech. Sci. 33 (5) (1991) 393–401. [4] D.Y. Yang, K.H. Kim, J.B. Hawkyard, Simulation of T-section profile ring rolling by the 3D rigid–plastic finite element method, Int. J. Mech. Sci. 33 (7) (1991) 541–550. [5] N. Kim, S. Machida, S. Kobayashi, Ring-rolling process simulation by three-dimensional finite element method, Int. J. Mach. Tools Manuf. 30 (4) (1990) 569–577. [6] Y.K. Hu, W.K. Liu, ALE finite element formulation for ring rolling analysis, Int. J. Numer. Meth. Eng. 33 (1992) 1217–1236. [7] T. Lim, I. Pillinger, P. Hartley, A finite-element simulation of profile ring rolling using a hybrid mesh model, J. Mater. Process. Technol. 80–81 (1998) 199–205. [8] K. Davey, M.J. Ward, An efficient solution method for finite element ring rolling simulation, Int. J. Numer. Meth. Eng. 47 (2000) 1997–2018. [9] D.J. Benson, An efficient, accurate, simple ALE method for nonlinear finite element programs, Comput. Meth. Appl. Mech. Eng. 72 (1989) 305–350. [10] K. Davey, M.J. Ward, A practical method for finite element ring rolling simulation using the ALE flow formulation, Int. J. Mech. Sci. 44 (2002) 165–190. [11] M.J. Ward, Practical models for ring rolling of railway wheels and tyres, Eng.D. Thesis, UMIST, 1999. [12] S. Kobayashi, S. Oh, T. Altan, Metal Forming and the Finite Element Method, Oxford University Press, Oxford, 1989. [13] S.I. Oh, Finite element analysis of metal forming processes with arbitrarily shaped dies, Int. J. Mech. Sci. 26 (1984) 165–176. [14] A.B. Chaudhary, K.J. Bathe, A solution method for static and dynamic analysis of three-dimensional contact problems with friction, Comput. Struct. 24 (1986) 855–873. [15] J.-L. Chenot, M. Bellet, The viscoplastic approach for the finite element modelling of metal-forming processes, in: P. Hartley, I. Pillinger, C. Sturgess (Eds.), Numerical Modelling of Material Deformation Processes—Research, Development and Applications, Springer, Berlin, 1992. [16] K. Davey, M.J. Ward, A successive preconditioned conjugate gradient method for the minimisation of quadratic and non-linear functions, IMACS J. Appl. Numer. Math. 35 (2000) 129–156. [17] A.G. Mamalis, J.B. Hawkyard, Johnson, Spread and flow patterns in ring-rolling, Int. J. Mech. Sci. 18 (1976) 11–16. [18] M. Salimi, Ph.D. Thesis, UMIST, 1988.