A Mathematical Model for Force Prediction in Single Point Incremental Sheet Forming with Validation by Experiments and Simulation

: Incremental sheet forming (ISF) is an emerging technology that has shown great potential in forming customized three-dimensional (3D) parts without the use of product-speciﬁc dies. The forming force is reduced in ISF due to the localized nature of deformation and successive forming. Forming force plays an important role in modeling the process accurately, so it needs to be evaluated accurately. Some attempts have been made earlier to calculate the forming force; however, they are mostly limited to empirical formulae for evaluating the average forming force and its different components. The current work presents a mathematical model for force prediction during ISF in a 3D polar coordinate system. The model can be used to predict forces for axis-symmetric cones of different wall angles and also for incremental hole ﬂanging. Axial force component, resultant force in the r - θ plane, and total force have been calculated using the developed mathematical model appearing at different forming depths. The cone with the same geometrical parameters and experimental conditions was modeled and simulated on ABAQUS, and ﬁnally, experiments were carried out using a six-axis industrial robot. The mathematical model can be used to calculate forces for any wall angle, but for comparison purposes, a 45 ◦ wall angle cone has been used for analytical, numerical, and experimental validation. The total force calculated from the mathematical model had a very high level of accuracy with the force measured experimentally, and the maximum error was 4.25%. The result obtained from the FEA model also had a good level of accuracy for calculating total force, and the maximum error was 4.89%.


Introduction
Incremental Sheet Forming (ISF) is a flexible manufacturing process because complex three-dimensional (3D) parts can be fabricated by using only specific tool path programming on the existing setup without the use of dedicated dies. It is widely used for small batches and customized production. Patented in 1967, the process has attracted the eyes of aerospace, automotive [1], and biomedical sectors [2] for manufacturing complex sheet metal components. The process was used by Iseki et al. [3] for small batch production of non-symmetrical shallow shells using path-controlled spherical rollers. Various shapes were produced using the process of ISF, and a bulging test was used to predict the forming For in situ measurement of forces appearing during the ISF process, a dynamometer is fixed on the tool, and it measures the forces appearing during forming. The forces acting in ISF are (a) the axial force (Fz), acting in the direction along the tool axis, (b) force in the radial direction (Fr), and (c) force along the tangential direction (Fθ). The vector diagram of forces appearing during the process is given in Figure 1. Recently, Kumar et al. [19] presented a comprehensive review of the state of the arts of forming force and their effects on various forming outputs. Duflou et al. [20] presented a model based on regression equations for force prediction during cone forming by SPIF as a function of step size, wall angle, tool diameter, and sheet thickness. Duflou et al. [21] performed a large set of experiments to form empirical relations for the calculation of force components for five different materials: AA3033, AA5754, DC01, AISI 304, and spring steel 65Cr2. Saidi et al. [22] conducted a set of experiments on sheets of AA1050 and 304 steel and validated their results, as predicted by Jeswiet et al. [23] and Aerens et al. [21]. Xiao et al. [24] measured forming forces by an in situ force monitoring system and studied the effect of forming temperature on the forming of an AA7075-T6 sheet. They found that the peak force was reduced to 1300 N from 1900 N when forming temperature was increased from 140° to 200 °C.
In addition to various experimental measuring methodologies, numerous attempts have been made to develop analytical and FEA models, considering various deformation mechanisms based on available literature for the calculation of forces and other parameters like strain, thickness, etc., during the process.
Flores et al. [25] used the FEA simulation of SPIF and studied the effects of constitutive laws on force prediction. They concluded that the selection of constitutive laws plays a vital role in force prediction, as the results were different for elastoplastic laws with isotropic and kinematic hardening. Henrard et al. [26] used FEA simulation to predict the forces for 20° and 60° cones using the Lagamine model, first-order brick element model for various simulation parameters, and inverse method for material data fit. They also conducted experiments for the validation of results obtained from FE modeling. Honarpisheh et al. [27] performed FEM analysis on ABAQUS to predict forces appearing during forming of Al/Cu bimetals and studied the effects of various parameters in ISF.
Iseki [28] used the equilibrium of the shell element to model bulge height, strain, and load during the process. He also used FEA analysis for validation of the results of the theoretical model. Silva et al. [29] proposed membrane analysis to consider the friction between the tool and sheet for analysis of the small plastic zone created during forming of symmetrical axis shapes in SPIF. They proposed that friction stress at the contact between tool and sheet can be assumed to be made of a couple of in-plane components, viz. meridional component ( ) and circumferential component ( ). They considered the thickness, meridional, and circumferential direction as principal directions and gave the following equations for the corresponding stresses: Duflou et al. [20] presented a model based on regression equations for force prediction during cone forming by SPIF as a function of step size, wall angle, tool diameter, and sheet thickness. Duflou et al. [21] performed a large set of experiments to form empirical relations for the calculation of force components for five different materials: AA3033, AA5754, DC01, AISI 304, and spring steel 65Cr 2 . Saidi et al. [22] conducted a set of experiments on sheets of AA1050 and 304 steel and validated their results, as predicted by Jeswiet et al. [23] and Aerens et al. [21]. Xiao et al. [24] measured forming forces by an in situ force monitoring system and studied the effect of forming temperature on the forming of an AA7075-T6 sheet. They found that the peak force was reduced to 1300 N from 1900 N when forming temperature was increased from 140 to 200 • C.
In addition to various experimental measuring methodologies, numerous attempts have been made to develop analytical and FEA models, considering various deformation mechanisms based on available literature for the calculation of forces and other parameters like strain, thickness, etc., during the process.
Flores et al. [25] used the FEA simulation of SPIF and studied the effects of constitutive laws on force prediction. They concluded that the selection of constitutive laws plays a vital role in force prediction, as the results were different for elastoplastic laws with isotropic and kinematic hardening. Henrard et al. [26] used FEA simulation to predict the forces for 20 • and 60 • cones using the Lagamine model, first-order brick element model for various simulation parameters, and inverse method for material data fit. They also conducted experiments for the validation of results obtained from FE modeling. Honarpisheh et al. [27] performed FEM analysis on ABAQUS to predict forces appearing during forming of Al/Cu bimetals and studied the effects of various parameters in ISF.
Iseki [28] used the equilibrium of the shell element to model bulge height, strain, and load during the process. He also used FEA analysis for validation of the results of the theoretical model. Silva et al. [29] proposed membrane analysis to consider the friction between the tool and sheet for analysis of the small plastic zone created during forming of symmetrical axis shapes in SPIF. They proposed that friction stress at the contact between tool and sheet can be assumed to be made of a couple of in-plane components, viz. meridional component (µ Φ σ t ) and circumferential component (µ θ σ t ). They considered the thickness, meridional, and circumferential direction as principal directions and gave the following equations for the corresponding stresses: Further, Chang et al. [30] used membrane analysis for the analytical modeling of the single-pass and multi-pass incremental sheet forming process (MPISF) and incremental hole flanging process (IHF). Bansal et al. [31] used the above equations for force calculation for single-stage and multi-stage ISF, modeling the contact area and validating with experiments and the formula given by Duflou et al. [20]. Similarly, Fang et al. [32] did analytical modeling of the process and validated it with numerical modeling, and also studied the fracture behavior of the sheet in ISF. Liu et al. [33] modeled the process in the 3D polar coordinate system and validated it experimentally. They further extended their work [34] and gave a comprehensive analytical model for the prediction of forming forces in three directions and studied the effect of step depth, tool diameter, and sheet thickness on forming forces occurring in ISF.
It can be concluded from the literature survey that although attempts have been made to develop a mathematical model of ISF, a comprehensive model for force prediction has not been established yet. Further, most of the models are empirical in nature, and they are used to calculate average forces during the process. Some other models consider only plane stresses in thickness, meridional, and circumferential direction and do not consider shear stresses. Very few works have been done considering all the six stress components to finally arrive at forces. In this current work, a simple mathematical model for the ISF process has been developed in a cylindrical coordinate system considering all the stress components. All force components have been calculated at different forming depths where the tool moves in different radii. The obtained results have been compared with the FEA model developed using the ABAQUS platform. Finally, experiments have been performed on aluminum alloy 6061 on a six-axis industrial robotic arm with the same experimental conditions as used for the analytical and FEA model. Different force components were measured using a dynamometer mounted on the robotic arm. The obtained results for various force components and total force have been compared.

