A Novel Algorithm for Thickness Prediction in Incremental Sheet Metal Forming

Incremental sheet metal forming characterized as increased flexibility and local plastic deformation is well suitable for low-production-run manufacturing and a new sample trial production of complex shapes. Thickness thinning is still an obstacle to the application of incremental forming. In this study, a novel mathematical algorithm based on a non-uniform rational B-spline (NURBS) surface was proposed and implemented which focuses on predicting and calculating the final thickness for arbitrary parts in incremental forming. In order to evaluate the validity of the proposed model, the finite element simulation and forming experiments of three kinds of parts, such as truncated cones, truncated pyramids and ellipsoid parts, were conducted. The thickness of theoretical prediction was compared with that of finite element simulation and experiment, and good agreements were obtained. The results show that the proposed model and the method are effective and robust for predicting the thickness of the formed parts in incremental sheet metal forming.


Introduction
Incremental sheet metal forming (ISMF) is an advanced rapid die-less flexible sheet forming technology for prototyping and low-production-run manufacturing of almost arbitrary 3D shapes. A variety of incremental forming (IF) technologies have been evolved, such as single-point incremental forming (SPIF), two-point incremental forming (TPIF) [1], double-sided incremental forming (DSIF) [2] and three opposite point incremental forming (TOPIF) [3]. The remarkable feature of the ISMF is that the sheet will be continuously thinned in the forming area [4]. During the ISMF process, thickness thinning will be induced as the localized plastic deformation is accumulated to produce the final part. So, thickness distribution is a critical factor for designing and controlling the incremental forming operations [5,6]. It is known that the sine law, verified on truncated cones, truncated pyramids and other simple formed parts is a concise model for prediction of thickness in the IF process [7][8][9][10][11][12]. Bambach presented a variant of the sine law by assuming that the sheet deforms along the normal direction of the smooth transient surface [13]. Malhotra et al. introduced a squeezing coefficient into the sine law for thickness prediction in DSIF [14]. Salem et al. presented an algorithm to predict and calculate the actual thickness by using the plane normal vector of the distorted element across the sheet [15]. Choi and Lee presented a mathematical prediction model of thickness distribution of the hybrid forming of ISMF and stretch forming based on the conventional sine law [16]. However, these methods have a great disadvantage in that they have a large number of calculations for the thickness prediction of parts with varied forming angles. Furthermore, they are not available for the thickness prediction of complicated parts with random shapes which are not represented directly by the mathematical function, such as F(x, y, z) = 0.
The present study focused on developing a general model based on NURBS surface, which can be used to predict and calculate the final thickness for arbitrary parts.
The freeform surface of an arbitrary part was reconstructed based on a non-uniform rational B-spline (NURBS) parametric surface. The first-order partial derivatives of the freeform surface were evaluated to solve the area. Based on the obtained area element, the thickness was calculated according to the principle of volume constancy in plastic deformation. In fact, the proposed model is an extension of the sine law. Taking the truncated cones with different forming angles, pyramid part and hemi-ellipsoid part for examples, experiments are performed to show that the proposed model can predict the final thickness for complicated parts with a random shape.

