M ODELLING C LOSED -D IE F ORGING O PERATIONS U SING T OTAL L AGRANGIAN S MOOTH P ARTICLE H YDRODYNAMICS

Total Lagrangian Smooth Particle Hydrodynamics (TLSPH) has been applied to a set of non-trivial, commercially interesting forging examples. Being a mesh-free method, TLSPH can conveniently simulate processes having large deformation and material separation. Test cases were designed that were characterized by large material ﬂows having large changes in grain connectivity. The implementation used, Smooth Mach Dynamics (SMD), provided tunable simulation parameters that enabled the simulation to optimally match each case. The results showed that the TLSPH/SMD has the potential to model the metal forging process efﬁciently without numerical instabilities. Each case studied required adaptation of the simulation parameters to optimize the results.


Motivation
Forging is an important manufacturing operation because it is economical and produces parts having high strength, toughness, and ductility as a result of the control of grain formation and orientation. The motivation for forging simulation is to predict the outcome of each stage of the forging workflow in order to mitigate problems such as: • incomplete die filling • production and cracking of flash • undesirable grain structure • cold shuts • fracture during forging • excessive required forging pressures • high residual stress • poor billet design • undesirable material flow • excess operations needed to create a part Ability to predict these outcomes can help optimize the design of preform dies, forging operations and process changes such as alterations of billet temperature and subsequent heat treatment operations.

Background
Understanding material deformation and flow patterns during forging can provide key insights to the design of forging operations. Forging defects can lead to premature failure of structures. Predicting the generation of defects will allow forging operations to be optimized.
Forging processes are characterised by rapid deformations under large forces that result in large plastic deformations. When the conventional finite element methods (FEM) are used to model large strain problems, high deformation inevitably leads to mesh distortion and possible instability in the simulation because the topology and shape of the elements have profound effects on the quality and accuracy of FE solutions [1]. Severely distorted elements may even lead to a singular stiffness matrix (non-positive Jacobian) and consequent analysis failure [2]. These issues are compounded in forging simulations, where very high displacements may produce mesh inter-penetration and entanglement resulting in the loss of connectivity and topological integrity of the finite element network. In such cases, re-meshing part or all of the geometry is the only alternative solution.
A mesh-less method, such as spherical particle hydrodynamics (SPH), offers many advantages for modelling continuous high deformation processes such as forging. SPH is a particle based numerical method for solving partial differential equations. The problem geometry is idealized into 'particles' that represent specific material volumes. The method uses a local interpolation from the neighboring particles to construct continuous field approximations that can be used to discretise the problem space. Even though SPH has been traditionally used for mod-elling fluid flows, SPH has been adapted to a wide variety of solid mechanics problems [3]. Since SPH does not require an explicit grid or mesh structure to represent the problem geometry it eliminates time-consuming effort of mesh generation and refinement through an implicit mechanism. It also avoids the need to address the difficulties associated with traditional mesh-based methods in maintaining the integrity and quality of the mesh under large deformations thereby reducing associated errors. The grid-free nature of SPH makes this method ideally suited to modelling material forging processes that involve large deformations and changes in topological connectivity. Moreover, SPH intrinsically handles surface contacts during forging and naturally predicts selfcontacts.
To address the specific needs of solid mechanics, the total-Lagrangian SPH (TLSPH) was devised [4,5]. In order for TLSPH to be viable, several problems had to be solved [6] to ensure its convergence and stability in order that large deformations having plastic flow could be accurately simulated [6]. In lieu of explicity remeshing, TLSPH handles remeshing implicitly by resetting of the material reference frame in response to a displacement metric. The implementation used in this study was implemented and packaged within the Lammps simulation platform [7] as the package USER-SMD.

Outline
The paper will review the theory associated with the TLSPH method, discuss the simulation method and related parameters, introduce the four forging case studies, present the results of these simulations and finally discuss issues found during the study.

Theory
The conservation equations will be solved by means of a total-Lagrangian approach which uses a reference configuration to compute stresses and accelerations. Both conservation and constitutive equations are expressed in terms of the reference coordinates X, which are taken to be the coordinates of the initial, undeformed configuration. A mapping between the two systems describes the motion at time t: x = φ(X, t).
(1) The displacement u (relative to the reference configuration X) is given by: The associated deformation gradient is given by: which is the transformation matrix that describes the rotation and stretch relative to the reference configuration.
The equations for mass, impulse, and energy conservation are given by: here, ρ is the mass density, P the first Piola-Kirchhoff stress tensor, e is the internal energy, and ∇ is the gradient or divergence operator and symbol T denotes the transpose of a matrix. The subscript 0 refers to a quantity evaluated in the reference configuration, while the absence of this subscript means that the current configuration is to be used. J is the determinant of the deformation gradient F.

