© AFM, EDP Sciences 2020

1 Introduction

Metal forging, considered to be one of the oldest known metal forming techniques to mankind ( ~5000BCE), is a process in which a block of metal having a simple shape (billet) is shaped under impact loading or pressure in between tools (die and punch) of desired shape, with or without the influence of temperature. Forging process usually produces little or no scrap and is able to produce the final geometry in a very short time. In addition, for a given weight, forged parts exhibit better mechanical and metallurgical properties and reliability than parts produced by casting or machining due to grain refinement and a favourable grain orientation. Thus, forged parts find large use in applications that require high performance and reliability such as the aerospace industry and in large production industries such as the automotive industry where a large number of identical parts with a high strength to weight ratio is required. With increasing complexity of the final shapes that need to be achieved, most industrial forging operations are quite complex and require two or more stages [1]. Hence, to increase the quality and the efficiency of production, it is necessary to have an optimal initial shape (preform).

Finite element simulation of metal forming processes started in the late 1960s, mainly in academic laboratories for 2D work-pieces [2]. A series of developments in the large deformation elasto-plastic analysis in the following decades lead to the development of simulation codes for industrial forming applications. A review of these developments and their application to metal forming processes have been presented in [3–5] and an extensive bibliography of finite element modelling of bulk material forming maybe found in [6]. All these developments were concerned with the direct analysis (forward analysis) of the metal forming problem using an incremental method, starting from a given initial shape, to simulate the large plastic deformation and the contact friction in a step-by-step method. Due to the complex material constitutive laws, the influence of contact and friction between the tool and work pieces, large deformations and the influence of temperature, metal forming problems are highly non-linear in nature and as a consequence the direct analysis requires large computational effort.

As the major objective in metal-forming processing is the optimum selection of process conditions in the design stage of processes given the material state and the geometry of the final product, conditions on the initial work-piece and possibly some restrictions on the processes, the required process conditions are required inputs for the direct analysis. The design of forming processes usually consists of the repeated trial and error use of direct modeling techniques where the initial shape forms are defined by qualitative guidelines; derived mainly from experience or experimental studies. The preform design and the die design problems can also be formulated under a rigorous mathematical basis by posing them as optimization problems and can be solved with forward analysis methods to find optimal solutions. Approaches based on sensitivity analysis [7], parameter-based B-spline optimization [8,9], parameter free shape optimization [10], meta-model assisted methods [11] and approaches based on neural networks [12,13] among others have been employed in this regard. Although simplistic in nature, the widespread use of the optimization techniques has been limited due to the computational expenses of the forward analysis for each design iteration. Another method proposed to solve this inverse shape finding problem has been the “Backward Tracing Method” [14]. Although quite useful in elastic deformation, the method is ill-posed in a plastic deformation problem and the inverse solution is not unique unless the loading history and plastic variables are known. Recursive formulations can be employed in this regard [15,16], but makes the method computationally expensive.

A simplified method for a single-step analysis between the known final shape and unknown initial shape named “Inverse Approach” (IA) was proposed in the context of modelling of sheet metal forming [17] and later applied to axisymmetrical cold-forging simulation [18]. Although the method could provide results quickly, the accuracy of the predicted stress and strains were directly related to the predicted initial shape used and relatively poor compared to the forward analyses. A multi-step version of the same method named Pseudo Inverse Approach (PIA) was introduced later in both sheet metal forming [19] and axisymmetric cold forging [20]. The PIA has been able to keep some of the advantages of IA with respect to computational effort while having a better accuracy. When deformation takes place at high temperatures, material properties can vary considerably with temperature. Heat is generated during a metal-forming process and if the dies are at a lower temperature, heat loss by conduction to the dies and by convection and radiation to the environment can cause severe temperature gradients within the work-piece. Since materials at elevated temperatures are usually rate sensitive, a complete analysis of hot forming operations requires a consideration of the effect of rate-sensitivity (viscoplasticity) of materials and the coupling of the metal-flow and thermal analyses [21,22]. The present work extends the PIA for the numerical modelling of hot-forging of 2D axisymmetric components by taking into consideration both viscoplastic and thermal behaviour so that it can be a suitable method to be used in optimization routines used for forging process design.

2 Pseudo Inverse Approach

2.1 PIA – outline