Mathematical Model and Methods for Thickness Calculation
In this work, the freeform surface of the forming part is described by using a NURBS surface which is a parametric polynomial surface.
A NURBS surface with the degrees of p and q along the parametric direction of u and v, respectively, can be defined as follows: where P(x, y, z) is the arbitrary point on the NURBS surface. i = 0, · · · , m, j = 0, · · · , n, P i,j are the control points, ω i,j are the weight coefficients, (m + 1) and (n + 1) are the amount of control points, and N i,p (u) and N j,q (v) are the B-spline basis functions of order p and q, respectively, defined over the knot vectors U = , v q+1 , · · · , v n , 1, · · · , 1 q+1   , which are calculated by: According to the first fundamental form of a surface, the surface area can be expressed as: where D is the parameter domain, and f u and f v are the partial derivatives of the surface f with respect to the parametric variable u and v, respectively. f u × f v is just the cross product, which can be calculated as follows: In differential notation, the surface area of the patch of the surface is: Then, we can obtain Due to volume constancy in plastic deformation, we obtain where t 0 and t represent the thickness of the original blank sheet and the actual formed sheet, respectively. Substitution of Equation (4) into Equation (5) gives: The algorithm presented in this study is rather general in predicting thickness. It holds for the thickness prediction of simple parts and complicated parts with a random shape. In fact, it can also be found that the proposed algorithm is the same as the sine law for the thickness prediction of parts with constant forming angle, such as conical shape and pyramid shape. In other words, its advantage is to predict and calculate the thickness of parts with a variable forming angle. In a few words, this algorithm is more convenient than the sine law to calculate the thickness once the surface equation of the desired component is obtained.
The procedure of the proposed mathematical modeling approach can be summarized as follows: 1.
Get the Cartesian or coordinates of the discrete points; 2.
Compute the first derivatives of the NURBS surface; 5.
Compute the cross product of two partial derivatives; 6.
Compute the final thickness.

Experimental Results and Discussion
Based on the above results, to verify the validity of the relation developed in Equation (6), the truncated cone part, the truncated pyramid part and the ellipsoid part have been designed and formed by a custom-designed and fabricated machine shown in Figure 1. Figure 1a shows the schematic of the mechanical structure. During the forming process, the blank holder driven by the servo motor and the sheet clamped on the blank holder rotate together. The forming tools driven by the stepper motors can move along the horizontal and vertical direction. The prototype of the deformation machine shown as Figure 1b was manufactured by Zhejiang University of Technology. This machine can be used to form the hyperbolic parts online by selecting and controlling different forming tools, benefiting from the design of its bilateral symmetrical structure. The maximum displacements of the forming head are 200 mm and 100 mm along the vertical and horizontal direction, respectively. The diameter of the annular blank holder is 160 mm, which determines the size of the blank sheet in the experiment. So, the diameter of the circular aluminum sheet used in experiments is determined to be 160 mm. The initial thickness of the aluminum sheet is 1 mm. The diameter of the forming tool with hemispherical head which is made of tool steel is 6 mm. In order to enable the final thickness to be as uniform as possible and obtain better surface quality, the oil lubrication was used during forming process in the following experiments. The material parameters of aluminum sheet and forming tool are shown in Table 1.   In order to better grasp the thickness distribution, the finite element simulation was performed using ANSYS/LS-DYNA. As the diameter of the aluminum sheet is 160 times than its thickness, the sheet can be regarded as a thin sheet in the simulation. The element type of the SHELL 163 was adopted for the aluminum sheet. Considering the influence of time and accuracy, the Belytschko-Wong-Chiang was used as the element formulation.  In order to better grasp the thickness distribution, the finite element simulation was performed using ANSYS/LS-DYNA. As the diameter of the aluminum sheet is 160 times than its thickness, the sheet can be regarded as a thin sheet in the simulation. The element type of the SHELL 163 was adopted for the aluminum sheet. Considering the influence of time and accuracy, the Belytschko-Wong-Chiang was used as the element formulation. Accordingly, the element type of SOLID 164 was adopted for the forming tool modeled as a rigid body. The Barlat model with the exponential anisotropic hardening assumption is defined as the constitutive model of aluminum sheet. The yield rule can be expressed as: where σ s = δε η is the yield stress, δ and η are the constants determined by experiment, ε is the strain. k is the Barrat constant, , α, β, λ and µ are the constants of the anisotropic material. σ 1 and σ 2 are the principal stress, τ 12 is the shear stress. Considering the use of oil lubrication in the experiment, the surface to surface contact in which the boundary friction with friction coefficient of 0.1 was introduced was defined between the aluminum sheet and the forming head. The outer edge of the aluminum sheet is constrained. The mesh size of the sheet metal deformation area and the forming tool are 1 mm and 2 mm, respectively.
In the experiments and simulations, some parameters including the coordinates of the starting point, radial feed rate, axial step down and rotational speed of the blank holder with sheet need to be set and input, according to the contour path. The forming parameters of the formed parts in this study are tabulated in Table 2. It should be noted that the coordination between sheet rotation and the feed of forming head needs to be planned and designed for pyramid and ellipsoid parts.

