Numerical-Experimental Study Regarding the Single Point Incremental Forming Process

The present paper proposes a numerical-experimental comparative study on the single point incremental forming process. A DC04 steel sheet with a thickness of 0.6 mm was used for both the numerical simulation using the finite element method and the experimental research. The type of trajectory used was a spiral trajectory and the finished part obtained was a truncated cone-shaped part. The analysis program used for simulation was Ls-Dyna. The simulations were performed in several variants: with a fixed mesh and with an adaptive mesh, using two different element formulations: 25 (Belytschko-Tsay formulation with thickness stretch) and -16 (fully integrated shell element modified for higher accuracy) and two contact types: automatic surface to surface (ASTS) and forming one way surface to surface (FOSS). The results of the numerical analysis and of the experimental research were focused on determining the major strain, minor strain, thickness reduction and forces at the end of the single point incremental forming process, as well as determining the processing time.


Introduction
Numerical simulation of metal forming processes has become a necessity in recent years when it comes to the analysis of the phenomenon or to the die design. While it has become somewhat common in the case of conventional metal forming processes, where there are programs dedicated to different processes (e.g. deep-drawing), in the case of unconventional deformation processes, such as the single point incremental forming process (SPIF), simulation using the finite element method is still a challenge. This is primarily due to the continuous changing of the contact between the blank sheet and the punch, the different loads to which the material is subjected, the different hardening laws that can be used, as well as the use of an adaptive or fixed mesh.
Duflou et al. and Behera et al. published two review papers in which they dedicated important chapters related to the numerical simulation using the finite element method of the SPIF process [1,2]. They highlighted the types of simulations, the types of elements used, hardening laws, flow criteria, the types of contacts, as well as the simulation programs used.
Rusu et al. applied the simulation process using the finite element method to study the conditions that favour the material failure during single point incremental forming. In the simulation they used the Von Misses flow criteria and the maximum plastic strain from the uniaxial tensile test as the failure criteria [3][4][5].
Shrivastava et al. built a finite element analysis based on the Johnson-Cook material model for aluminium sheets that were preheated prior to being manufactured by SPIF. They used the Altair HyperWorks software to analyse the effect that preheating may have on the forming behaviour during SPIF. Based on the numerical simulation, they analysed the results obtained for the Von Misses stress, thickness distribution and part precision [6].
Gupta and Jeswiet simulated the single point incremental forming process by means of implicit and explicit simulation techniques, but also used shell, solid or solid-shell finite elements. They also used the mass scaling and the sub-structuring of the finite element network [7].
Two other papers analysed the meshing techniques used for the finite element simulation of the SPIF process [8]. Sena et al. used various remeshing techniques for solid-shell elements that combine the advantages of shell formulation with the hexahedral solid topology to increase the accuracy of the analysis and to limit the computational time. The analysis program used was Lagamine, an implicit analysis program.
Husain et al. simulated the material behaviour during the hole flanging by SPIF of rolled AA1060 aluminium sheet. They used a combination of shell and solid elements and analysed the thickness variation along the flange wall [9]. Another paper also analyses the thickness variation but in the case of the SPIF of pyramidal frustum shapes for DC04 steel [10].
Popp et al. also used the numerical simulation to analyse one of the newest applications of the single point incremental forming, namely implant manufacturing. They determined by numerical simulation the variation of forces, stresses and main deformation in the SPIF process [11].
Suresh et al. used the Ls-Dyna software to analyse the multi-stage SPIF process for a conical part with an 85 0 degrees wall angle. They analysed the part accuracy, strain and thickness distribution and compared them to the experimental results. They concluded that with the increase of the complexity of the model and implicitly of the complexity of the punch trajectory, the accuracy of the results decreases [12].
Benedetti et al. conducted a comparative study of experimental and finite element analysis on two different alloys, a micro-alloy steel and an aluminium alloy with completely different mechanical properties and different anisotropy coefficients. The results obtained following the finite element analysis were compared to the experimental results based on digital image correlation in five points on the part [13].
Based on the state-of-the-art analysis in the field of simulation using the finite element method of the SPIF process, the authors in the present paper focused on conducting a comparative study on the use of shell elements with different formulations, the opportunity of using an adaptive or fixed mesh, as well as the opportunity of using different types of contact. The results obtained were compared to the results obtained experimentally.