Classical direct methods used for numerical simulation of the forging process are based on a Lagrangian/Updated Lagrangian formulation where the equilibrium equations are based on an undeformed/known equilibrium configuration. Since linearity needs to be maintained in the evolving contact boundary conditions and the force boundary conditions between the deformed and undeformed configurations, it implicitly brings a restriction on the maximum step sizes that can be used even though no assumptionsare made in the strain measures used. This is the primary reason for relatively small step sizes being used even in finite deformation problems based on the Lagrangian/Updated Lagrangian approaches. A way to overcome this restriction would be to solve the equilibrium equations based on the current deformed configuration following a Eulerian-type approach. The primary difficulty in this kind of an approach in a general problem is the prediction of a suitable deformed configuration which is close to the actual deformed configuration. The forging process, being an inverse-shape-finding problem, is an interesting conundrum in this context since, unlike general problems, the target final deformed shape is known at the beginning of the analysis. By using information from both the required target shape and the initial billet shape, PIA forms a hybrid scheme that easily predicts a close-enough deformed configuration that can be used in a large-step Eulerian-type approached. The salient steps in the PIA can be enumerated as below.

• 1.

  Create a suitable mesh on the target final deformed shape. Since the target shape of the part is generally known at the onset, any suitable meshing program can be used to create a mesh on the target shape.

• 2.

  Generate a suitable initial shape using a single step IA or from the thumb-rules used by forging process designers.

• 3.

  Transfer the mesh created on the target deformed shape to the generated initial shape using mesh-morphing techniques. In the case of an axisymmetric component, the position of the boundary nodes of the target mesh can be interpolated on the boundary of the initial shape created in step 2. Using the displacements of these boundary nodes as Dirichlet boundary conditions, we can solve a linear elastic problem to find the positions of the interior nodes of the target mesh in the given initial shape. Similar methods maybe used in the case of 3D components as well.

• 4.

  For a given time-step, geometrically interpolate the possible deformed configuration and the tool positions at that time step.

• 5.

  Use free surface correction to ensure that the incompressibility constraint and boundary conditions are satisfied, and generate a kinematically admissible configuration from the geometrically interpolated configuration.

• 6.

  Solve for the non-linear mechanical and thermal equilibrium using the kinematically admissible configuration to obtain the deformed configuration at that time step.

• 7.

  Repeat steps 4, 5, 6 till the entire time-range of deformation has been accounted for.

2.2 Generation of intermediate configurations

The creation of intermediate shapes in the Pseudo Inverse Approach is based on a geometrical interpolation between the initial/intermediate shapes and the target shape. This imposes a restriction that the finite element meshes used in both the target shape and the initial shape should be consistent (equal number of nodes and the same element connectivity). This restriction is related only to the geometric proportional interpolation that is being used presently. In the interest of avoiding computational complexity mesh refinement has not been considered, but it maybe easily incorporated into the method very easily, albeit at the cost of increasing the computational complexity and time as is usually the case with classical mesh refinement methods.

The initial shape is either assumed or found out from a single-step inverse approach. The shape obtained from the single step inverse approach has been shown to have a good approximation for shapes and strains even though the stress predictions are not quite accurate. The finite element mesh of the desired part is then mapped onto the initial shape. This is done by mapping the contour nodes of the mesh M^fin of the desired shape onto the contour C⁰ of the initialshape. This can be efficiently done using mesh morphing algorithms in use. The interior nodes of the mesh M⁰ can then be determined by a linear solution of the mechanical equilibrium with the imposed displacements on the boundary nodes as the boundary conditions.

Once the conforming meshes have been created in the initial shape, the intermediate configuration can be created by using a geometric proportional interpolation of the form, $$x^{t+\Delta t}=x^t+\eta \left(x^{\text{fin}}-x^t\right)$$(1)

where x^t+Δt is the nodal position of the configuration at time t + Δt, x^t is the nodal position of the configuration at time t, x^fin is the nodal positions of the final configuration and η is a constant between [0,1] and is dependent on the total no. of steps that is considered in the simulation. More details concerning the generation of the intermediate configurations maybe found in [23].

2.3 Free surface correction and boundary conditions