Case 1: Truncated Cone Shape
A cone as a 2 × 2 degree NURBS surface with 9 × 3 control points shown in Figure 2 can be constructed as: Two kinds of truncated cones with forming angles of 76 • and 53. can be constructed as: where the knot vectors are × = 0,0,0,0.25,0.25,0.5,0.5,0.75,0.75,1,1,1 and × = 0,0,0,1,1,1 . The array of weights is , =   The two truncated cones with the different forming angles were formed by using the aluminum blank sheet with thickness of t 0 = 1 mm. The step down of ∆z = 0.5 mm and ∆z = 0.75 mm shown in Table 2 were used in the tests, respectively. During the forming process, the oil lubricant was painted on the sheet to improve the surface quality by reducing the friction caused by the forming tool. The formed cones are shown in Figure 4.
In order to evaluate the proposed theoretical model, the comparison of theoretical, simulated and experimental thickness was performed. A micrometer with the accuracy of 0.01 mm was used to measure the thickness of the experimental sheet. In order to reduce the errors as much as possible, the formed parts are cut using wire cutting to obtain the measured sample, and then the average value of five measurements is taken as the actual thickness. The simulated thickness is extracted from the thickness fields based on LS-PrePost, which is solved using the finite element method (FEM). The theoretical thickness is calculated using the presented algorithm according to Equation (6). The corresponding thicknesses for the truncated cones are plotted in Figure 5 as a comparison. Figure 5a shows the thickness profile of the truncated cone with a forming angle of 76 • away from the central axis at x = 0 mm along the radial direction. The bottom surface of the truncated cone located within the range of x = 0-20 mm, which belongs to the unformed area, gradually moves downward as a whole, when the forming is carried out. The inclined wall located within the range of x = 20-40 mm is the formed area. The area over x = 40 mm, of which the outer edge is clamped by blank holder, belongs to the unformed area. The thickness of the formed region is basically unchanged because of the constant forming angle. The average thickness of experiment and simulation are 0.950 mm and 0.955 mm, respectively. The result predicted by this presented algorithm is 0.97 mm for the truncated cone with the forming angle of 76 • , which is the same as the result of sin(76 • ) = 0.97 mm calculated using the traditional sine law. The thickness thinning was found in the transition regions between the unformed and formed areas. The center and outer edge of the sheet basically maintain the original thickness. The same situation was found in Figure 5b. The average thicknesses of the sheet in the formed area located within the range of x = 10-30 mm are 0.788 mm and 0.797 mm, respectively, in the experiment and simulation. The predicted result is 0.80 mm for the truncated cone with the forming angle of 53.2 • , which is equivalent to the result of sin(53.2 • ) = 0.8 mm calculated using the sine law. The two truncated cones with the different forming angles were formed by using the aluminum blank sheet with thickness of t0 = 1 mm. The step down of Δz = 0.5 mm and Δz = 0.75 mm shown in Table 2 were used in the tests, respectively. During the forming process, the oil lubricant was painted on the sheet to improve the surface quality by reducing the friction caused by the forming tool. The formed cones are shown in Figure  4.  In order to evaluate the proposed theoretical model, the comparison of theoretical, simulated and experimental thickness was performed. A micrometer with the accuracy of 0.01 mm was used to measure the thickness of the experimental sheet. In order to reduce the errors as much as possible, the formed parts are cut using wire cutting to obtain the the errors as much as possible, the formed parts are cut using wire cutting to obtain the measured sample, and then the average value of five measurements is taken as the actual thickness. The simulated thickness is extracted from the thickness fields based on LS-PrePost, which is solved using the finite element method (FEM). The theoretical thickness is calculated using the presented algorithm according to Equation (6). The corresponding thicknesses for the truncated cones are plotted in Figure 5 as a comparison.  Figure 5a shows the thickness profile of the truncated cone with a forming angle of 76° away from the central axis at x = 0 mm along the radial direction. The bottom surface of the truncated cone located within the range of x = 0-20 mm, which belongs to the unformed area, gradually moves downward as a whole, when the forming is carried out. The inclined wall located within the range of x = 20-40 mm is the formed area. The area over x = 40 mm, of which the outer edge is clamped by blank holder, belongs to the