Materials and method
The purpose of this paper is to determine the optimal analysis variant to obtain the smallest possible percentage error, as well as the shortest possible analysis time. For this reason, the following 3 parameters were used as variables: the mesh type (a fixed mesh and an adaptive mesh), the formulation type for the shell elements (formulation 25 and formulation -16) and the contact type (automatic surface to surface -ASTS and forming one-way surface to surface -FSTS). This resulted in a complex factorial experimental model of type 23, i.e. a total of 8 analyses. Table 1 shows all the analyses performed. The principal strains (both logarithmic strains) and thickness reduction (engineering strain), the variation of forces on the three directions of the coordinate axes, as well as the calculation time, were determined for each analysis. Regarding the mesh type, for the fixed mesh variant, the semi-finished product was firstly discretized in such a way that the length of the side of an element measured about 3 mm, followed by a refinement on the central area, where the punch moves. Thus, the final length of the element reached 1 mm. The discretized geometric model for this case is shown in Figure 1 (initial stage and final stage). Due to the manner in which the elements are connected, the network consisted of a number of 18175 nodes connected in a number of 18005 elements. The case of the adaptive mesh variant also started with a mesh having a large element size, of 4 mm. Following the meshing, a network containing 4096 elements and 4225 nodes was obtained. The discretized geometric model for this case is shown in Figure 2 (initial stage and final stage). The Ls-Dyna analysis program provides the user with two indicators based on which the semi-finished product can be refined: the change of the angle between two adjacent finite elements or the thinning of the part. On account of the specifics of the matter, the second indicator was chosen. Level 2 was chosen as a refined level, which allows the division of each finite element into a maximum of 16 elements. Following the adaptive remeshing, the element reached a final length of about 1 mm also. From the shell formulation point of view, two of the most common formulations were chosen: formulation 25 and formulation -16. Formulation 25 (Belytschko-Tsay with thickness stretch), which has the option of adopting a linear strain through thickness, was selected due to its recognized shorter analysis times. Formulation -16 (fully integrated shell elements), which uses the Bathe-Dvorkin transverse shear treatment, was chosen in view of the fact that it does not strongly distort the shell elements during the simulation.
In terms of the contact types, two of the most frequently used contact types in Ls-Dyna were selected: ASTS requires defining the thickness including that of the rigid parts, while FOSS is easier to model as the thickness definition is not required for the rigid parts and is especially designed for the numerical simulation of the deep-drawing process. Regardless of the contact type, it was materialized using a static coefficient of friction μ = 0.08. The time increment used was 3e-06 seconds.
Most pre-processors that come with software packages allow the analytical or discrete modelling of rigid (non-deformable) bodies. The Marc program can be mentioned among those which use the analytical modelling, while those which use the discrete modelling are: Abaqus, Hyperform and Ls-Dyna. In this case the active elements are also discretely described, not analytically.
Blank sheets were modelled as solid deformable parts, while the punch, die and blank holder were modelled as rigid, non-deformable parts. The blank sheet is placed on the die and is fastened to the blank holder by means of a retaining force of 5000 N, uniformly applied on its surface. The impossibility of the semi-finished product sliding through the clearance between the blank holder and the die is thus ensured.
The die is generated from a number of 1841 elements and 1978 nodes, respectively. The punch was generated in the shape of a spherical surface, positioned in point contact with the blank sheet. The punch contains 599 elements and 602 nodes, respectively. The blank holder is modelled as an annulus area which is composed of 1042 elements connected at 1174 nodes.
The part, meshed as a deformable solid body, is meshed using Thin Shell 163 elements. No boundary conditions which would cancel degrees of freedom were imposed on the nodes of the semi-finished product that are located on the circumference, this being achieved by means of the retaining ring and the retaining force. Nor are other boundary conditions explicitly prescribed, such as imposed nodal displacements or force loads in nodes, but are achieved during the running of the program as a result of the imposing of no-penetration conditions of the nodes of the discretized semi-finished product.
The blank sheets were 250 x 250 mm in size, with a thickness of 0.6 mm. The trajectory used was a frustum of cone trajectory with a large base of 85 mm and the total height of the part was 40 mm (Figure 3). A punch with an 8 mm diameter was used and the vertical step was 0.75 mm.  Since the blank sheet material was a DC04 steel, material 37, which represents a transversely anisotropic elastic-plastic material that uses Hill 48 yield function. was chosen as the type of material associated with its elements. The modulus of elasticity considered is E = 73160 MPa and the transverse contraction ratio ν = 0.28. The yield stress determined by the imposed uniaxial tensile test is σc = 152.6 MPa. The introduction in the program of the true yield stress-true plastic strain curve was chosen in order to define the plastic zone ( Figure  4). The strength coefficient K = 561 MPa and the hardening exponent n = 0.224 were also determined on the basis of the uniaxial tensile test. The value of the anisotropy coefficient is R = 1.99. To reduce the time for analysis, a mass scaling was performed, which led to an increase in mass by about 10%.
Given that the breaking of the part was out of question, the criterion of the occurrence of a break was not taken into consideration.
The Kuka Kr210 Robot was used to conduct the experimental research. The main strains and thinning were determined by the digital image correlation method, using the Aramis optical measurement system.
In order to be able to measure the main strains, the blank sheet is vertically positioned, while the cameras of the optical system are located on the opposite side of the manufacturing area. Before the forming process the blank sheet was painted with a matt white enamel and after the paint dried a mesh of black graphite points was applied to the surface of the blank sheet. The optical system measured the displacement of each black point and calculated the principal strains.
The PCB 261A13 force sensor was used to determine the forces in the process, together with the Quantum X MX840B digital acquisition system and the CMD 600 digital charge amplifier. The force sensor allows the measurement of forces up to 44 kN in the direction perpendicular to the blank sheet plane (z direction) and 19 kN in the directions in the blank sheet plane (x and y directions).
Due to the fact that the forces required for the single point incremental forming of steel sheet are significant values, the sensor was calibrated using the Instron 5587 testing machine in the range of 0 … 3 kN for the direction perpendicular to the blank sheet plane (z direction), and 0 … 1.5 kN for the directions in the blank sheet plane (x and y direction), respectively. The force sensor is mounted on the Kuka Kr210 robot arm, having the punch clamping system attached to it.