One limitation of this geometrically interpolated configuration is that the boundary conditions (contact conditions) with respect to the tool positions may not be respected on the contour. To correct this limitation and to obtain a kinematically admissible solution that is close to the statically admissible solution, the boundary nodes of the workpiece that penetrate the tool surface are projected back onto the tool surface using a node-to-element projection and subsequently a linear solution is carried out to determine the new positions of the internal nodes. Assuming rigid tools, it is trivial to compute the tool position corresponding to the time-step under consideration. Since these nodes are projected onto the tool surface, these are implicitly assumed to be in contact with the tool. But in a practical situation this need not be true and the contact conditions at the surface nodes needs to be verified. For this, we make use of the free-surface method. In a metal-forming operation simulation, for any configuration Ω, the free-surface condition can be expressed as, $$\sigma n=t\le 0\text{\hspace{1em}on}\partial \Omega$$(2)

where t is the traction force in the outward normal direction from the workpiece surface. This means that there is a compressive loading on the boundary in contact with the tool and no loading on the surface that is free from contact with the tool surface. To verify that the free surface conditions are respected, the nodal forces are computed on the geometrically interpolated/corrected configuration. Using u₁, u₂ and u₃ to denote the components of the nodal displacement in the normal and tangential directions respectively, the following cases can be considered.

• 1.

  If t ≥ 0, then the node belongs to the free surface and the boundary conditions on the node can be represented as u₁ ≠ 0; u₂≠0; u₃≠0.

• 2.

  If t < 0, then the node is in contact with the tool surface (compressive force) and the boundary conditions on the node can be represented as u₁ = 0; u₂≠0; u₃≠0.

The linear solution is done iteratively till the contact conditions are also satisfied and the kinematically admissible configuration is obtained.

Once a kinematically admissible configuration is obtained, the statically admissible equilibrium configuration that satisfies both the thermal and mechanical equilibrium is found out using the mathematical formulation of the PIA. For maintaining the accuracy, verification of the contact conditions using the free surface method is conducted during the mechanical equilibrium computations also.

2.4 Mathematical formulation and Strain Increments

Using the well-known mathematical formulations of the non-linear deformation process [8,24], the weak-form of the mechanical equilibrium using a penalty formulation for the incompressibility constraint can be expressed in a compact form as, $${\displaystyle {\int }_{\Omega }\sigma }:{\epsilon }^{*}d\Omega +{\displaystyle {\int }_{\Omega }K}{\epsilon }_v{\epsilon }_v^{*}d\Omega -{\displaystyle {\int }_{\Omega }\rho }b.u^{*}d\Omega =0$$(3)

and the thermal equilibrium can be expressed as, $$\begin{array}{l}{\displaystyle {\int }_{\Omega }k}\nabla T\nabla T^{*}d\Omega +{\displaystyle {\int }_{\Omega }\rho }c\dot{T}{T}^{*}d\Omega \\ -{\displaystyle {\int }_{\Omega }r}{T}^{*}d\Omega -{\displaystyle {\int }_{\partial \Omega }k}\nabla T{T}^{*}d\Gamma =0.\end{array}$$(4)

By using the description of motion in [24], the deformation gradient tensor can be defined as, $$F=\frac{\partial x}{\partial X}={\nabla }_{0}x=\left[I+\frac{\partial u}{\partial X}\right]$$(5)

or using the stretch tensors, the polar decomposition of F as, $$F=RU=VR$$(6)

and the inverse deformation gradient tensor can be defined as, $$F^{-1}=\frac{\partial X}{\partial x}=\nabla X=\left[I-\frac{\partial u}{\partial x}\right].$$(7)

The left and right Cauchy-Green tensors can then be defined as, $$B=F{F}^T=VR{R}^T{V}^T=V^2$$(8) $$C=F^{T}F=U^T{R}^{T}RU=U^2.$$(9)

Hence, the inverse of the left Cauchy-Green tensor gives that, $$B^{-1}=F^{-\text{T}}{F}^{-1}$$(10)

By doing a spectral decomposition, equation (10) can be expressed as, $$B^{-1}=M\left[\begin{array}{ccc}{\lambda }_1^{-2}& 0& 0\\ 0& {\lambda }_2^{-2}& 0\\ 0& 0& {\lambda }_3^{-2}\end{array}\right]M^{\text{T}}$$(11)

The Seth-Hill family of generalized strains [25] defines the logarithmic strain increment in a time-step as, $$\Delta \epsilon =\text{ln}(V).$$(12)

