Issue 
Mechanics & Industry
Volume 21, Number 2, 2020



Article Number  202  
Number of page(s)  12  
DOI  https://doi.org/10.1051/meca/2020002  
Published online  05 February 2020 
Regular Article
Thermoviscoplastic numerical modeling of metal forging process by Pseudo Inverse Approach
^{1}
GRESPI, Université de Reims ChampagneArdenne, UFR Sciences, Campus Moulin de la Housse,
51100
Reims, France
^{2}
ESI Group, Parc ICADE,
3 bis Rue Saarinen,
94528
Rungis CEDEX, France
^{*} email: boussad.abbes@univreims.fr
Received:
3
June
2019
Accepted:
2
January
2020
Constant demands of lightweighting have forced many industries to resort to manufacturing practices that promise a component with a higher strengthtoweight ratio. Hot forging is one such method used to produce parts using difficulttoform materials as well as to achieve complex geometries. Although numerical methods provide an efficient means to predict the material yield and the stress/strain states of the product at different stages of forming and classical methods are accurate enough to provide a suitable representation of the process, they tend to be computationally expensive. This limits their use in process optimization studies. The Pseudo Inverse Approach (PIA) developed in the context of 2D axisymmetric cold forming, provides a quick and fairly accurate estimate of the stress and strain fields in the final product for a given initial shape. In this work, the PIA is extended to include the thermal and viscoplastic effects associated with the hot forging process. The results are compared with commercially available software based on the classical approaches to show the efficiency and the limitations of PIA.
Key words: Pseudo Inverse Approach / hot forging / thermoviscoplastic / finite element method
© 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 workpieces [2]. A series of developments in the large deformation elastoplastic 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 stepbystep 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 nonlinear in nature and as a consequence the direct analysis requires large computational effort.
As the major objective in metalforming 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 workpiece 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], parameterbased Bspline optimization [8,9], parameter free shape optimization [10], metamodel 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 illposed 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 singlestep 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 coldforging 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 multistep 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 metalforming 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 workpiece. Since materials at elevated temperatures are usually rate sensitive, a complete analysis of hot forming operations requires a consideration of the effect of ratesensitivity (viscoplasticity) of materials and the coupling of the metalflow and thermal analyses [21,22]. The present work extends the PIA for the numerical modelling of hotforging 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 Euleriantype 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 inverseshapefinding 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 closeenough deformed configuration that can be used in a largestep Euleriantype 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 thumbrules used by forging process designers.
 3.