Total Lagrangian SPH Formulation
In the SPH numerical method, the interpolated value of a function A at any position x can be calculated from the discrete values A j on the surrounding particles by: Here, the sum is over all particles j within a radius 2h of x. W (x, h) is a spline based interpolation or smoothing kernel with radius 2h and V j is the volume fraction of the support particle j.
At the beginning of each time step, a shape matrix is computed as the sum of the dyadic product of two vectors, ∇W i (X ij ), the gradient of the kernel function and X j − X i , the displacement of support particle j relative to particle i in the reference (ie. material) configuration.
Here, W i (X ij ) is a kernel function which depends only on the scalar distance between particles i and j. V j denotes the supporting particle volume, and S i the support region of the particle i. The shape matrix reflects the distribution of the particles in the support region surrounding the particle i within the reference configuration and is used to normalize the deformation to a spherical shape so that first order consistency is maintained. This constraint is expressed: The deformation gradient and its time derivative are then obtained as: where x and v are the coordinates and velocities in the current configuration. The velocity gradient in the current configuration is then obtained as L i ←Ḟ i F −1 i . The Cauchy stress σ is ultimately obtained by time integration of a constitutive model which relates the Cauchy stress rate to the strain rate tensor D i ← 1 2 (L i + L T i ). To obtain the Cauchy stress, one first needs the polar decomposition of the deformation gradient F: which is used to remove the rigid rotation from the strain rate tensor: From there, one can then use the constitutive relations to calculate the Cauchy stress σ i .
The first Piola-Kirchhoff stress, P i = Jσ i F −T i is used to link the material frame to the current frame where the deformation related forces are calculated as: In addition to the deformation related force, there is a viscous force which becomes significant at higher velocities. This is given as: where , c ij the speed of sound between particles i and j, and α, β the shear and bulk viscosities respectively. Finally, the net force acting on particle i is:

Hourglass Control
Within the region neighboring any given particle, a degeneracy of zero energy motions can occur, the so-called hourglass effect. To eliminate this effect, an ideal particle separation based on the initial configuration is proposed: As the particles deviate from this ideal separation, a corrective force can be applied. This corrective force is designed to minimize the error metric: To enforce the hourglass condition and maintain proper angular momentum, the corrective force: is applied.

Constitutive Model
The stress is broken into its pressure p and deviatoric σ d components.
The shear modulus for an isotropic material is given as: In this study, a simple piecewise linear constitutive model is used where stress varies linearly with strain up to the yield point (in proportion with the elastic modulus) and varies linearly past the yield point (in proportion to the strain hardening modulus).

Plasticity Model
The radial return algorithm of Wilkins [8] is used to calculate the elastoplastic deformation. A trial deviatoric stress is first calculated assuming an initial elastic response. This gives an increment in plastic strain: Here, σ vm is the von Mises stress and σ y is the current yield stress. The increment in the plastic strain is: The yield stress increment is given as: where H is the hardening modulus. The updated deviatoric stress σ d is given by: where r s is the radial scale factor given as:

Algorithm
The algorithm for describing the stress driven motion of SMD particles given in [9] is given in Appendix A containing the steps to reset the material reference frame along with the radial return algorithm used to handle plastic deformation [9].

Forging Cases
The test cases have been chosen to address select problem areas mentioned in the introduction such as die filling and fracture. The simplest billet geometries (blocks, cylinders) were chosen to simplify the interpretation of material flow patterns.

Sliding Yoke
The sliding yoke was chosen due to its non-trivial topology and variation in depth and breadth along its long dimension. The regions between the legs of the yoke and around the slider pin are expected to be regions of high shear. Here, the particle topology will change and require a reset of the reference state.
The part is forged normal to its thickness (the smallest dimension), similar to an engine connecting rod and forged from a rectangular billet with a volume slightly larger than the part itself to allow for complete die filling. The billet is oriented lengthwise to the part and has a width slightly smaller than the part itself. The center of the finish die has a raised region that forms the space between the legs of the yoke fig. 2a. Preliminary simulations showed that this feature would cause an excessive separation (ie. fracture) of the billet away from the die cavity toward the flash region of the die. To address this, the hollowed section was smoothed to create a preform die fig. 2b.