Hence, the principal logarithmic strains can be given by, $$\Delta {\epsilon }_{princ}=\left\{\begin{array}{c}\Delta {\epsilon }_1\\ \Delta {\epsilon }_2\\ \Delta {\epsilon }_3\end{array}\right\}=\left\{\begin{array}{c}\text{ln}({\lambda }_1)\\ \text{ln}({\lambda }_2)\\ \text{ln}({\lambda }_3)\end{array}\right\}$$(13)

and the global logarithmic strains can be given by, $$\Delta \epsilon =M\Delta {\epsilon }_{princ}{M}^{\text{T}}$$(14)

Following the principle of logarithmic strains, the total strain increment can be split into an elastic, viscoplastic and thermal component as, $$\Delta \epsilon =\Delta {\epsilon }_e+\Delta {\epsilon }_{vp}+\Delta {\epsilon }_{th}$$(15)

where, $$\Delta {\epsilon }_{th}=\text{ln}(1+{\alpha }_{th}\Delta T)I.$$(16)

2.5 Stress integration

Since an Eulerian-type formulation of the equilibrium equations are used, the equations need to be solved in the current deformed conditions. Since boundary conditions in the current deformed condition is readily known, the difficulty lies in computing the Cauchy-stress in this current deformed configuration. To have an appropriate description of the state, a proper stress integration needs to be done. Using an incremental approach for the stresses, the stress at the current deformed configuration at time t + Δt can be found as, $${\sigma }^{t+\Delta t}={\sigma }_{t+\Delta t}^t+\Delta \sigma .$$(17)

In the case of infinitesimally small deformation, we can consider that ${\sigma }_{t+\Delta t}^t\simeq {\sigma }^t$, but the same is not valid in the case of finite deformation. Appropriate methods need to be used for evaluating ${\sigma }_{t+\Delta t}^t$. Using the definition of the Piola-Kirchhoff stress tensor and the push-back and pull-forward operators, we can see that, $${\sigma }_{t+\Delta t}^t=\frac{1}{det(F)}F{\sigma }^t{F}^{\text{T}}.$$(18)

Using the elasticity tensor H^e, the stress increment can be found out as, $$\Delta \sigma =H^e\Delta {\epsilon }_e=H^e(\Delta \epsilon -\Delta {\epsilon }_{vp}-\Delta {\epsilon }_{th}).$$(19)

From the consistency model for viscoplasticity [26–28], the yield function for a rate-dependent material considering isotropic constitutive law and Von-Mises plasticity can be written as, $$f(\sigma ,{\overline{\epsilon }}_{vp},{\stackrel{.}{\overline{\epsilon }}}_{vp},T)=\overline{\sigma }-{\sigma }_Y({\overline{\epsilon }}_{vp},{\dot{\overline{\epsilon }}}_{vp},T)=0.$$(20)

The consistency model also allows us to define the viscoplastic strain increment as, $$\Delta {\epsilon }_{vp}=\Delta \lambda \frac{\partial f}{\partial \sigma }.$$(21)

Thus, from equations (19), (16) and (21), equation (17) can be expressed as, $${\sigma }^{t+\Delta t}\text{ \&}={\sigma }_{t+\Delta t}^t+H^e(△\epsilon -\frac{△\lambda }{{\overline{\sigma }}^{t+\Delta t}}P{\sigma }^{t+\Delta t}-\text{ln}(1+{\alpha }_{th}\Delta T)I)={\left(I+\frac{△\lambda }{{\overline{\sigma }}^{t+\Delta t}}{H}^{e}P\right)}^{-1}\times \left({\sigma }_{t+\Delta t}^t+H^e(△\epsilon -\text{ln}(1+{\alpha }_{th}\Delta T)I)\right)$$(22)

Equation (22) can be solved efficiently using the Return Mapping Algorithm [29] or by a Direct Scalar Algorithm [30,31] extended to the thermo-viscoplastic regime.

3 Simulation results

In order to validate the performance of the proposed PIA in thermo-viscoplastic regimes, comparison has been made with results obtained from reference simulations. For ensuring a none-to-one comparison with the reference results, target final meshes in PIA is taken to be the deformed meshes obtained from the reference simulation. A Johnson-Cook viscoplastic model [32] is employed in the examples considered.