Results
For the experimental study, 3 rectangular specimens measuring 250 x 250 mm length were prepared and 3 parts were manufactured using the same trajectory and the same punch. The results obtained were statistically processed and have been observed to respect the normal distribution. For the median values, figures 5 ... 10 present the results for the major strain, minor strain, thinning and the three forces on the directions of the coordinate system. It should be noted that the z direction coincides with the axis of the punch. The quantitative results of the measured quantities are also presented in Table 2 in order to allow the comparative study with the simulations using the finite element method. The distribution of the major strain is observed to be relatively uniform on the conical wall, while the minor strain has maxima positioned on the generators of the conical surface. The maximum value obtained for the major strain is ε1 = 0.742 mm/mm, and for the minor strain ε2 = 0.071 mm/mm. It is observed that the thinning has a distribution similar to that of the major strain, with a maximum value of t = 56.03%. The force that develops in the direction of the punch axis (Fz) increases constantly, with local maxima, until it reaches approximately 50% of the total part height, after which it starts oscillating around this value.    Only the variation of the force in the z direction is slightly dissimilar, as it no longer presents those local maxima, but this is due to the frequency with which the results of the ASCII files were saved in the program. A lower frequency was chosen to save these results (which also include the forces) in view of the possibility to thus reduce the CPU times, already large enough. However, the analysis of figures 11 ... 22 and of the results presented in Table 2 leads to the observation that for the two cases, the results obtained are closest to the real, experimental ones. In fact, for case A5 (-16 shell formulation, local refined mesh and forming one way surface contact) the results obtained are the best in terms of the analysis accuracy.
Case A7 (-16 shell formulation, adaptive remeshing and forming one way surface contact) also has results that are very close in the case of main strain and thinning, and slightly dissimilar in the case of forces. Basically, the only difference between the two analyses is represented by the mesh. The higher percentage error in the calculation of forces is due to the fact that in case A7, which uses the adaptive remeshing as strategy, when the program detects that a remeshing is necessary, it no longer calculates all the time increments from the beginning, but only after the mesh was modified. For case A7, the percentage error for the major strain is 0.4%, for the minor strain 7%, for the forces in the x and y directions 9% and 10%, respectively, and for the force in the z direction 2%. For case A5, the percentage error for the major strain measures 3%, for the minor strain 8.4%, for the forces in the x and y directions 3% and 2.5%, respectively, and for the force in the z direction 0.4%. For the case having the biggest differences compared to the real, experimental situation, namely case A2, the percentage error for the major strain is 19%, for the minor strain 156%, for the forces in the x and y directions 33% and for the force in the z direction 39%.
The maximum time for analysis is obtained in case A4 (115613 seconds), while the minimum time is obtained in case A5 (53043 seconds).
The analysis of the sequence of figures 23 … 29, in fact of the Pareto chart for each type of measured quantity, can easily lead to the conclusion that all three analysed parameters for the main strains and thinning significantly influence the results obtained. In terms of the values of the forces, only the shell formulation and the mesh type are observed to significantly influence the results. This is obvious, as both contact types are very frequently used in the numerical simulation of metal forming processes.  The analysis of figure 29 shows that the total processing time is influenced by the shell formulation and the mesh type, as well as by the interaction between the two factors. However, it can be seen that the times required to run the analyses containing shell formulation number 25 are almost double the times required in the case of shell formulation -16.

Conclusions
The numerical simulation using the finite element method ensures the possibility of estimating the experimental results with sufficient precision, with the mention that everything must be analysed in view of the compromise between the accuracy of the analysis and the total processing time.
Following the comparison of the analytical results determined by the finite element method with the experimental results, the use of a forming one way surface contact type, a shell formulation -16 and an adaptive remeshing is recommended provided that an estimation with high precision is desired for the main strain and thinning in the SPIF process, in contrast to using a fixed mesh with local refinement if a high precision estimation of the force values in the SPIF process is desired. When using these parameters, the time allotted for the processing of the analysis is also reasonable.