Seamless Pipe (Reverse Extrusion)
The seamless pipe is awkward to model using FEM techniques due to the continuous changes in grain topology at the shear zone near the outer edge of the punch. The meshless approach eliminates these remeshing steps.
The die was modeled as two concentric cylinders, each open at one end fig. 3. The billet is cylindrical and the forging path along the axis of the two die cylinders. As the forging progresses, material is compressed and initially displaced outward followed by flow in the opposite direction of the punch motion. Preliminary simulations showed that a preform die was not neccesary.

Four arm offset joint
This part, fig. 4, is similar to a part modeled in previously [10] and will be reproduced for consistency. It is expected to exhibit frequent particle topology changes along with non-trivial material flows. As the material navigates the constricted and expanding regions, one expects perturbations in velocity and flowrate.
The combined two punch motion is collinear with the two aligned arms of the part fig. 4b. The operation requires the closure of the upper and lower die sections followed by the two cylindrical punches sliding toward the center. In this case one is expected to initially see compression along the long axis of the cylindrical billet. As the two punches move toward each other, the billet will buckle and fill the two orthogonal arms of the die.

Crankshaft Unit
This case has similarities to the yoke case, but with a wider variation of thickness in the direction of the punch motion. As a result, it can test the extent of die filling as a function of initial billet placement.
A cross plane crankshaft can be forged using a compound die with two forging planes oriented at right angles fig. 5b. Each plane has two die units each shaping a single crank throw. To handle the cross plane case, a full length billet stradling each crank across the main bearings is needed, otherwise one would rely on the fusion of pairs of oxidized surfaces. The construction of this billet will be addressed in later work.
Each crankshaft throw has a large depth variation and is expected to produce large material deformations fig. 6. The action of the forging press is along the shortest dimension of the part ie. normal to the plane of the crank Figure 6: Comparison of finish die and preform die.
throw. Due to the large hollowed out section in the center of the crank throw, a preform die, similar to the yoke case, is employed fig. 6d.

Simulation Strategy
Each part and and its corresponding die geometry was modeled and meshed using SALOME [11] and the die meshs exported as .stl files describing a 2D triangular surface mesh. The forging simulation was done using Lammps [7] with the SMD user package [9]. Each simulation instance had a corresponding Lammps script file where the material, simulation parameters and simulation control logic were defined. The die geometry was imported to Lammps and the billet geometry defined within Lammps itself. The material used was steel at different temperatures: three cases at room temperature and the 4arm case at 1000 • C. During each trajectory, the positions, velocities and relevant SMD properties (such as stress and strain) were written to a Lammps dump file. Other properties such as punch displacement and punch reaction force were computed and written out to data files to be graphed. Initial parameter settings and billet sizes were adjusted according to any issues that occured during the simulation. The volume of the part was used to estimate the initial billet size. The billet geometry was matched to the die geometry as closely as possible to prevent excess flash. In two cases, preform dies were modeled to address problematic material flows observed when finish dies were directly applied to the raw billet. This work will report the simulation of the preform dies since this is where most of the deformation occured.  Table 1 identifies the relevant material parameters. Here, steel was used where the modulus of elasticity and the yield stress varied with temperature according to [12]. Table 2 provides the material parameters by case.

Simulation Parameters
A set of tunable SMD parameters were used to control aspects of the simulation. These parameters are defined in tab. 3 and the settings given, by case, in tab. 4.