3.1 2D axisymmetric screw

The target shape and the initial shape of the workpiece and the tools used can be seen in Figure 1.

For the purposes of the simulation, the workpiece is discretised using 878 nodes and 1578 elements (Type CAX3T, Abaqus). The tools are considered to be rigid and the contact between the tools and the workpiece are assumed to be friction-less.

The workpiece is considered to be made of 12Kh18N10T steel with the following properties [33,34].

Density, ρ = 7930 kg m⁻³, Young’s Modulus, E = 195 GPa, Poisson’s ratio, ν = 0.28, Specific heat capacity, c = 462 J kg⁻¹ K⁻¹, Conductivity, k = 25 Wm⁻¹ K⁻¹, Melting temperature, T_m = 1300 °C, initial temperatureof the workpiece, T = 1000 °C. The workpiece material is considered to follow a Johnson-Cook viscoplastic model with parameters A = 196 MPa, B = 615.5 MPa, n = 0.7005, C = 0.04071, m = 1.479 and ${\dot{\overline{\epsilon }}}_0=0.0001$ s⁻¹. The tools are considered to be iso-thermal with a constant temperature of T_d =750 °C while the ambient temperatureT_a is considered to be 20 °C. The workpiece material is considered to have a heat transfer coefficient h_wd = 310 Wm⁻²K⁻¹ with respect to the tools and h_wa = 50 Wm⁻²K⁻¹ with respectto the ambient environment. Radiation effects and inelastic heat fraction is not taken into consideration for the initial analysis.

The punch is given a movement of 4.4 mm in the negative Y-direction while the bottom die is considered to be stationary. The process time is considered to be 1 second and the results obtained from ABAQUS^® Explicit and PIA are compared below. For the purpose of having a one-to-one comparison of the results, the target shape used in PIA is the final deformed shape obtained from the reference simulation.

As can be seen from the results presented in Figures 2 and 3, the results obtained from PIA compare very favourably with the reference results. The average percentage error in the results maybe computed by considering that the $error\%=\frac{(Referencevalue-PIAresult)}{Referencevalue}$. The error inthe maximum value of the computed stress is 6.35%, in maximum equivalent plasticity is 11.65% and in maximum temperaturevalue is 0.05%. It is also seen that the PIA captures the evolution of the stress, strain and temperature quite well. This is illustrated in Figure 4 where we compare the evolution of the equivalent plastic strain, Von-Mises Stress and temperature at node 236, which is the node where we have the maximum stress/equivalent plastic strain.

One important effect in the hot forging process is the conversion of plastic work into thermal energy. If we consider a Taylor-Quinney coefficient (inelastic heat fraction) of 0.9, we get results as seen in Figure 5.

As can be immediately inferred from the results presented in Figures 5 and 6, the temperature values predicted from the PIA is much higher than the reference results. A proportionate effect is seen in the predicted stress values also due to the effect of the predicted temperature. Also it is seen that for convergence, a higher number of steps (30 steps) is required, mainly due to the effects of the temperature. The error in the maximum value of the computed stress is 9.27%, in maximum value of equivalent plasticity is 11.25% and in maximum temperature value is 1.12%. This higher error is due to the effect of the large step sizes used in PIA. Since we consider a larger step size, the effect of the plastic work is more prominent in that step and gets diffused throughout the body. To verify this conclusion, the same analysis was made by considering a larger number of steps (80 steps) for the PIA and the results of the same can be seen in Figure 7.

It is immediately evident that the results obtained from PIA with a larger number of steps compare more favourably with the reference results and confirms the hypothesis that the difference between the reference results and the results obtained from PIA with a much lower number of steps in Figure 6 is indeed due to the large step size effect. The error in the maximum value of the computed stress is 0.74%, in maximum equivalent plasticity is 3.28% and in maximum temperature value is 0.23%.

The effect of the step sizes is better captured by considering the evolution of the field values at a node for the different cases. We consider 3 different case for the PIA including 15 steps, 30 steps and 80 steps and the evolution of the field variables is plotted in Figure 8.