Case 2: A truncated Pyramid Part
A truncated pyramid as a NURBS surface of degree 1 in each of u and v with 5 × 3 control points can be constructed as  The first derivatives over the designed pyramid surface are shown in Figure 7. The truncated pyramid, which was formed using aluminum sheet with thickness of 1 mm and axial step down 0.3 mm from Table 2, is presented in Figure 8. The first derivatives over the designed pyramid surface are shown in Figure 7. The truncated pyramid, which was formed using aluminum sheet with thickness of 1 mm and axial step down 0.3 mm from Table 2, is presented in Figure 8. The truncated pyramid, which was formed using aluminum sheet with thickness of 1 mm and axial step down 0.3 mm from Table 2, is presented in Figure 8. Figure 9 shows the comparison among experimental, simulated and theoretical predicted thickness distributions for the truncated pyramid part. The average thicknesses of the inclined wall in the formed area within the range of 10-30 mm are 0.801 mm and 0.816 mm, respectively, in the experiment and simulation. The results predicted by this presented method and calculated by the traditional sine law are both 0.80 mm for the truncated pyramid part with the forming angle of 53.2 • , which indicates that the result of the presented algorithm is completely consistent with the sine law for parts with a constant forming angle. In addition, compared with Figure 5, it can be seen that the smaller the forming angle, the greater the thickness reduction, and the thinner the final sheet, that is, the actual thickness is basically determined by the forming angle. The thicknesses of the  Figure 9 shows the comparison among experimental, simulated and theoretical predicted thickness distributions for the truncated pyramid part. The average thicknesses of the inclined wall in the formed area within the range of 10-30 mm are 0.801 mm and 0.816 mm, respectively, in the experiment and simulation. The results predicted by this presented method and calculated by the traditional sine law are both 0.80 mm for the truncated pyramid part with the forming angle of 53.2°, which indicates that the result of the presented algorithm is completely consistent with the sine law for parts with a constant forming angle. In addition, compared with Figure 5, it can be seen that the smaller the forming angle, the greater the thickness reduction, and the thinner the final sheet, that is, the actual thickness is basically determined by the forming angle. The thicknesses of the inclined wall region are basically the same for the truncated cone and pyramid with the forming angle of 53.2°.    Figure 9 shows the comparison among experimental, simulated and theoretical predicted thickness distributions for the truncated pyramid part. The average thicknesses of the inclined wall in the formed area within the range of 10-30 mm are 0.801 mm and 0.816 mm, respectively, in the experiment and simulation. The results predicted by this presented method and calculated by the traditional sine law are both 0.80 mm for the truncated pyramid part with the forming angle of 53.2°, which indicates that the result of the presented algorithm is completely consistent with the sine law for parts with a constant forming angle. In addition, compared with Figure 5, it can be seen that the smaller the forming angle, the greater the thickness reduction, and the thinner the final sheet, that is, the actual thickness is basically determined by the forming angle. The thicknesses of the inclined wall region are basically the same for the truncated cone and pyramid with the forming angle of 53.2°.  As shown in Figure 9, we can conclude that: (1) the experimental and simulated results are basically consistent with the theoretical prediction results; (2) the experimental thickness is basically the same as the original thickness in the unformed area and the experimental thickness is basically constant in the inclined wall area; and (3) the experimental thinnest thickness is 0.76 mm at 18 mm away from the center; the corresponding simulated and theoretical predicted thicknesses are 0.7923 mm and 0.8 mm, respectively. The deviations from the predicted thicknesses are 0.04 mm, and the corresponding errors are 5.26%.