Presentation of results
The results are generated from the simulation trajectories using the Open Visualization Toolkit (OVITO) [13] and typically show the evolution of deformation and other parameters in space and time for each case. In the figures to follow, rainbow color coding is used where blue hues indicate a lower value while red hues indicate a higher value. For instance, strain values will vary within the range [−

Preform Die
The yoke system showed good die filling with a reasonable amount of flash. A lateral cross section about midway along the legs of the yoke fig. 1a was featured with respect to particle deformation and other properties. The evolution of a diagonal layer of particles was followed in fig. 7. The trajectory showed that the outer diagonal would align with the legs of the yoke and the center diagonal would flatten out fig. 7d. As the die closes, the orientation/elongation of the particles at the die surface becomes apparent.  The sectional strain is shown in fig. 8 at an intermediate stage of the forging process. The lateral strain xx is maximized near the midline of the billet where the central region is being carved out. The vertical strain yy is minimized near the midline of the billet. The shear yx zone appears at the corners where the central die section displaces the region between the yoke legs.
The plastic strain rate fig. 9 is also maximized in this central region where material is being displaced toward the sides. Initially it is concentrated at the surface and then spreads ubiquitously throughout the part as the die reaches the bottom of its stroke fig. 9. As the upper die descends, the strain rate increases on the lower billet as it is pressed against the lower die fig. 9a. From there the strain rate spreads throughout the billet and becomes confined to the midline 9b. As the base of the forging stroke is reached, the strain rate distributes throughout as the billet completely fills the die cavity fig. 9c.
The global distribution of properties is shown in fig. 10. The strain xx , normal to the major axis, given in fig. 10a, shows highest strain in the region between the legs of the yoke as the material is displaced to fill the legs of the part. The strain yy , in the direction of the punch motion, shows

Reverse Extrusion of Seamless Pipe
The evolution of a cylindrical material layer in the seamless pipe is shown in fig. 12. The layer is compressed and squeezed outward and then upward along the outer part of the die cylinder fig. 12d. The divergence of the grains from the center can be easily seen, fig. 12, as the material is displaced outward. The SMD implementation computes the shape and orientation of each particle fig. 14. This information allows the structure of the particle grains and their orientation to be rendered as an ellipse. In fig. 14 one sees the grains below the punch being compressed while those on the extruded pipe walls being elongated and aligned with the die wall.
Force vs displacement is shown in fig. 15. In the early stages there is a spike in force as the material is upset. As the punch descends the force increases and near the bottom there is a notable upward inflection where the last bit of material is squeezed out. The force becomes slightly negative as the punch rebounds.

Four Armed Extrusion
To visualize the results for the four arm part, two clipping planes (at right angles) are used to carve out a quadrant of the extrusion process fig. 16, showing the axial and lateral cross sectional evolution of the billet. The die is shown in a purple color to indicate the boundary of material flow fig. 16. The initial stages of extrusion exhibit a buckling phase as the billet is longitudinally compressed (not shown). This compression gives way to the lateral plastic flow of material into the lateral arms of the die fig. 16b. As the material flows through the constricted regions, the billet narrows fig. 16c. Eventually the material fills the die cavity fig. 16d.
The spatial-temporal distributions of selected billet properties are shown in fig. 17. The strain along the punch axis fig. 17a and the associated particle velocity fig. 17b indicate peak values in the constricted regions and lower  Step 125 Step 125 Figure 19: Final xx , yy of crank throw shown within preform die.
values where the radius of the part increases. In the lateral arms, v y and yy become more uniform in the direction of flow normal to the puch axis.
The nearest neighbor plot fig. 17c shows highest values at the center of the die as expected, but also shows a decrease in the center regions where the die has a larger radius, probably due to the expansion as material flows into these regions. The plasticity rate fig. 17d is maximized at the entrance to the die and at its outer surface, areas of flow expansion and high shear respectively. This is expected as the snapshot is taken late in the forging process as the material is filling out the die. The y-axis aligned grain orientation fig. 17e depicts how the particles are aligned with the punch axis. Grain alignment is maximized in the constricted regions and near the outer surface of the die. The plasticity fig. 17f correlates well with the plasticity rate. It is seen to be more coherent along the punch axis and conversely in the lateral directed die arm.
The force vs displacement is shown in fig. 18. In comparison to the seamless pipe extrusion fig. 15, it appears to be more noisy. The forging process is characterized by stages that result in non-uniform flows such as when material is expanding into thicker regions or where it is buckling just prior to changing direction. Once these bottlenecks were overcome, the punch force decreased. In the latter extrusion phases, the force increases to allow the material to fill the die cavity fully.

Preform Die
The forged crankunit is shown in the context of the die mesh surface fig. 19. The strain in the direction of die motion is shown in fig. 19a and the strain in the direction of the longest part dimension in fig. 19b. The strain, xx , is maximized at the counterweights opposite to the connecting rod bearing whereas the strain, yy , is maximized toward the rod bearing. This is expected because the material will tend to flow in the direction normal to the punch motion in the shallower part of the die cavity.
The evolution of a billet layer (normal to the punch motion) is shown in fig. 20. The billet is shown with a    fig. 21 where the surface is painted with the strain along the main axis. The increase in the strain around the central bearing and throughout the region between the counterweights is clearly visible fig. 21d. The evolution of zz throughout the interior of the crankunit, in its thickest part, is given in fig. 22. Here the effects around the central bearing and near the surface between the counterweights are further emphasized 22d.
The effect of the forging process on the thickest features of the part are highlighted in fig. 23 where xx is shown. One sees that the strain in this region is aligned with the axis of die motion. Finally the strain along the longest part dimension (direction of crank offset) is shown fig. 24. The strain is maximized near the parting line and toward the connecting rod bearing where the die is most shallow.

Discussion
Overall the simulations provided stable and credible results for the cases under study. The predicted deformations helped in identifying problem areas such as billet fracturing and incomplete die filling (data not shown) so that the test cases could be optimized. In some cases, however, the simulation parameters needed to be modified to adapt to the nuances of each case.

Effect of Simulation Parameters
Update frequency and contact stiffness had the greatest impact on the trajectories. The contact stiffness was increased to compensate for the stiffness of steel so the billet would not penetrate the die surface. The reference frame update and smoothing kernel radius influenced the out-put of the highest plasticity cases, namely the seamless pipe extrusion and the four arm extrusion (discussed below). Decreased particle size was also investigated (data not shown), but not surprisingly, was found to have a cubic impact the on simulation time.

Effect of Update Frequency
In order to handle large plastic deformations and maintain simulation stability [14], the initial configuration of the system was reset when the displacement of any one particle exceeded a tunable displacement criteria. The reset criteria occured when the displacement between any two particles exceeded 50 percent of the smoothing radius. This was notable in regions of high shear where the particle registry (contact topology) would change. When this threshold was exceeded, the initial configuration for all billet particles was reset to the current positions. The update would affect the cumulative strain but (by design) leave the stress and plastic deformation history intact. The part geometry and configuration of the forging system also had an impact on these updates. For instance, the yoke and crank-throw experienced few updates, while the seamless pipe and four arm part experienced several updates; the material flow in these latter cases being more complex. The reverse extrusion had shear bands where particles were disconnected from their neighbors by the punch action, similar to a cutting tool. Here, the effects were most apparent on the inner wall of the extruded pipe. The flow pattern within the four arm die exhibited diverted flows and flow through constrictions. Within these regions the particles would shear past each other so the configuration would frequently be reset.
When the particles were reset, the strain would also be reset. To get a sense for the overall strain, the Lammps [7] platform and the SMD implementation [6] were altered to compute a cumulative strain based on the initial state of the trajectory. Even though this strain was invalid around the areas where the reset criteria triggered, it nonetheless gave an indication of the overall deformation. To facilitate the display, the determinant of the deformation matrix F was output (in the dump data) and used to filter out particles [13] where the determinant was either negative or assumed excessive values.

Effect of Physical Parameters
The higher temperature case required lower average die forces as expected, compare figs. 15, 18. This can be understood from the data given in [12], which shows a 20 fold reduction in billet stiffness for the higher temperature. Also apparent was the greater material flow due to the lower yield stress.

Summary of Effects of Die Geometry
The simulations where the range of punch motion relative to the maximum billet extent was smallest seemed to be the best behaved. In these cases, the number of reference state resets was minimal. When the punch motion was along the longest axis of the billet, however, the reference state had to be reset many times. Also, it was noted in the 4arm case, that the billet would undergo large compressive deformation before actually flowing into the lateral die spaces. This was surprising considering the relative incompressibility of steel, even at high temperatures. This anomaly is probably due to a simplistic constitutive model and could be remedied by modeling incompressibility.

Conclusion
The Lammps implementation of Smooth Mach Dynamics was used to simulate four representative non-trivial metal forging operations. The cases covered a wide spectrum of material deformation and forging die geometry. The Total Lagrangian method provided credible predictions of the deformation process and proved useful in identifying problems occuring in the design of forging processes.

Appendix A
The algorithm for computing equations of motion of SMD particles due to stress forces from [9] is given here along with the steps needed for reference frame reset (ie. implicit remeshing). while t ≤ t total do run simulation for time interval 3: for all i do loop over all particles 4: v i ← v i + 1 2mi ∆tf i compute velocity 5: x i ← x i + ∆tv i compute new position 6: for all j ∈ S i do loop over all particles j in neighborhood of particle i 7: r ij ← |x i − x j | Compute distance from particle i to j 8: r 0 ij ← |X i − X j | Compute reference distance from particle i to j 9: if r ij − r 0 ij > T hresHold then If distance exceeds reference distance by tunable threshold 10: U pdateF lag ← 1 Set flag to update material frame 11: if U pdateF lag then 12: for all i do loop over all particles 13: X i ← x i Update reference positions 14: for all i do loop over all particles 15: F i ← 0 initialize deformation tensor 16:Ḟ i ← 0 initialize deformation rate tensor 17: K i ← 0 initialize shape normalization tensor 18: for all j ∈ S i do loop over all neighbors of particle i 19: update deformation tensor i from particle j 20:Ḟ i ←Ḟ i + V 0 j v ij ⊗ ∇W (X ij ) update deformation rate tensor i from particle j 21: update shape normalization tensor i from particle j for all j ∈ S i do loop over all neighbors of particle i 35: t ← t + ∆t increment time