As can be seen from the plots, the analysis using 15 steps in PIA does not get completed as the simulation crashes at around 93% due to the predicted temperature exceeding the melting temperature. The simulation with 30 steps gets completed but the temperaturefield is over-predicted due to the large step-size effect. This in-turn affects the prediction of the stress field and consequently we see a lower stress when using just 30 steps. The simulation with 80 steps in PIA mitigates the step-size effect and hence gives results that have a very good co-relation with that of the reference results. The equivalent plastic strain predicted by PIA in all 3 cases remain pretty much constant. This underlines the fact that the step-size effects due to the use of large steps is limited to the temperature field prediction and the consequent effect on the stress field prediction. This effect will not play a role in cold-forging simulations using PIA.

Fig. 1 Forging of a screw: (a) target shape, (b) initial billet and tools.

Fig. 2 Reference results obtained from ABAQUS®Explicit- Axisymmetric Screw.

Fig. 3 Results obtained from PIA – Axisymmetric Screw (15 steps).

Fig. 4 Evolution of field values at Node 236 – Axisymmetric Screw. (a) Equivalent Plastic Strain; (b) Von-Mises Stress; (c) Temperature.

Fig. 5 Reference results obtained from ABAQUS®Explicit- Axisymmetric screw with Inelastic heat fraction.

Fig. 6 Results obtained from PIA – Axisymmetric screw with Inelastic heat fraction (30 steps).

Fig. 7 Results from 80 steps PIA – Axisymmetric Screw with Inelastic heat fraction.

Fig. 8 Evolution of field values at Node 236 – Axisymmetric Screw including inelastic work. (a) Equivalent Plastic Strain, (b) Von-Mises Stress, (c) Temperature.

3.2 2D axisymmetric wheel

As can be seen with the model in Section 3.1 there are differences in the temperature and stress fields when the inelastic heat fraction is taken into account. The effect of contact and the large step sizes cannot be easily identified in this example since the area where an evolution of contact happens is relatively small. To investigate the effect of the contact evolution and the subsequent changes in the solution, another example of an axisymmetric wheel is taken. The target shape of the wheel is shown in Figure 9. The initial shape of the workpiece used and the tools used for the process can be found in Figure 10.

For the purposes of the study, the workpiece is discretised by 3123 nodes and 5999 elements (Type CAX3T, Abaqus). The tools are considered to be rigid and friction-less contact is assumed on the contact surfaces. The workpiece is considered to be made of aluminium with the following properties [35].

Density, ρ = 2810 kg m⁻³, Young’s Modulus, E = 71.7 GPa, Poisson’s ratio, ν = 0.33, Specific heat capacity, c = 960 J kg⁻¹ K⁻¹, Conductivity, k = 205 Wm⁻¹ K⁻¹, Melting temperature, T_m = 477 °C, initial temperatureof the workpiece, T = 400 °C. The workpiece material is considered to follow a Johnson-Cook viscoplastic model with parameters A = 517 MPa, B = 405 MPa, n = 0.41, C = 0.0075, m = 1.1 and ${\dot{\overline{\epsilon }}}_0=0.000161$ s⁻¹. The tools are considered to be iso-thermal with a constant temperature of T_d =250 °C while the ambient temperatureT_a is considered to be 25 °C. The workpiece material is considered to have a heat transfer coefficient h_wd = 310 Wm⁻² K⁻¹ with respect to the tools, h_wa = 50 Wm⁻² K⁻¹ with respectto the ambient environment.

The punch is given a movement of 23.7 mm in the negative Y-direction while the bottom die is considered to be stationary and the process time is considered to be 1 second and Taylor-Quinney coefficient for the inelastic heat fraction κ = 0.9.

Figures 11–13 show the results obtained from the reference simulation and Figures 14–16 show the results obtained from PIA. A comparison of the reference results and the results obtained by PIA show that the results obtained from PIA are very close to those obtained in the reference simulation.

The average error in the computed Von-Mises stress is 2.11%, in the computed Equivalent Plastic Strain is 10.99% and in the temperaturevalue is 0.11%. As is quite evident from the comparison plots, the PIA is able to provide a good approximation of the field values accounting for the step-size effects.

Fig. 9 Target shape of the wheel.

Fig. 10 Tools and workpiece used for forging the 2D axisymmetric wheel.

Fig. 11 Reference results – Axisymmetric Wheel: Equivalent Plastic Strain.

Fig. 12 Reference results – Axisymmetric Wheel: Von-Mises Stress.