Case 3: A Hemi-Ellipsoid Part
The construction of an ellipsoid as a NURBS surface with the degree of 2 in each of u and v is similar to that for the cone.
With 9 × 3 control points can be constructed as: A hemi-ellipsoid shown in Figure 10 was designed by:

Case 3: A Hemi-Ellipsoid Part
The construction of an ellipsoid as a NURBS surface with the degree o and is similar to that for the cone.

32
+ 24 + 15 = 1, 0 ≤ ≤ 15 The first derivatives over the designed hemi-ellipsoid surface are shown It can also be expressed as: The first derivatives over the designed hemi-ellipsoid surface are shown in Figure 11. The hemi-ellipsoid, which was formed using Al1060 sheet with original thickness of 1 mm and the step-down 1 mm shown in Table 2, is presented in Figure 12.  Materials 2022, 15, 1201 13 of 16 Figure 11. The first derivatives over the designed hemi-ellipsoid surface.
The hemi-ellipsoid, which was formed using Al1060 sheet with original thickness of 1 mm and the step-down 1 mm shown in Table 2, is presented in Figure 12.  Figure 13 shows the thickness distributions of the experimental part and the theoretical result for the hemi-ellipsoid part. The curvature of an ellipsoidal surface is variable. Within the range of x = 0-30 mm along the major axis and y = 0-24 mm along the minor axis, the thickness of the sheet located in the formed region is varied. In the transition area where the length of the major axis exceeds 32 mm and the length of the minor axis exceeds 24 mm, the final thickness decreases slightly, because the deformation occurs in this area. In the small area near the center, the final thickness is basically close to the original thickness, indicating that the sheet in the area is mainly translated downward during the forming process, and is not directly processed by the forming tool. The minimum thicknesses along the major axis and minor axis occur near the outer edge of the ellipsoid. It should be noted that the final thickness of the formed hemi-ellipsoid part at any position can be predicted by the presented algorithm in this paper. However, The hemi-ellipsoid, which was formed using Al1060 sheet with original thickness of 1 mm and the step-down 1 mm shown in Table 2, is presented in Figure 12.  Figure 13 shows the thickness distributions of the experimental part and the theoretical result for the hemi-ellipsoid part. The curvature of an ellipsoidal surface is variable. Within the range of x = 0-30 mm along the major axis and y = 0-24 mm along the minor axis, the thickness of the sheet located in the formed region is varied. In the transition area where the length of the major axis exceeds 32 mm and the length of the minor axis exceeds 24 mm, the final thickness decreases slightly, because the deformation occurs in this area. In the small area near the center, the final thickness is basically close to the original thickness, indicating that the sheet in the area is mainly translated downward during the forming process, and is not directly processed by the forming tool. The minimum thicknesses along the major axis and minor axis occur near the outer edge of the ellipsoid. It should be noted that the final thickness of the formed hemi-ellipsoid part at any position can be predicted by the presented algorithm in this paper. However, Figure 12. Formed hemi-ellipsoid part. Figure 13 shows the thickness distributions of the experimental part and the theoretical result for the hemi-ellipsoid part. The curvature of an ellipsoidal surface is variable. Within the range of x = 0-30 mm along the major axis and y = 0-24 mm along the minor axis, the thickness of the sheet located in the formed region is varied. In the transition area where the length of the major axis exceeds 32 mm and the length of the minor axis exceeds 24 mm, the final thickness decreases slightly, because the deformation occurs in this area. In the small area near the center, the final thickness is basically close to the original thickness, indicating that the sheet in the area is mainly translated downward during the forming process, and is not directly processed by the forming tool. The minimum thicknesses along the major axis and minor axis occur near the outer edge of the ellipsoid. It should be noted that the final thickness of the formed hemi-ellipsoid part at any position can be predicted by the presented algorithm in this paper. However, it is very difficult to calculate the thickness of the parts with the varying curvature by using the traditional sine law.
it is very difficult to calculate the thickness of the parts with the varying curvature by using the traditional sine law. As shown in Figure 13, we can conclude that: (1) the experimental and theoretical thickness profiles are basically consistent; (2) the final thickness differs significantly with the alteration of the forming angle; and (3) the experimental thinnest thicknesses are 0.580 mm and 0.465 mm, respectively, at 30 mm and 22.5 mm away from the center along the direction of major and minor axes; the corresponding theoretical predicted thicknesses are 0.6208 mm and 0.5106 mm, respectively. The deviations between the experimental and theoretical thickness are 0.0408 mm and 0.0456 mm, respectively, and the corresponding errors are 7.03% and 9.81%, respectively.
In order to have a further understanding of the presented model, the maximum deviation and root mean square error (RMSE) between the predicted, experimental and simulated thicknesses shown in Figure 5, Figure 9 and Figure 13 are summarized, calculated and listed in Table 3.  As shown in Figure 13, we can conclude that: (1) the experimental and theoretical thickness profiles are basically consistent; (2) the final thickness differs significantly with the alteration of the forming angle; and (3) the experimental thinnest thicknesses are 0.580 mm and 0.465 mm, respectively, at 30 mm and 22.5 mm away from the center along the direction of major and minor axes; the corresponding theoretical predicted thicknesses are 0.6208 mm and 0.5106 mm, respectively. The deviations between the experimental and theoretical thickness are 0.0408 mm and 0.0456 mm, respectively, and the corresponding errors are 7.03% and 9.81%, respectively.
In order to have a further understanding of the presented model, the maximum deviation and root mean square error (RMSE) between the predicted, experimental and simulated thicknesses shown in Figures 5, 9 and 13 are summarized, calculated and listed in Table 3. Table 3. The maximum deviation and RMSE between the predicted, experimental and simulated thicknesses. The following points can be seen from Table 3.