Mathematical Model Development
For solving the mathematical equations, the cylindrical 3D coordinate has been chosen as tool motion to give a conical shape that can be easily replicated in this coordinate system. The center of the undeformed sheet has been chosen as the origin of the global coordinate system. The element has been taken in a local 3D polar co-ordinate system which is at a forming depth h and of length rdθ in the radial direction at an angle θ in the circumferential direction on the surface of the sheet in contact with the tool in the local polar coordinate system, as shown in Figure 2a,b. It is worth noting that the three directions-radial (rdirection), circumferential (θ-direction), and axial (z-direction)-may not necessarily be the principal directions. Hence all the six-stress component has been considered.
A few assumptions listed below have been made to develop the model:

1.
All body forces have been neglected during the analysis; 2.
τ θz has been neglected throughout the analysis; 3.
Normal stress in radial and z-directions are distributed in the same ratio as that of increments provided in this direction for forming a fixed wall angle cone. ∂r ∂z = tan α, σ rr σ zz = tan α, where α is the wall angle of the formed cone. The wall angle and step depths can be understood in Figure 3a For a given radius in which the tool is moving, τ rθ varies only along θ direction. If the considered element is under static equilibrium, the equilibrium equations in the differential form can be written directly [35]. The free-body diagram of the considered element is given in Figure 3b. Neglecting body forces, the equilibrium equations in the respective directions are:  A few assumptions listed below have been made to develop the model: 1. All body forces have been neglected during the analysis; 2. has been neglected throughout the analysis; 3. Normal stress in radial and z-directions are distributed in the same ratio as that of increments provided in this direction for forming a fixed wall angle cone. = tan , = tan , where α is the wall angle of the formed cone. The wall angle and step depths can be understood in Figure 3a; 4. For a given radius in which the tool is moving, varies only along direction. If the considered element is under static equilibrium, the equilibrium equations in the differential form can be written directly [35]. The free-body diagram of the considered element is given in Figure 3b. Neglecting body forces, the equilibrium equations in the respective directions are: In radial direction In circumferential direction In z direction  The major shear component has been assumed to act only along the radial direction only. Hence has been neglected throughout the model. Further, from membrane analysis, it can be said that Using the membrane analysis (4) and Equation (1) becomes Membrane analysis (4) and the assumption 2, Equation (3) becomes For any incremental sheet forming, the fabrication of 3D shapes is done by making the tool move over the sheet and providing increments to the tool so that it can perform fabrication in different passes. For fabricating a conical wall, the increments are provided in vertical and horizontal directions. In the current work, the increments are provided in radial and the vertically downward direction. These increments are generally small in comparison to r, and if the given assumption holds, then it can be assumed that, on neglecting spring back In radial direction ∂σ rr ∂r In circumferential direction ∂τ rθ ∂r In z direction ∂τ rz ∂r The major shear component has been assumed to act only along the radial direction only. Hence τ θz has been neglected throughout the model. Further, from membrane analysis, it can be said that Processes 2023, 11, 1688 6 of 22 Using the membrane analysis (4) and Equation (1) becomes Membrane analysis (4) and the assumption 2, Equation (3) becomes For any incremental sheet forming, the fabrication of 3D shapes is done by making the tool move over the sheet and providing increments to the tool so that it can perform fabrication in different passes. For fabricating a conical wall, the increments are provided in vertical and horizontal directions. In the current work, the increments are provided in radial and the vertically downward direction. These increments are generally small in comparison to r, and if the given assumption holds, then it can be assumed that, on neglecting spring back z = r tan α, where z is the increment in the vertical direction and r is the increment in radial direction. If these are small, then ∂r tan α = ∂z (10) Using assumption 3, it can be directly said that normal stress in radial and z-directions are distributed in the same ratio as that of increments.
Taking differentials on both sides, it can be said that From (7) and (8), it can be said that Using Equation (9), Equation (6) further converges to Taking integrals on both sides of the above relation From the above analysis, it is evident that for one complete cycle, the stress in the radial direction is constant. This analysis was performed on ABAQUS with the same parameters as used for the analytical model, and considering the element when the instantaneous radius of the undeformed part is 82.5 mm, the RMS value of the S 11 was found to be approximately equal to 140 MPa as shown in Figure 3.
Since σ rr depends on instantaneous circle radius only, hence From the equilibrium equation in the r-direction, we have ∂σ rr ∂r Using Equations (18) and (20) ∂ ∂r We can directly say this using Equation (17) σ In the circumferential direction, using assumption 4, hence the equation in the circumferential direction converges to 1 r From Equations (23) and (25), we have From Von-Mises criteria, we have Processes 2023, 11, 1688 8 of 22 The model is valid for all wall angles; however, for validation purposes, a 45 • wall angle has been chosen because of its simplicity.
If we consider the strain to be acting from bending and stretching only [34], we have from volume conservation Total bending strain can be given by, ln r tool R m Whereas total tensile strain can be given as, ε stretching = ln t 0 t Total strain can be given by Taking the simple power law for annealed sheet, we have − σ = Kε n . From tensile test we have K = 205 MPa, and n = 0.2 (σ rr ) 2 From the analytical model developed, it can be said that the stress state appearing during the process is complex, but a few simple conclusions can be summarized: (a) The largest stress component is normal stress on the theta plane along the tool movement direction. Appreciable shear stress appears in the direction perpendicular to the tool motion along the radial and z planes; (b) From the Equation (28), it is clear that