Fig. 13 Reference results – Axisymmetric Wheel: Temperature.

Fig. 14 PIA results 25 steps – Axisymmetric Wheel: Equivalent Plastic Strain.

Fig. 15 PIA results 25 steps – Axisymmetric Wheel: Von-Mises Stress.

Fig. 16 PIA results 25 steps – Axisymmetric Wheel: Temperature.

3.2.1 Effect of time steps and inelastic heat fraction

To study the effects of the number of steps and the inelastic heat fraction in the prediction of the temperature and the consequent changes in the predicted Von-Mises stress and the equivalent plastic strains, the evolution of the field values at a node inside the element is done for different load cases as detailed in Table 1.

The field values at node number 965 is extracted since it is the node with the highest temperature values in almost all of the case studied. The results from the Reference simulation and PIA is presented in Figures 17–19. The solid lines represent the solution obtained from the reference simulations while the dotted lines represent the solution obtained from PIA.

From Figure 17 it is apparent that the equivalent plasticity values predicted by both the reference simulation and PIA arevery similar except in the case when the process time is very small. This deviation could also be due to the effect of the explicit simulation employed for the reference results. The biggest differences between the results obtained from PIA and the reference simulations can be observed in the Von-Mises Stress values (Fig. 18) and the Temperature values (Fig. 19). These plots clearly demonstrate the effect of the large steps used in PIA. The conductive heat transfer has a highereffect than the convective heat transfer in the example that we have studied. Since PIA uses large step sizes, the elements on the surface are considered to be in contact for a larger time interval than when a smaller step-size is used. Hence this reflects as a lower temperature at the node. This effect is most prominently observed in the results obtained for case 4 where there is no volumetric heat source considered and there is only heat loss from the billet to the surrounding environment and the tools. Another consequence of the large step sizes is the larger effect of the inelastic heat fraction which has observed in the previous example of the axisymmetric screw as well. In the present example, this effect can be clearly observed in Case 1 and Case 2, while in Case 3 the effect of the inelasticheat fraction is overshadowed by the effect of heat loss due to thermal conductance at the contact interfaces.

The accuracy of the PIA with respect to the different cases that has been studied can be seen in Table 2. It can be noted that the accuracy of prediction of the equivalent plasticity has a lower effect compared to the accuracy of the prediction of the temperaturefields. Of particular significance is the fact that in the first three cases where the average error in the temperature fields is ~0.4%, the average error in the computed stress is only around 2–3%. But, when the average error in the predicted temperature fields rises to ~1%, the average error in the predicted Von-Mises stress values rises sharply to ~5%. Thus it can be concluded from both the axisymmetric examples that have been considered that, even though PIA is able to predict the deformed shapes well, in the case of models where the effect of temperature on the material properties is high, proper care must be taken in terms of the step sizes to capture the thermal effect and the consequent effects in the stress.

Fig. 17 Evolution of field values at node 965 – Equivalent Plastic Strain.

Fig. 18 Evolution of field values at node 965 – Von-Mises Stress.

Fig. 19 Evolution of field values at node 965 – Temperature.

3.3 Computation time

The PIA code that has been developed during the course of this work is predominantly an academic effort and hence it is not optimized version that can be directly used to compare with commercially available software that has been completely optimized. Moreover, the reference results used in this work are obtained from ABAQUS^® Explicit which uses an explicit formulation compared to an implicit formulation used in PIA. Nevertheless, the computation times of both the methods for the different cases that we have studied in the previous section is compared in Figure 20. The computation time corresponds to the single-core performance using an Intel^® Core™2 Quad CPU Q9300 @ 2.5 GHz with 8 GB system memory and full double precision. Usually in an explicit simulation, the stable step-size increments is calculated based on the element sizes and hence can cause very high computation times in cases where the process time is high. Mass-scaling is a common technique used in practice to keep computation times low and is used in the reference simulations to keep computation times low. A semi-automatic mass scaling is used in our case so that the target time increment in each step is 1e-06.

From the results it is clear that PIA gives a clear advantage over classical incremental methods, but for the proper quantification and comparison of the computational time advantage of using PIA, a comparison should be made with a commercially available software based on the implicit analysis like FORGE^®NxT. Due to lack of resources, the same has not been attempted in this work.

Fig. 20 Computation time – with mass scaling.

4 Conclusion