Transfer the mesh created on the target deformed shape to the generated initial shape using meshmorphing 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 timestep, 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 nonlinear 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 timerange 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 singlestep 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^{0} of the initialshape. This can be efficiently done using mesh morphing algorithms in use. The interior nodes of the mesh M^{0} 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, (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 nodetoelement 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 timestep 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 freesurface method. In a metalforming operation simulation, for any configuration Ω, the freesurface condition can be expressed as, (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_{1}, u_{2} and u_{3} 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_{1} ≠ 0; u_{2}≠0; u_{3}≠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_{1} = 0; u_{2}≠0; u_{3}≠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 wellknown mathematical formulations of the nonlinear deformation process [8,24], the weakform of the mechanical equilibrium using a penalty formulation for the incompressibility constraint can be expressed in a compact form as, (3)
and the thermal equilibrium can be expressed as, (4)
By using the description of motion in [24], the deformation gradient tensor can be defined as, (5)
or using the stretch tensors, the polar decomposition of F as, (6)
and the inverse deformation gradient tensor can be defined as, (7)
The left and right CauchyGreen tensors can then be defined as, (8) (9)
Hence, the inverse of the left CauchyGreen tensor gives that, (10)
By doing a spectral decomposition, equation (10) can be expressed as, (11)
The SethHill family of generalized strains [25] defines the logarithmic strain increment in a timestep as, (12)
Hence, the principal logarithmic strains can be given by, (13)
and the global logarithmic strains can be given by, (14)
Following the principle of logarithmic strains, the total strain increment can be split into an elastic, viscoplastic and thermal component as, (15)
2.5 Stress integration
Since an Euleriantype 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 Cauchystress 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, (17)
In the case of infinitesimally small deformation, we can consider that , but the same is not valid in the case of finite deformation. Appropriate methods need to be used for evaluating . Using the definition of the PiolaKirchhoff stress tensor and the pushback and pullforward operators, we can see that, (18)
Using the elasticity tensor H^{e}, the stress increment can be found out as, (19)
From the consistency model for viscoplasticity [26–28], the yield function for a ratedependent material considering isotropic constitutive law and VonMises plasticity can be written as, (20)
The consistency model also allows us to define the viscoplastic strain increment as, (21)
Thus, from equations (19), (16) and (21), equation (17) can be expressed as, (22)
Equation (22) can be solved efficiently using the Return Mapping Algorithm [29] or by a Direct Scalar Algorithm [30,31] extended to the thermoviscoplastic regime.
3 Simulation results
In order to validate the performance of the proposed PIA in thermoviscoplastic regimes, comparison has been made with results obtained from reference simulations. For ensuring a nonetoone comparison with the reference results, target final meshes in PIA is taken to be the deformed meshes obtained from the reference simulation. A JohnsonCook 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 frictionless.
The workpiece is considered to be made of 12Kh18N10T steel with the following properties [33,34].
Density, ρ = 7930 kg m^{−3}, Young’s Modulus, E = 195 GPa, Poisson’s ratio, ν = 0.28, Specific heat capacity, c = 462 J kg^{−1} K^{−1}, Conductivity, k = 25 Wm^{−1} K^{−1}, Melting temperature, T_{m} = 1300 °C, initial temperatureof the workpiece, T = 1000 °C. The workpiece material is considered to follow a JohnsonCook viscoplastic model with parameters A = 196 MPa, B = 615.5 MPa, n = 0.7005, C = 0.04071, m = 1.479 and s^{−1}. The tools are considered to be isothermal 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^{−2}K^{−1} with respect to the tools and h_{wa} = 50 Wm^{−2}K^{−1} 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 Ydirection 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 onetoone 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 . 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, VonMises 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 TaylorQuinney 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 overpredicted due to the large stepsize effect. This inturn 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 stepsize effect and hence gives results that have a very good corelation 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 stepsize 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 coldforging 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) VonMises 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) VonMises 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 frictionless 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^{−3}, Young’s Modulus, E = 71.7 GPa, Poisson’s ratio, ν = 0.33, Specific heat capacity, c = 960 J kg^{−1} K^{−1}, Conductivity, k = 205 Wm^{−1} K^{−1}, Melting temperature, T_{m} = 477 °C, initial temperatureof the workpiece, T = 400 °C. The workpiece material is considered to follow a JohnsonCook viscoplastic model with parameters A = 517 MPa, B = 405 MPa, n = 0.41, C = 0.0075, m = 1.1 and s^{−1}. The tools are considered to be isothermal 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^{−2} K^{−1} with respect to the tools, h_{wa} = 50 Wm^{−2} K^{−1} with respectto the ambient environment.
The punch is given a movement of 23.7 mm in the negative Ydirection while the bottom die is considered to be stationary and the process time is considered to be 1 second and TaylorQuinney 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 VonMises 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 stepsize 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: VonMises 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: VonMises 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 VonMises 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 VonMises 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 stepsize 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 VonMises 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.
Test cases – specifications.
Fig. 17
Evolution of field values at node 965 – Equivalent Plastic Strain. 
Fig. 18
Evolution of field values at node 965 – VonMises Stress. 
Fig. 19
Evolution of field values at node 965 – Temperature. 
Accuracy of PIA v/s Reference results.
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 singlecore 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 stepsize increments is calculated based on the element sizes and hence can cause very high computation times in cases where the process time is high. Massscaling 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 semiautomatic mass scaling is used in our case so that the target time increment in each step is 1e06.
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 hotforging simulations. The core objective of this paper has been to demonstrate that the PIA, previously developed in the context of coldforging, can be successfully extended to handle the thermoviscoplastic multiphysiscs that is a characteristic of hotforging 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 stepsizes do play an important part in the accuracy of the obtained results. Hence, automatic timestepping 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
σ^{t}: Cauchy stress tensor at time t in configuration Ω^{t}
: Cauchy stress tensor at time t in configuration Ω^{t+Δt}
Δε: Logarithmic strain increment
k: Coefficient of conductivity
c: Coefficient of specific heat
r: Volumetric body heat source
F: Deformation gradient tensor
C: Right CauchyGreen stress tensor
B: Left CauchyGreen stress tensor
Δε_{princ}: Strain increment in principal directions
α_{th}: Coefficient of linear thermal expansion
References
 T. Altan, G. Ngaile, G. Shen, eds., Cold And Hot Forging: Fundamentals And Applications. ASM International, 2004. [Google Scholar]
 J.L. Chenot, M. Bernacki, P.O. Bouchard, L. Fourment, E. Hachem, E. Perchat, in 11th International Conference on Technology of Plasticity, ICTP 2014, edited by M.K.I. Ishikawa (Nagoya, Japan, 2014), vol. 81 [Google Scholar]
 J.H. Cheng, N. Kikuchi, An analysis of metal forming processes using large deformation elasticplastic formulations, Comput. Methods Appl. Mech. Eng. 49, 71–108 (1985) [Google Scholar]
 C. Bohatier, J.L. Chenot, Finite element formulations for nonsteadystate large viscoplastic deformation, Int. J. Numer. Methods Eng. 21, 1697–1708 (1985) [Google Scholar]
 P. Hartley, I. Pillinger, Numerical simulation of the forging process, Comput. Methods Appl. Mech. Eng. 195, 6676–6690 (2006) [Google Scholar]
 J. Mackerle, Finite element modelling and simulation of bulk material forming: a bibliography (1996–2005), Eng. Comput. 23, 250–342 (2006) [Google Scholar]
 S. Badrinarayanan, N. Zabaras, A sensitivity analysis for the optimal design of metalforming processes, Comput. Methods Appl. Mech. Eng. 129, 319–348 (1996) [Google Scholar]
 L. Fourment, J.L. Chenot, Optimal design for nonsteadystate metal forming processes I. shape optimization method, Int. J. Numer. Methods Eng. 39, 33–50 (1996) [Google Scholar]
 L. Fourment, J.L. Chenot, T. Balan, Optimal design for nonsteadystate metal forming processes ii. application of shape optimization in forging, Int. J. Numer. Methods Eng. 39, 51–65 (1996) [Google Scholar]
 P. Steinmann, P. Landkammer, A fast approach to shape optimization using the inverse FEM, in Material Forming ESAFORM 2014, vol. 611 of Key Engineering Materials. Trans Tech Publications (2014) 1404–1410 [Google Scholar]
 M. Ejday, L. Fourment, Metamodel assisted evolutionary algorithm for multiobjective optimization of nonsteady metal forming problems, Int. J. Mater. Form. 3, 5–8 (2010) [CrossRef] [Google Scholar]
 C. Ciancio, T. Citrea, G. Ambrogio, L. Filice, R. Musmanno, Design of a high performance predictive tool for forging operation, Proc. CIRP 33, 173–178 (2015) [CrossRef] [Google Scholar]
 D.M. D’Addona, D. Antonelli, Neural network multiobjective optimization of hot forging, Proc. CIRP 67, 498–503 (2018) [CrossRef] [Google Scholar]
 J. Park, N. Rebelo, S. Kobayashi, A new approach to preform design in metal forming with the finite element method, Int. J. Mach. Tool Des. Res. 23, 71–79 (1983) [CrossRef] [Google Scholar]
 S. Germain, P. Landkammer, P. Steinmann, On a recursive formulation for solving inverse form finding problems in isotropic elastoplasticity, Advanced Modeling and Simulation in Engineering Sciences (2014), vol. 1, p. 10 [Google Scholar]
 P. Landkammer, S. Germain, P. Steinmann, On inverse form finding for orthotropic plasticity, Comput. Assist. Mech. Eng. Sci. 20, 337–348 (2014) [Google Scholar]
 Y.Q. Guo, J.L. Batoz, J.M. Detraux, P. Duroux, Finite element procedures for strain estimations of sheet metal forming parts, Int. J. Numer. Methods Eng. 30, 1385–1401 (1990) [Google Scholar]
 A. Halouani, Y. Li, B. Abbes, Y.Q. Guo, An efficient axisymmetric element based on inverse approach for cold forging modeling, Eng. Lett. 18 (2010) [Google Scholar]
 W. Gati, Y.Q. Guo, H. Naceur, J.L. Batoz, Approche pseudo inverse pour estimation des contraintes dans les pieces embouties axisymetriques, Revue Eur. Elements Finis 12, 863–886 (2003) [Google Scholar]
 A. Halouani, Y.M. Li, B. Abbes, Y.Q. Guo, Simulation of axisymmetrical cold forging process by efficient pseudo inverse approach and direct algorithm of plasticity, Finite Elements Anal. Des. 61, 85–96 (2012) [CrossRef] [Google Scholar]
 N. Rebelo, S. Kobayashi, A coupled analysis of viscoplastic deformation and heat transfer. I: Theoretical considerations. Int. J. Mech. Sci. 22, 699–705 (1980) [CrossRef] [Google Scholar]
 N. Rebelo, S. Kobayashi, A coupled analysis of viscoplastic deformation and heat transfer. II: Applications. Int. J. Mech. Sci. 22, 707–718 (1980) [CrossRef] [Google Scholar]
 A. Halouani, Y.M. Li, B. Abbès, Y.Q. Guo, A pseudo inverse approach with kinematically admissible intermediate configurations for the axisymmetrical cold forging modelling, Adv. Mater. Res. 399, 1832–7 (2011) [CrossRef] [Google Scholar]
 J. Bonet, R.D. Wood, Nonlinear Continuum Mechanics for Finite Element Analysis. Cambridge University Press (2008) [CrossRef] [Google Scholar]
 R. Hill, On constitutive inequalities for simple materials—i, J. Mech. Phys. Solids 16, 229–242 (1968) [Google Scholar]
 W.M. Wang, L.J. Sluys, R. de Borst, Viscoplasticity for instabilities due to strain softening and strainrate softening, Int. J. Numer. Methods Eng. 40, 3839–3864 (1997) [Google Scholar]
 W. Wang, L. Sluys, Formulation of an implicit algorithm for finite deformation viscoplasticity, Int. J. Solids Struct. 37, 7329–7348 (2000) [Google Scholar]
 A. Carosio, K. Willam, G. Etse, On the consistency of viscoplastic formulations, Int. J. Solids Struct. 37, 7349–7369 (2000) [Google Scholar]
 J.C. Simo, R.L. Taylor, A return mapping algorithm for plane stress elastoplasticity, Int. J. Numer. Methods Eng. 22, 649–670 (1986) [Google Scholar]
 Y.M. Li, B. Abbes, Y.Q. Guo, Two efficient algorithms of plastic integration for sheet forming modeling, J. Manufactur. Sci. Eng. 129, 698 (2007) [CrossRef] [Google Scholar]
 B. Abbes, Y.M. Li, A. Halouani, Y.Q. Guo, Fast plastic integration algorithms for forming process simulations, in Material Forming ESAFORM 2014, vol. 611 of Key Engineering Materials. Trans Tech Publications (2014) 1336–1343 [Google Scholar]
 G. Johnson, W. Cook, A constitutive model and data for metals subjected to large strains, high strain rates, and high temperatures, in Proceedings of the 7th International Symposium on Ballistics, The Hague (1983) 541–547 [Google Scholar]
 S.V. Stankus, I.V. Savchenko, A.V. Baginskii, O.I. Verba, A.M. Prokop’ev, R.A. Khairulin, Thermal conductivity and thermal diffusivity coefficients of 12kh18n10t stainless steel in a wide temperature range, High Temperature 46, 731–733 (2008) [CrossRef] [Google Scholar]
 A. Sobolev, M. Radchenko, Use of johnson–cook plasticity model for numerical simulations of the SNF shipping cask drop tests, Nucl. Energy Technol. 2, 272–276 (2016) [CrossRef] [Google Scholar]
 E. Corona, G.E. Orient, An evaluation of the johnsoncook model to simulate puncture of 7075 aluminum plates, tech. rep., Sandia National Laboratories, 2014 [Google Scholar]
Cite this article as: A.E. Thomas, B. Abbès, Y.M. Li, F. Abbès, J.L. Duval, Thermoviscoplastic numerical modeling of metal forging process by Pseudo Inverse Approach, Mechanics & Industry 21, 202 (2020)
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) VonMises 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) VonMises 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: VonMises 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: VonMises 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 – VonMises 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.