Contact Area Model
For contact area estimation, the approach of Bansal et al. [31] has been implemented. The contact area has been divided into two parts. The area directly below the tool in the region below the top of the undeformed sheet and the area above the undeformed sheet, as shown in Figure 4a,b. The area above the undeformed sheet leaves contact with the sheet tangentially in the circumferential direction. Let section XX be cut, and the top sectioned view is shown in Figure 4c.
From Figure 5b, it can be directly said that Processes 2023, 11, x FOR PEER REVIEW 9 of 22 region below the top of the undeformed sheet and the area above the undeformed sheet, as shown in Figure 4a,b. The area above the undeformed sheet leaves contact with the sheet tangentially in the circumferential direction. Let section XX′ be cut, and the top sectioned view is shown in Figure 4c.   as shown in Figure 4a,b. The area above the undeformed sheet leaves contact with the sheet tangentially in the circumferential direction. Let section XX′ be cut, and the top sectioned view is shown in Figure 4c.  The area of portion AB can be directly estimated by considering the concept of solid angle. The area can be given by The area of the portion above the top of the undeformed surface can be given by integrating the perimeters of the strips of width dθ taken at an angle θ the length of which will be rβ where r = R sin θ As θ goes from ϕ to α, β goes from π 2 to 0.
Taking the variation of β to be linear with θ, we have ( From (39) and (40), we obtain The total area of contact can be given by Putting R = 5 mm, ∆z = 0.5, α = π 4 and ϕ = 5π 36 , A T = 19.05. Considering the whole contact area as a rectangle of (l × b) area. l = R(ϕ + α) = 6.10 mm and b = 3.11 mm, the length of the portion above the top of the undeformed sheet is Total projected area in r − plane = A r = (∆Z + R(cos ϕ − cos α))b sin α = 3.28 mm 2 (46) Total projected area in θ − plane = A θ = tl= 4.52 mm 2 (47) Different forces acting on the sheet in the contact region can be represented by the force diagram as given in Figure 1.
Once the relations for various force components were established, the calculation for forces using different developed equations was done on MATLAB for different instantaneous radii. The whole process was modeled on ABAQUS, and average forming forces were evaluated from the FEA model for the same instantaneous radii. The values obtained from the analytical model were compared with FEA and experimental results, as discussed in subsequent sections.

Simulation Study
The simulation analysis on the formed cone was carried out on ABAQUS. The finite element model was made on the software, which consisted of a rigid tool with a tool diameter of 10 mm, a deformable aluminum alloy 6061 sheet of thickness 1.05 mm, and a flange to hold the sheet. The geometry of the model used for simulation is given in Figure 5.
The properties of the material used in the simulation were obtained from the uniaxial tensile test and are given in Table 1. For simulation, a true stress-true strain curve up to maximum load was obtained from the tensile test. Since a large strain is obtained in ISF for post-necking values of true stress and true strain, the following method was adopted. The log-log plot of true stress and the true strain was made till necking. Thereafter, a linear relation was assumed between the logarithm of true stress and the logarithm of true strain from necking to fracture. The strain hardening exponent in this region can be evaluated by finding the slope of this linear portion of the plot from necking to fracture. In this way, the true stress-true strain data was obtained and fed in ABAQUS till fracture. The engineering and true stress-strain curve obtained after tensile testing is given in Figure 6.

Simulation Study
The simulation analysis on the formed cone was carried out on ABAQUS. The finite element model was made on the software, which consisted of a rigid tool with a tool diameter of 10 mm, a deformable aluminum alloy 6061 sheet of thickness 1.05 mm, and a flange to hold the sheet. The geometry of the model used for simulation is given in Figure  5.
The properties of the material used in the simulation were obtained from the uniaxial tensile test and are given in Table 1. For simulation, a true stress-true strain curve up to maximum load was obtained from the tensile test. Since a large strain is obtained in ISF for post-necking values of true stress and true strain, the following method was adopted. The log-log plot of true stress and the true strain was made till necking. Thereafter, a linear relation was assumed between the logarithm of true stress and the logarithm of true strain from necking to fracture. The strain hardening exponent in this region can be evaluated by finding the slope of this linear portion of the plot from necking to fracture. In this way, the true stress-true strain data was obtained and fed in ABAQUS till fracture. The engineering and true stressstrain curve obtained after tensile testing is given in Figure 6.  The isotropic hardening model was used for analysis as the sheet was annealed. Square meshing was done on the plate and tool with an average mesh size of 6.96 mm.
Loading and Boundary condition: Penalty contact between the tool and sheet was taken with the coefficient of friction between the tool and the sheet to be 0.2, which was obtained from the wear test. The boundary of the sheet was encastred, and the blank was made rigid and encastred. The tip of the tool was chosen as a reference point and was given corresponding motions as per the conical geometry. The tool path was generated on MATLAB, and the corresponding positions were fed in ABAQUS to obtain the conical shape. The simulation was run, and dynamic analysis was done for stresses, strains, and tool forces arising during the tool motion, the results of which are given in Figure 7.
The isotropic hardening model was used for analysis as the sheet was annealed. Square meshing was done on the plate and tool with an average mesh size of 6.96 mm.
Loading and Boundary condition: Penalty contact between the tool and sheet was taken with the coefficient of friction between the tool and the sheet to be 0.2, which was obtained from the wear test. The boundary of the sheet was encastred, and the blank was made rigid and encastred. The tip of the tool was chosen as a reference point and was given corresponding motions as per the conical geometry. The tool path was generated on MATLAB, and the corresponding positions were fed in ABAQUS to obtain the conical shape. The simulation was run, and dynamic analysis was done for stresses, strains, and tool forces arising during the tool motion, the results of which are given in Figure 7. As can be seen from Figure 7a, the Von-Mises stress is maximum in the direct contact region between the undeformed sheet and the tool. The middle region of the cone is strained to the maximum value, as can be seen in Figure 7b. Further, it can also be seen from Figure 7c that the middle region of the cone has the minimum thickness showing that this region has undergone the maximum amount of thinning, which also affirms the distribution as seen in Figure 7b. A path along the central meridional plane was chosen in the cone after deformation and strain occurred along the path, and the thickness of the sheet along the path was plotted, which is shown in Figure 8a,b. The simulation analysis was used to find out the forces appearing on the tool during the deformation. The tool force has been assumed to have two components. The axial As can be seen from Figure 7a, the Von-Mises stress is maximum in the direct contact region between the undeformed sheet and the tool. The middle region of the cone is strained to the maximum value, as can be seen in Figure 7b. Further, it can also be seen from Figure 7c that the middle region of the cone has the minimum thickness showing that this region has undergone the maximum amount of thinning, which also affirms the distribution as seen in Figure 7b. A path along the central meridional plane was chosen in the cone after deformation and strain occurred along the path, and the thickness of the sheet along the path was plotted, which is shown in Figure 8a The isotropic hardening model was used for analysis as the sheet was annealed. Square meshing was done on the plate and tool with an average mesh size of 6.96 mm.
Loading and Boundary condition: Penalty contact between the tool and sheet was taken with the coefficient of friction between the tool and the sheet to be 0.2, which was obtained from the wear test. The boundary of the sheet was encastred, and the blank was made rigid and encastred. The tip of the tool was chosen as a reference point and was given corresponding motions as per the conical geometry. The tool path was generated on MATLAB, and the corresponding positions were fed in ABAQUS to obtain the conical shape. The simulation was run, and dynamic analysis was done for stresses, strains, and tool forces arising during the tool motion, the results of which are given in Figure 7. As can be seen from Figure 7a, the Von-Mises stress is maximum in the direct contact region between the undeformed sheet and the tool. The middle region of the cone is strained to the maximum value, as can be seen in Figure 7b. Further, it can also be seen from Figure 7c that the middle region of the cone has the minimum thickness showing that this region has undergone the maximum amount of thinning, which also affirms the distribution as seen in Figure 7b. A path along the central meridional plane was chosen in the cone after deformation and strain occurred along the path, and the thickness of the sheet along the path was plotted, which is shown in Figure 8a,b. The simulation analysis was used to find out the forces appearing on the tool during the deformation. The tool force has been assumed to have two components. The axial The simulation analysis was used to find out the forces appearing on the tool during the deformation. The tool force has been assumed to have two components. The axial components of the force F zz and the resultant of in-plane X and Y components of forces is also equal to the resultant of tangential and radial force components.
As can be seen from Figure 9, the largest force component arising during the deformation is in the axial direction, as was reported by Duflou et al. [20] as well. Additionally, it was observed that the largest axial force appearing in the axial direction was found to be 2629.11 N. The mean force was found to be 1512.24 N. The variation of in-plane force, that is, the resultant of the forces appearing in the plane of the undeformed sheet with time, has been shown in Figure 9b, and the force diagram for the same has been given in Figure 9c. The largest value corresponding to in-plane forming force was found to be 2172.31 N, and the average force during forming was found to be 1109.3 N. Finally, the variation of total forming force with time has been given in Figure 9d. The maximum forming force appearing on the tool was found to be 2747.46 N, and the average tool force was found to be 1812 N.
components of the force Fzz and the resultant of in-plane X and Y components of forces is also equal to the resultant of tangential and radial force components.
As can be seen from Figure 9, the largest force component arising during the deformation is in the axial direction, as was reported by Duflou et al. [20] as well. Additionally, it was observed that the largest axial force appearing in the axial direction was found to be 2629.11 N. The mean force was found to be 1512.24 N. The variation of in-plane force, that is, the resultant of the forces appearing in the plane of the undeformed sheet with time, has been shown in Figure 9b, and the force diagram for the same has been given in Figure 9c. The largest value corresponding to in-plane forming force was found to be 2172.31 N, and the average force during forming was found to be 1109.3 N. Finally, the variation of total forming force with time has been given in Figure 9d. The maximum forming force appearing on the tool was found to be 2747.46 N, and the average tool force was found to be 1812 N.

Experimental Validation
Once the analytical model was developed and force calculations were done, a simulation was carried out for preliminary validation of obtained forces analytically. Finally, experiments were carried out on a six-axis industrial robotic arm. 1.05 mm thick aluminum 6061 alloy sheets were taken for experimentation. The material composition is given in Table 2.

Experimental Validation
Once the analytical model was developed and force calculations were done, a simulation was carried out for preliminary validation of obtained forces analytically. Finally, experiments were carried out on a six-axis industrial robotic arm. 1.05 mm thick aluminum 6061 alloy sheets were taken for experimentation. The material composition is given in Table 2. Tensile tests and Erichsen cup tests were carried out for the determination of material properties before conducting the experiments. The sample was considered to be isotropic. The tensile test was carried out on the 100 kN INSTRON 8801 model. ASTM/E8/E8M standards were followed for the making of the specimens [36]. The sample after the test is shown in Figure 10.
Tensile tests and Erichsen cup tests were carried out for the determination of material properties before conducting the experiments. The sample was considered to be isotropic. The tensile test was carried out on the 100 kN INSTRON 8801 model. ASTM/E8/E8M standards were followed for the making of the specimens [36]. The sample after the test is shown in Figure 10. Once the test was conducted, the results were calculated using the data obtained from the test, which are summarized in Table 3. Before selecting a material, it is very crucial to have an idea of its formability. For the formability prediction of the sheet to be formed, the Erichsen cup test setup was used. The diameter of the used indenter was 20 mm. The main scale division was 1 mm, and the circular scale division was 50/5 MSD. The measured depth (HE) for all the domes was measured, which is given in Table 4. The domes are made as a result of indentation made by the indenter, as shown in Figure 11.  Once the test was conducted, the results were calculated using the data obtained from the test, which are summarized in Table 3. Before selecting a material, it is very crucial to have an idea of its formability. For the formability prediction of the sheet to be formed, the Erichsen cup test setup was used. The diameter of the used indenter was 20 mm. The main scale division was 1 mm, and the circular scale division was 50/5 MSD. The measured depth (H E ) for all the domes was measured, which is given in Table 4. The domes are made as a result of indentation made by the indenter, as shown in Figure 11. After obtaining the mechanical parameters from the uniaxial tensile test and Erichsen cup test, experiments were performed for which the chosen parameters are given in Table 5. The parameters were obtained using surface regression models. Dimension of the undeformed workpiece was taken as 250 mm × 250 mm × 1.05 mm. The sheet was deformed into a cone with a radius of the top circle 96 mm. The whole setup, the forming tool with the tool dynamometer, and the deformed axis-symmetric cone are shown in Figure 12a-c.
Before selecting a material, it is very crucial to have an idea of its formability. For the formability prediction of the sheet to be formed, the Erichsen cup test setup was used. The diameter of the used indenter was 20 mm. The main scale division was 1 mm, and the circular scale division was 50/5 MSD. The measured depth (HE) for all the domes was measured, which is given in Table 4. The domes are made as a result of indentation made by the indenter, as shown in Figure 11. Figure 11. Domes made after the Erichsen cup test. Figure 11. Domes made after the Erichsen cup test.  After obtaining the mechanical parameters from the uniaxial tensile test and Erichsen cup test, experiments were performed for which the chosen parameters are given in Table  5. The parameters were obtained using surface regression models. Dimension of the undeformed workpiece was taken as 250 mm × 250 mm × 1.05 mm. The sheet was deformed into a cone with a radius of the top circle 96 mm. The whole setup, the forming tool with the tool dynamometer, and the deformed axis-symmetric cone are shown in Figure 12ac.

Thickness Distribution
The instantaneous thickness has a role to play in Equation (18). Hence the state of stress and the force are affected by the instantaneous sheet thickness, so the value of the sheet thickness was obtained from the FEA model and was validated by measuring the values after the experiments.  After obtaining the mechanical parameters from the uniaxial tensile test and Erichsen cup test, experiments were performed for which the chosen parameters are given in Table  5. The parameters were obtained using surface regression models. Dimension of the undeformed workpiece was taken as 250 mm × 250 mm × 1.05 mm. The sheet was deformed into a cone with a radius of the top circle 96 mm. The whole setup, the forming tool with the tool dynamometer, and the deformed axis-symmetric cone are shown in Figure 12ac.

Thickness Distribution
The instantaneous thickness has a role to play in Equation (18). Hence the state of stress and the force are affected by the instantaneous sheet thickness, so the value of the sheet thickness was obtained from the FEA model and was validated by measuring the values after the experiments.
Once the experiments had been carried out, the thickness of the sheet was measured. The measurement was done by a micrometer with conical tips of least count 0.001 mm. The sheet was found to have the minimum thickness in the middle of the cone. The trend was also shown in previous work by Ambrogio et al. [37], who showed that the thickness prediction by the sine law is most accurate in the middle region of the deformed sheet.
The conical region was divided into seven sub-regions (0-6), as shown in Figure 13a, and the thickness in the different regions was measured. Region-0 is the un-deformed region where thickness is 1.05 mm, and Region-6 is the region where the material gets accumulated because the tool moving along the sheet drags the material with it, and finally, it gets accumulated in the lowest region of the sheet; hence, the sheet thickness exceeds the original sheet thickness. However, Region-1-5 is the region in which sheet thinning occurs. The thickness plot for ISF in FEA and experiments are given in Figure 13c. For thickness prediction using FEA analysis, a path was chosen in the meridional direction, as shown in Figure 13b. The thickness variation obtained from FEA is shown in Figure 13c. All three thickness variations, obtained experimentally by FEA and by sine law [4], are plotted along the meridional path, as shown in Figure 13. The sheet was found to have the minimum thickness in the middle of the cone. The trend was also shown in previous work by Ambrogio et al. [37], who showed that the thickness prediction by the sine law is most accurate in the middle region of the deformed sheet. The conical region was divided into seven sub-regions (0-6), as shown in Figure 13a, and the thickness in the different regions was measured. Region-0 is the un-deformed region where thickness is 1.05 mm, and Region-6 is the region where the material gets accumulated because the tool moving along the sheet drags the material with it, and finally, it gets accumulated in the lowest region of the sheet; hence, the sheet thickness exceeds the original sheet thickness. However, Region-1-5 is the region in which sheet thinning occurs. The thickness plot for ISF in FEA and experiments are given in Figure 13c. For thickness prediction using FEA analysis, a path was chosen in the meridional direction, as shown in Figure 13b. The thickness variation obtained from FEA is shown in Figure 13c. All three thickness variations, obtained experimentally by FEA and by sine law [4], are plotted along the meridional path, as shown in Figure 13.

Force Measurement
For the measurement of forces, the cone was divided into seven zones. Zone 0 was taken as the zone of the upper undeformed cone, the middle region was taken as Zones 1-5 in which near uniform thinning occurs, and the sheet thinning followed the sine law, and Zone 6 was the region between the lower undeformed sheet and the cone, as shown in Figure 13a. The force was measured during deformation using a drill tool dynamometer. The setup is shown in Figure 12a,b. The capacity of the dynamometer is 500 kgf having a strain gauge based 350 Ω bridge sensor. The force was measured in X, Y, and Z directions, and the resultant of the forces in X and Y directions was calculated. Three forces were defined as Fz, Fp (resultant of forces in X and Y directions), and total force (Ft). Average values of Fz, Fp, and Ft were calculated for one complete cycle and were compared with the values obtained analytically. For comparison with the forces obtained from ABAQUS analysis, the average of forces for one complete circular cycle was calculated. Thus, all the components of forces were obtained from the analytical model, as well as the FEA model, using ABAQUS, and final validation was done by experimental results. The corresponding values of forces were plotted against the values of the instantaneous radius of the cone.

Variation of Axial Force
This force acts in the vertical direction. The force was calculated via an analytical model as well as the FEA model. Both were validated using final experimentation. The

Force Measurement
For the measurement of forces, the cone was divided into seven zones. Zone 0 was taken as the zone of the upper undeformed cone, the middle region was taken as Zones 1-5 in which near uniform thinning occurs, and the sheet thinning followed the sine law, and Zone 6 was the region between the lower undeformed sheet and the cone, as shown in Figure 13a. The force was measured during deformation using a drill tool dynamometer. The setup is shown in Figure 12a,b. The capacity of the dynamometer is 500 kgf having a strain gauge based 350 Ω bridge sensor. The force was measured in X, Y, and Z directions, and the resultant of the forces in X and Y directions was calculated. Three forces were defined as F z , F p (resultant of forces in X and Y directions), and total force (F t ). Average values of F z , F p , and F t were calculated for one complete cycle and were compared with the values obtained analytically. For comparison with the forces obtained from ABAQUS analysis, the average of forces for one complete circular cycle was calculated. Thus, all the components of forces were obtained from the analytical model, as well as the FEA model, using ABAQUS, and final validation was done by experimental results. The corresponding values of forces were plotted against the values of the instantaneous radius of the cone.

Variation of Axial Force
This force acts in the vertical direction. The force was calculated via an analytical model as well as the FEA model. Both were validated using final experimentation. The peak and the average value of the axial forces obtained in all three cases are given in Table 6. The axial force increased with decreasing radius, as observed in the analytical model. The force in the axial direction is given by F zz = σ zz A z + z rz A r . The values of forces were calculated for different radii of the undeformed sheet using the formula, and the force was found to increase with the decreasing radius. The variation of the axial force with decreasing radii was almost linear. The same trend was seen in the FEA model, and experimentally, it was seen to increase with decreasing radii, but the shape of the plot departed from linearity slightly in the middle zone, as shown in Figure 14a,b. This can be due to sheet thinning, which can give rise to inaccurate sheet thickness at different points for a given radius. Nevertheless, the model gives a very good approximation of axial forming force during deformation. The error function for the error of the forces obtained analytically was defined as peak and the average value of the axial forces obtained in all three cases are given in Table  6. The axial force increased with decreasing radius, as observed in the analytical model. The force in the axial direction is given by = + . The values of forces were calculated for different radii of the undeformed sheet using the formula, and the force was found to increase with the decreasing radius. The variation of the axial force with decreasing radii was almost linear. The same trend was seen in the FEA model, and experimentally, it was seen to increase with decreasing radii, but the shape of the plot departed from linearity slightly in the middle zone, as shown in Figure 14a,b. This can be due to sheet thinning, which can give rise to inaccurate sheet thickness at different points for a given radius. Nevertheless, the model gives a very good approximation of axial forming force during deformation. The error function for the error of the forces obtained analytically was defined as Percentage error = × 100 (51) The maximum error was found to be 4.92%. The force trends, as observed in ABAQUS, also closely matched the data obtained experimentally. The maximum error obtained in this case was 5.52%. The error in the axial force with respect to the experimental force calculated at every value of radius is shown in Section 6.4.

Variation of Resultant Force in r-θ Plane
The force in the r-θ plane has two components. For the mathematical model, the two components of forces evaluated were the ones in the r and θ planes. However, the components of the forces measured by the dynamometer and the forces observed in the ABAQUS analysis are in X and Y directions. So, the resultant of Fr and Fθ was calculated, The maximum error was found to be 4.92%. The force trends, as observed in ABAQUS, also closely matched the data obtained experimentally. The maximum error obtained in this case was 5.52%. The error in the axial force with respect to the experimental force calculated at every value of radius is shown in Section 6.4.

Variation of Resultant Force in r-θ Plane
The force in the r-θ plane has two components. For the mathematical model, the two components of forces evaluated were the ones in the r and θ planes. However, the components of the forces measured by the dynamometer and the forces observed in the ABAQUS analysis are in X and Y directions. So, the resultant of F r and F θ was calculated, which was equal to the resultant of F X and F Y . Finally, the force obtained analytically, and that obtained in FEA was compared with the experimental force. The peak and the average value of the resultant force in r-θ obtained in all three cases are given in Table 7. The plot of forces obtained in the r-θ plane from the analytical model, FEA, and experiment is given in Figure 15. The force obtained analytically had a good level of accuracy with experimental force. The largest error was found to be 8.47%, and the force obtained on FEA was found to have a maximum deviation of 4.61% from the values measured experimentally. The error in the axial force with respect to the experimental force calculated at every value of radius is shown in Section 6.4.
which was equal to the resultant of FX and FY. Finally, the force obtained analytically, and that obtained in FEA was compared with the experimental force. The peak and the average value of the resultant force in r-θ obtained in all three cases are given in Table 7. The plot of forces obtained in the r-θ plane from the analytical model, FEA, and experiment is given in Figure 15. The force obtained analytically had a good level of accuracy with experimental force. The largest error was found to be 8.47%, and the force obtained on FEA was found to have a maximum deviation of 4.61% from the values measured experimentally. The error in the axial force with respect to the experimental force calculated at every value of radius is shown in Section 6.4.

Variation of Total Force
In the analytical model, the total force appearing on the forming tool-sheet interface was calculated by the formula = + . Once the force was calculated, it was compared with the total force obtained experimentally and the one obtained in the FEA model. The peak and the average value of the total force obtained in all three cases are given in Table 8.

Variation of Total Force
In the analytical model, the total force appearing on the forming tool-sheet interface was calculated by the formula F t = F p 2 + F z 2 . Once the force was calculated, it was compared with the total force obtained experimentally and the one obtained in the FEA model. The peak and the average value of the total force obtained in all three cases are given in Table 8. The plot of total force obtained analytically, by the FEA model, and experimentally is given in Figure 16. The computed axial force from the mathematical model has a good level of accuracy as far as the total force is concerned. The maximum error in the calculated force from the total experimental force was found to be 4.25%. The plot of total force obtained analytically, by the FEA model, and experimentally is given in Figure 16. The computed axial force from the mathematical model has a good level of accuracy as far as the total force is concerned. The maximum error in the calculated force from the total experimental force was found to be 4.25%. The total force was also obtained from the FEA model on ABAQUS. The force obtained from the FEA model was also in close agreement with the total force obtained experimentally. The maximum error in obtained total force in the FEA model from the experimental model was found to be 4.89%, as shown in Figure 17c.

Error in Calculated and Experimentally Determined Forces
The forces calculated from the analytical model were compared with forces obtained by the FEA model and experimentally. For validation purposes, the true value of force was taken as the one which was obtained experimentally. The values of the forces were obtained from the analytical model for different instantaneous radii, and the corresponding values of the forces were also obtained from ABAQUS data. Both the values were compared with experimental values. So, validation of the mathematical model was done by comparing the values of the forces with the experimental values.
The error function was defined as percentage error = × 100 For axial forces and the total forces, the error was maximized in Zone 2, i.e., the middle region of the cone. However, for the resultant force in the r-θ plane, the error was maximum in Region-1, which decreased, and then the error was increased in Region-2 before finally decreasing in Region-3. The total force was also obtained from the FEA model on ABAQUS. The force obtained from the FEA model was also in close agreement with the total force obtained experimentally. The maximum error in obtained total force in the FEA model from the experimental model was found to be 4.89%, as shown in Figure 17c.
The plot of total force obtained analytically, by the FEA model, and experimentally is given in Figure 16. The computed axial force from the mathematical model has a good level of accuracy as far as the total force is concerned. The maximum error in the calculated force from the total experimental force was found to be 4.25%. The total force was also obtained from the FEA model on ABAQUS. The force obtained from the FEA model was also in close agreement with the total force obtained experimentally. The maximum error in obtained total force in the FEA model from the experimental model was found to be 4.89%, as shown in Figure 17c.

Error in Calculated and Experimentally Determined Forces
The forces calculated from the analytical model were compared with forces obtained by the FEA model and experimentally. For validation purposes, the true value of force was taken as the one which was obtained experimentally. The values of the forces were obtained from the analytical model for different instantaneous radii, and the corresponding values of the forces were also obtained from ABAQUS data. Both the values were compared with experimental values. So, validation of the mathematical model was done by comparing the values of the forces with the experimental values.
The error function was defined as percentage error = × 100 For axial forces and the total forces, the error was maximized in Zone 2, i.e., the middle region of the cone. However, for the resultant force in the r-θ plane, the error was maximum in Region-1, which decreased, and then the error was increased in Region-2 before finally decreasing in Region-3.

Error in Calculated and Experimentally Determined Forces
The forces calculated from the analytical model were compared with forces obtained by the FEA model and experimentally. For validation purposes, the true value of force was taken as the one which was obtained experimentally. The values of the forces were obtained from the analytical model for different instantaneous radii, and the corresponding values of the forces were also obtained from ABAQUS data. Both the values were compared with experimental values. So, validation of the mathematical model was done by comparing the values of the forces with the experimental values.
The error function was defined as percentage error = F calculated − F experimental F experimental × 100 For axial forces and the total forces, the error was maximized in Zone 2, i.e., the middle region of the cone. However, for the resultant force in the r-θ plane, the error was maximum in Region-1, which decreased, and then the error was increased in Region-2 before finally decreasing in Region-3.
As far as the FEA model is concerned, the validation was done by comparing the forces obtained experimentally. The error for all three forces followed the same trend. The percentage error vs. instantaneous radius followed a bell curve pattern. It was maximized in Region-2 and was minimal in Region-1 and Region-3. The error curves are given in Figure 17a-c.

Conclusions
This paper presents an analytical model to predict the forming forces during incremental sheet forming. The modeling was done in a polar 3D coordinate system. All six stress components were taken into consideration initially before neglecting τ θz . Analytical formulae were developed for all the stress components. The contact area was modeled analytically before finally developing the equations for all forces. The three force components for which the equations were developed were F r , F θ , and F z . The result of the three components is the total force F t . The resultant of F r and F θ is the resultant force in the r-θ plane, represented as F p . The material used in the mathematical model was AA6061. The process was modeled in FEA using the same experimental parameters, and the force components F z and F p were calculated before finally calculating the total force (F t ). Finally, the process was carried out on a six-axis industrial robot, and forces were measured using a tool dynamometer. The data obtained from the analytical model and FEA were validated by the experimental results. The model is for all wall angle cones; however, in this work, the validation has been done for a 45 • cone. The results obtained from the analytical model and FEA were in good agreement with the experimental result. The maximum error in the analytically calculated total force was found to be 4.25%. Similarly, the results obtained from the FEA model were also in agreement with the experimental results. The maximum error in evaluated total force from the FEA model was found to be 4.89% with respect to the experimental results. It can be concluded that the mathematical model developed can be used for a force calculation to a good level of accuracy.  Abbreviations σ zz Normal stress in axial direction σ rr Normal stress in radial direction σ θθ Normal stress in circumferential direction τ rθ Shear stress in radial plane in circumferential direction = τ θr τ rz Shear stress in radial plane in axial direction = τ zr τ zθ Shear stress in axial plane in circumferential direction = τ θz f Coefficient of kinetic friction t o Initial thickness of the sheet t Instantaneous thickness of the sheet r Radius of current pass (Instantaneous radius) R t Tool radius R m Radius of curvature of neutral axis of the taken element