RMSE between Experiment and Simulation
(1) For all formed parts, the maximum deviation is only 0.0456 mm, and the corresponding percentage is 9.81%, which indicates that the deviation of prediction is basically small. It can be concluded that the thickness of theoretical prediction matches well with that of experiment and simulation within the acceptable deviation. (2) With the forming part deeper and the surface feature more complex, the prediction error will be larger. (3) The maximum thinning often occurs in the transition region between the formed area and the unformed area. (4) The maximum RMSE between the experimental and simulated thicknesses which occurs for the formed truncated pyramid is 0.0336 mm. The maximum RMSE between the predictions and experimental data, finite element results are 0.0499 mm and 0.0539 mm, respectively. The results show that the simulated thickness is closer to experimental thickness than that of the prediction, and the proposed model in this paper can effectively predict the thickness of parts in incremental sheet metal forming.

Conclusions
In summary, a general algorithm to predict and calculate the final thickness for arbitrary shape parts in incremental sheet metal forming based on a NURBS surface is presented. A detailed analysis is accomplished for the truncated cones, pyramid and ellipsoid parts. The experimental tests based on the customized deformation machine and finite element simulations using ANSYS/LS-DYNA are conducted. The comparison of the thickness of three kinds of formed parts shows that the maximum deviation and maximum RMSE are 0.0456 mm (9.81%) and 0.0539 mm, respectively. It is also found that the maximum deviation in the case of a constant forming angle is smaller than that of a variable forming angle. With the increase in depth and feature complexity, the deviation increases. The prediction results agree well with the experimental measurements and finite element results, which verifies that the proposed mathematical modeling approach can well predict the thickness distributions of the parts manufactured by using incremental sheet metal forming. Compared with the traditional sine law, this presented method can predict and calculate the thickness profile of arbitrary complex shape parts with variable curvature surfaces such as spherical and ellipsoidal parts in the design stage. This work provides a promising method to prevent thickness thinning during sheet metal forming.