From the examples that we have studied, we can conclusively show that an accurate modelling of thermal effects is of paramount importance in hot-forging simulations. The core objective of this paper has been to demonstrate that the PIA, previously developed in the context of cold-forging, can be successfully extended to handle the thermo-viscoplastic multi-physiscs that is a characteristic of hot-forging and hence, enhance its suitability to be employed in industrial cases. From the results obtained, we can see that PIA is capable of predicting the temperature fields and the stress/strain fields with a high degree of accuracy at a highly reduced computational effort. However, it is to be noted that the chosen step-sizes do play an important part in the accuracy of the obtained results. Hence, automatic time-stepping strategies that can ensure a good balance between efficiency and accuracy needs to be developed to fully realize the potential of PIA. Nevertheless, PIA promises to be an exciting frontier that can lead to practical optimization simulations of the hot and cold forging processes at a much reduced computational cost in comparison with classical approaches.

Nomenclature

u: Displacement

T: Temperature

σ: Cauchy stress tensor

σ^t: Cauchy stress tensor at time t in configuration Ω^t

${\sigma }_{t+\Delta t}^t$: Cauchy stress tensor at time t in configuration Ω^t+Δt

Δε: Logarithmic strain increment

ε^*: Virtual strain tensor

ε_v: Volumetric strain

K: Penalty coefficient

ρ: Density

b: Volumetric body force

k: Coefficient of conductivity

c: Coefficient of specific heat

r: Volumetric body heat source

F: Deformation gradient tensor

R: Orthogonal Rotation tensor

U: Right stretch tensor

V: Left stretch tensor

C: Right Cauchy-Green stress tensor

B: Left Cauchy-Green stress tensor

λ_i: Principal Eigen values

Δε_princ: Strain increment in principal directions

I: Identity matrix

σ_Y: Yield stress

$\overline{\sigma }$: Equivalent stress

$\overline{\epsilon }$: Equivalent strain

$\dot{\overline{\epsilon }}$: Equivalent strain rate

α_th: Coefficient of linear thermal expansion

Δλ: Plastic increment

References

All Tables

All Figures

	Fig. 1 Forging of a screw: (a) target shape, (b) initial billet and tools.
In the text

	Fig. 2 Reference results obtained from ABAQUS®Explicit- Axisymmetric Screw.
In the text

	Fig. 3 Results obtained from PIA – Axisymmetric Screw (15 steps).
In the text

	Fig. 4 Evolution of field values at Node 236 – Axisymmetric Screw. (a) Equivalent Plastic Strain; (b) Von-Mises Stress; (c) Temperature.
In the text

	Fig. 5 Reference results obtained from ABAQUS®Explicit- Axisymmetric screw with Inelastic heat fraction.
In the text

	Fig. 6 Results obtained from PIA – Axisymmetric screw with Inelastic heat fraction (30 steps).
In the text

	Fig. 7 Results from 80 steps PIA – Axisymmetric Screw with Inelastic heat fraction.
In the text

	Fig. 8 Evolution of field values at Node 236 – Axisymmetric Screw including inelastic work. (a) Equivalent Plastic Strain, (b) Von-Mises Stress, (c) Temperature.
In the text

	Fig. 9 Target shape of the wheel.
In the text

	Fig. 10 Tools and workpiece used for forging the 2D axisymmetric wheel.
In the text

	Fig. 11 Reference results – Axisymmetric Wheel: Equivalent Plastic Strain.
In the text

	Fig. 12 Reference results – Axisymmetric Wheel: Von-Mises Stress.
In the text

	Fig. 13 Reference results – Axisymmetric Wheel: Temperature.
In the text

	Fig. 14 PIA results 25 steps – Axisymmetric Wheel: Equivalent Plastic Strain.
In the text

	Fig. 15 PIA results 25 steps – Axisymmetric Wheel: Von-Mises Stress.
In the text

	Fig. 16 PIA results 25 steps – Axisymmetric Wheel: Temperature.
In the text

	Fig. 17 Evolution of field values at node 965 – Equivalent Plastic Strain.
In the text

	Fig. 18 Evolution of field values at node 965 – Von-Mises Stress.
In the text

	Fig. 19 Evolution of field values at node 965 – Temperature.
In the text

	Fig. 20 Computation time – with mass scaling.
In the text