Characterization of shape and dimensional accuracy of incrementally formed titanium sheet parts

Single Point Incremental Forming (SPIF) is a relatively new process that has been recently used to manufacture medical grade titanium sheets for implant devices. However, one limitation of the SPIF process may be characterized by dimensional inaccuracies of the final part as compared with the original designed part model. Elimination of these inaccuracies is critical to forming medical implants to meet required tolerances. In this study, a set of basic geometric shapes were formed using SPIF to characterize the dimensional inaccuracies of grade 1 titanium sheet parts. Response surface functions are then generated to model the deviations at individual vertices of the STL model of the part as a function of geometric shape parameters such as curvature, depth, wall angle, etc. The generated response functions are further used to predict dimensional deviations in a specific clinical implant case. The predicted deviations show a reasonable match with the actual formed shape and are used to generate optimized tool paths for minimized shape and dimensional inaccuracy. Further, an implant part is then made using the accuracy characterization functions for improved accuracy. The results show an improvement in shape and dimensional accuracy of incrementally formed titanium medical implants.


Introduction
Titanium is the material of choice in class III medical implants due to its biological inertness, strength, lightweight nature, bio-compatibility, and low-cost production [1]. Forming titanium into desired implant shapes within a specific time-frame is therefore of fundamental importance to clinical practice. To enable this, one relatively new manufacturing technique that has come forth is incremental sheet forming (ISF). Typically, in ISF, a hemispherical tool is used to deform a flat sheet in steps following a toolpath tailored to the geometry to be formed on a computer-numeric controlled (CNC) machine. ISF can be done in many different ways, as shown in Fig. 1. Variants include the use of a counter-support tool and use of dies, either full or partial. A number of efforts have been made to manufacture implants and supports for different parts of the human body using ISF such as the skull [2][3][4][5], knee [6], face [7], and ankle support [8] using ISF.
Process optimization efforts in the field of incremental forming have been focused on three key aspects: enhanced process limits, improved accuracy, and uniformity in sheet thickness [10]. The generation of intelligent toolpath strategies has been the key to improved part manufacture in all of these aspects. Process limits were enhanced using multi-step toolpaths as illustrated by Duflou et al. [11]. The use of helical toolpaths as shown by Skjødt et al. [12] and Cao et al. [13] helps eliminate scarring caused by the tool stepping down in a contouring toolpath. Bambach et al. [14] proposed the use of generating toolpaths on compensated part geometries to improve the accuracy, while Li et al. [15] have shown that a multi-stage process can result in improved thickness distribution in incrementally formed parts.
Despite a number of efforts to make medical implant shapes using ISF [2][3][4][5][6][7][8], making these parts with high accuracy has been a problem. Even though Behera et al. [2] tried forming implant shapes with high forming angles with improved accuracy, the part made of titanium grade II failed. The work of Göttmann [3] illustrated the ability to form implants with a maximum deviation at the edges just less than 2 mm using two point incremental forming (TPIF). Despite these efforts, no definite characterization of freeform surfaces and titanium implants made by ISF is currently available. Some accuracy characterization techniques are available. These include the use of multivariate adaptive regression splines (MARS) within a feature-based framework for predicting the behavior of simple features and feature interactions [16] and a local geometry matrix to predict spring back [17].
Earlier, Verbert et al. proposed the use of a feature-based approach to form parts with high accuracy [18]. In this approach, four basic geometric features were identified based on principal curvatures in the part, viz.: planar, ruled, freeform, and ribs. While this study provided an overall schematic for carrying out part compensation in order to optimize the toolpaths needed to form the part accurately, it did not cover the specific steps of identifying the relevant geometric parameters and error correction functions necessary for each feature type. Likewise, the work of Behera et al. [16,19] was limited to studies on planar and ruled features and aluminum and low carbon steel alloys. Micari et al. [20] and Essa et al. [21] outlined various process strategies to improve accuracy in incremental forming. An in-process online correction strategy was laid out by Rauch et al. [22] which was limited to correcting the depth accuracy of the parts. Lu et al. [23] showed that the use of critical edges in generating toolpaths can improve surface quality, forming time, and geometric accuracy in specific cases. Despite a number of efforts, these works did not provide any methods for predicting inaccuracy in freeform implant shapes. Furthermore, titanium is a material not covered by current accuracy models.
To overcome the limitations of the current accuracy characterization techniques, an effort is made in this work to generate accuracy response surfaces for freeform shapes. This is done by studying the accuracy behavior of ellipsoidal shapes formed using single point incremental forming (SPIF), which is one of the process variants of ISF. The major and minor axes of the ellipsoids are used as parameters in the characterization models generated using MARS. To account for the effect of presence of multiple features in a part, a mixed model using an index generated from principal curvatures of points in the part is also proposed. These models are then used to predict accuracy behavior of new implant geometries and the predicted behavior is then used to compensate for the inaccuracy of formed parts. All parts made in this research are formed using uni-directional contouring tool paths with a uniform scallop height between successive contours, where scallop height is a parameter that determines surface quality as shown in Fig. 2.

Accuracy characterization methodology
The accuracy of a part formed by incremental forming is typically determined by measuring the same with metrology tools such as a laser scanner or a coordinate measuring machine. After carrying out this measurement, a point cloud showing the a use of counter support or two point incremental forming (TPIF), b full die, and c partial die [9] representing the part can be generated and this point cloud can be meshed and then compared with a mesh representing a nominal model obtained from the computer-aided design (CAD) model of the part. The measurement process with a point cloud gives the coordinates of the formed part, which is a large data set, often as high as 100,000-500,000 points for a single part. Hence, the deviations with respect to the CAD model also form a large dataset. These deviations can thus be modeled as function of geometrical parameters for points on the surface of the nominal model. However, to do this, a robust, statistical tool for high dimensional data is needed.
Some tools that are currently available and used in modeling of high dimensional data include generalized additive models (GAM), multivariate adaptive regression splines (MARS), minimax probability machine regression (MPMR), and least square support vector machine (LSSVM) [24]. Of these, MARS has already been used to generate models of accuracy behavior in planar and ruled features made with specific materials such as AA 3103 and DC01 [19]. This technique was also applied in the current study. Figure 3 shows a schematic of the accuracy characterization methodology. Training set CAD models are used within a CAM software such as Siemens Unigraphics NX to generate an uncompensated toolpath. This toolpath may need to be post processed into the machine tool format and then fed to a CNC machine tool used for incremental sheet forming. The formed part from the ISF machine tool is then scanned with a metrology tool such as a co-ordinate measuring machine (CMM) or a laser scanner to generate a point cloud, which is then fed to a metrology software such as GOM Inspect. This software compares the training set CAD model to the mesh generated from the point cloud generated from the CMM to yield an accuracy data file. This accuracy file is fed to a custom STL processing software for incremental sheet forming built in Visual C# to carry out the current study. This software detects features in the part, using the criteria discussed later in Section 2.3, and also calculates geometrical parameters such as wall angle, principal curvatures, etc. which are discussed in Section 2.2 depending on the type of detected feature. The accuracy data is then linked to the training set CAD model using KDTrees. KDTrees are multidimensional binary search trees which can carry out quick spatial comparisons between two data sets using an associative searching technique [25]. Each vertex in the nominal training set CAD model is linked to a vertex on the measured CAD model closest to the nominal vertex. This is necessary as the geometrical parameters are calculated for points on the nominal model while the accuracy data file consists of deviations for measured data points and hence, the linking process links accuracy data to geometrical parameters. This linked data set is then exported to a data file which is fed to the statistical software "R" for generating accuracy response surface using MARS.  In this study, MARS models are generated as a sum of basis spline functions which are chosen using a forward pass and backward pruning step. The basis functions take one of three forms: (a) constant, (b) hinge functions of the type max(0, x-k) or max(0, k-x) where k is constant called the "knot", (c) product of hinge functions. To generate the MARS models, the statistical package "Earth" available within the software "R" is used [26]. It was important to find out suitable parameters to characterize the accuracy, and this is described in Section 2.1.

Model parameters
To build regression models, geometric and process-related variables need to be selected that potentially affect the accuracy behavior of the formed part. In prior work for planar and ruled features, it has been shown that the distance to the feature borders affects the accuracy behavior [26]. In the case of freeform features, the obvious problem is the lack of a generic defining distance in the horizontal plane of the backing plate as freeform surfaces do not have an immediately obvious symmetry and as such a defining distance can be problematic in characterizing the accuracy. However, when we consider the case of a cranial implant, we observe that the shape of the implant can be thought of as being close to that of an ellipsoid (Fig. 4a), characterized by a major axis and minor axis. This observation was further supplemented by experimental results, discussed later in Section 3, which showed that the accuracy in the direction of the major axis was different from the accuracy in the direction of the minor axis at a specific cross-sectional depth. It may also be noted that for ruled features, the maximum principal curvature is a model parameter, while the minimum principal curvature is zero with a tolerance value and hence, not a model parameter. In freeform surfaces, both principal curvatures have a finite value and hence, both can be included as parameters in the modeling. The accuracy behavior of ruled features that have close resemblance with a cranial implant is modeled by designing extruded drafts of ellipses, as shown in Fig. 4b.
The following parameters are thus used to predict the 3D deviation from the CAD model at a point: 2. total vertical length of the feature at the vertex (d v =(A+B) in Fig. 4), 3. major axis length on slicing the feature at the depth of the point (l max =C) 4. minor axis length on slicing the feature at the depth of the point (l min =D) 5. maximum curvature at the point, k max 6. minimum curvature at the point, k min 7. wall angle at the vertex (in radians), α 8. angle of the tool movement with respect to the rolling direction of the sheet (in radians), ω Of these, the minimum curvature for ruled surfaces by definition is zero with a tolerance and hence is not a parameter in the modeling for ellipse drafts.

Experimental details
A hemispherical tool of radius 5 mm was used for all the tests along with a feedrate of 2 m/min. A soft low melting-point paste lubricant, Rocol RTD Compound, was applied during the process. Four ellipsoids were used for training the models for freeform features with major and minor axis diameters 110× 60, 110×70, 110×90, and 90×60 (all dimensions in millimeters). Likewise, twelve ellipse drafts representing ruled features were made using the same major and minor axis diameters and wall angles of 15°, 30°, and 45°for each diameter combination. All parts were made in grade 1 titanium alloy of thickness 0.5 mm. A backing plate with an elliptical cross-section corresponding to the dimensions of the top contour of the part was used for each test. A contouring tool path with constant scallop height of 0.05 mm was used for forming the parts. The formed parts were unclamped and measured with a 3D coordinate measuring machine to generate point clouds representing the formed part shape. In the current work, the sheets used are thin and so the deformation on unclamping is significant. Hence, it is important to develop models for the net effect of deviations due to plastic deformation while forming and deviations due to unclamping.

Feature detection thresholds
The file format used within this research consists of triangulated representation of the part's surface known as the STL (stereolithography) format. The detection of features relevant to incremental forming requires the segmentation of the part's surface based on curvature calculations performed at individual vertices. These curvature calculations are done following the procedure outlined by Cohen-Steiner [27,28]. The curvature tensor at a vertex v can be calculated as: where |A| is the surface area of the spherical zone of influence of the tensor and β(e) is the signed angle between the normal vectors of the STL facets connected by the edge e. The weight for the contribution by each edge is given by the factor e∩A. The normal at each vertex is estimated as the eigenvector of Λ(v) evaluated by the eigenvalue of minimum magnitude. The other two eigenvalues, k min and k max , provide estimates of the minimum and maximum principal curvatures at v. Using these principal curvatures, the part can be split up into features using the following classification criteria [29]: Planar feature: k p min = 0±ε p and k p max = 0±ε p , where ε p is a small number that can be tuned for identifying planar features Ruled feature: k p min = 0±ε r and k p max = X, where X is a positive non-zero variable. Another possible case is where k p min = X and k p max = 0±ε r , where X is a negative non-zero variable. ε r is a small number that can be tuned for identifying ruled and freeform features. Freeform feature: k p min = Y±ε r and k p max = X±ε r , where X and Y are non-zero variables such that X ≤ ρ max and Y≥ρ min , where ρ max and ρ min are threshold values for distinguishing freeform and rib features. Rib feature: k p min ≤ρ min and/or k p max ≥ρ max Figure 5 shows examples of these features. A ruled feature is defined by a generatrix curve that is swept along a directrix line to generate the surface, as in the cone shown in Fig. 5b. A double curved surface where a generatrix is swept along another curve creates a freeform surface, as shown in Fig. 5c. It is important to set appropriate thresholds so as to get a usable segmentation of the part for building response surface models for accuracy.

Model for freeform ellipsoidal parts
The accuracies of the formed ellipsoids (shown in Fig. 6) are listed in Table 1. It can be seen that changing the diameters of the major and minor axes affects the accuracy of the part. The smallest part shows the highest over forming. This can be attributed to the low wall angles in the part, which are usually responsible for over forming. The largest part with the highest wall angles (top contour of 110×90 mm) shows exclusively under forming. This is due to two reasons: (a) ellipsoids are essentially positive curvature freeform surfaces and positive curvature tends to under form [9], (b) the biggest part has high wall angles in the initial forming steps and high wall angles in a positive curvature region are known to under form [9]; this leads to the lower depths also showing under forming being a continuation of the top surface.
The MARS model was trained with accuracy data from these tests resulting in the following model. where, e f =deviation at STL vertex of a freeform feature and the remaining abbreviations are the same as in Section 2.1.

Model for ellipse draft ruled surfaces
Nine of the 12 ellipse draft ruled surface parts were chosen as training sets for generating the models. These were the parts made with major and minor axis diameter combinations of 110 × 60, 110 × 90, and 90 × 60. The remaining three were used later for verifying the model's validity for a new combination of major and minor axis diameters. Table 2 lists the accuracies of these parts. It is observed that at low wall angles, the formed part shows large deformations after unclamping, which is absent at higher wall angles, 30°and 45°. In general, where e r =deviation at STL vertex of a ruled feature.

Generalized model for a part with mixed curvatures
For a part with mixed curvatures, feature detection with thresholds used for detecting the ellipsoidal parts as freeform features results in the part being detected as a mixture of ruled and freeform surfaces as shown in Fig. 7a. This creates small ruled features surrounded by a larger freeform feature. To carry out compensation of the part, the vertices in the ruled features would then be compensated with a different error correction function than the freeform feature. This would introduce discontinuities in the predicted and compensated surfaces. The problem is observed even if the thresholds are changed to those needed to detect the part largely as a ruled surface by using the same thresholds that were used for detecting the ellipse drafts, as shown in Fig. 7b. Another technique to correct the part would be to use a network of features and use a single compensation function for the feature interactions [30]. However, modeling this feature interaction would require more experimental tests and also need to cover the different locations where the ruled  features occur. Furthermore, due to the small surface area of the ruled features, the feature interaction error prediction equations may not capture the correct accuracy behavior. Hence, a new approach to model the accuracy of parts with mixed curvatures was used in this study. In this approach, the part is detected as a single feature and an index that evaluates the extent to which the feature is a freeform feature is calculated. This index is evaluated as: where δ is the extent to which the feature is a freeform feature, κ min part ð Þ is the mean of minimum curvatures of the vertices in the part, κ min ellipsoids ð Þis the mean of the minimum curvatures of the vertices in the ellipsoid training sets used for generating Eq. (2), and κ min ruled ð Þ is the mean of the minimum curvatures of the vertices in the ellipse draft training sets used for generating Eq. (3). Now, the deviations at individual vertices of the part are calculated as : where, e m is the predicted deviation of a vertex in a mixed feature consisting of both ruled and freeform surface vertices, and e f and e r are calculated using Eqs. (2) and (3) (5) was then also used to predict the shape of the implant to find if it improves the prediction provided by Eq. (2). The geometry of the cranial implant is shown in Fig. 8.

Model validation results
The results of the validation are presented in Table 3. Here, the predicted meshes are compared with the mesh representing the point clouds from the actual formed part from experiments and the prediction error is thereby found as the deviations between these two meshes. It is seen that the training test cases show a mean deviation close to zero for three of the four ellipsoids, while the implant shows a mean deviation of −0.82 mm. The error in the prediction of the formed shape of the implant is [−1.0 mm, 1.0 mm], which is not as good as the prediction for the ellipsoids (see example prediction for ellipsoids in the sections shown in Fig. 9). The reason for the Fig. 8 Cranial implant geometry in isometric view, x-section along AA′ and y-section along BB  Fig. 11a, is that the ellipsoid is still not a perfect representation of the curvatures in the cranial plate and the major axis and minor axis dimensions obtained for the cranial plate are only an approximation of the material deformation in the case of an ellipsoid. However, it would be reasonable to say that the model in Eq. (2) is a good starting point for prediction of positive curvature freeform surfaces and this model can possibly be improved further either by choosing different geometrical parameters for the model or using more training sets for more complex geometries. This would however reduce the robustness of this methodology. The predictions for the ellipse drafts was done for parts with major axis diameter of 110 mm and minor axis diameter of 70 mm and wall angles of 30°, 40°, and 45°, which constitutes parts outside of the training sets used for generating the model. It is seen that Eq. (3) predicts the accuracy of ellipse drafts with mean deviations of −0.77, −0.12, and 0.29 mm respectively for parts with 30°, 40°, and 45°wall angles. Figure 10 shows the predicted sections for the part with wall angle of 40°.
Using Eq. (4), δ is evaluated as 0.73 and the prediction accuracy for the cranial implant improves to [−0.85 mm, 0.14 mm] by using the mixed model (see section in Fig. 11b).

Compensation technique
The compensation of the parts is carried out by translating individual vertices in the nominal CAD model of the part normal to the part geometry by a magnitude equal to a compensation factor multiplied with the predicted deviation at the point. This follows the strategy outlined by Bambach et al [14]. Three different compensation factors were tried out for compensation using Eq. (2), 0.7, 1 and 1.3 and the best among these factors, 0.7 was used for compensation using Eq. (5).

Accuracy of compensated implant
The results of the compensation are presented in Table 4 along with the result for a part made without compensation (compensation factor 0). Using Eq. (2), it can be seen that the part with the best accuracy is realized with a compensation factor of 0.7. A color plot of the part accuracy and that of the   Fig. 12 to illustrate the utility of using the model in Eq. (2) in forming parts with higher accuracy. It can be seen that with increasing compensation, the over forming in the part systematically increases from 0.57 mm with a compensation factor of 0.7-0.96 mm with a compensation factor of 1.3. In contrast, the under forming decreases from 0.53 mm to 0.25 mm. This indicates that increasing the compensation factor creates an over-compensated surface, which is outside the nominal CAD model, thereby enhancing the over forming and reducing the under forming, which is caused primarily due to spring back of the material.
It is also noted that in the uncompensated part, the color plot in Fig. 12a reveals that the zone of significant under forming seen in reddish-yellow tone corresponds to the area with the highest wall angles. This is on expected lines, as previous experiments have shown that low wall angle regions show either low under forming or over forming, while high wall angle regions are prone to spring back of the material [9]. The color plot shown in Fig. 12b shows that this zone of under forming is well compensated using the compensated tool path generated by the prediction given by Eq. (2), as evident by the green patch in the center of the part.
Using Eq. (5) and a compensation factor of 0.7, the mean deviation improves from 0.05 to 0 mm, although the maximum deviation in the part increases from 0.53 to 0.74 mm. Figure 13 shows the manufactured part and a color plot of the accuracy. From this result, it can be concluded that the better predictions obtained by using Eq. (5) help in reducing the average deviations in the part geometry. It is again evident from the color plot that the zone of high wall angles in the center of the part is mostly well compensated and is in green.
There are two small areas in reddish-yellow tone in this plot, which shows that the prediction was not very accurate in these regions. However, the improved mean deviation indicates that the mixed MARS model in Eq. (5) is able to account for the variations in the principal curvatures in the part better than the model specifically built for freeform features in Eq. (2).

Limitations and discussion on the developed models
It is noteworthy that within this study, only positive curvature parts were considered for developing the response surface models. This may pose a limitation for directly applying the developed models for parts that also constitute of negative curvature regions. In parts with a mix of positive and negative curvature freeform regions, a model for negative curvature features can be developed and using the approach provided in this work to develop a mixed-MARS model, a similar mixed model can be generated.
Yet another limitation of the current models is for parts with wall angles close to the failure limit for the material. In the current study, the parts considered were far away from failure  and hence, the compensated surfaces were also within the failure limits. As shown in earlier works [26,30], when the compensated surface is beyond the failure limit, toolpaths generated using MARS models may lead to sheet thinning and failure. In such cases, it is useful to adopt other tool path strategies such as morphing [31]. Furthermore, this method will work well when the part to be formed has a similarity to that of an ellipsoid. While the ellipsoidal shape was seen to be a very generic method of characterizing a number of cranial implant shapes, the application of these models to parts which are distinctly different from an ellipsoid should be done with caution. This is particularly important when sections taken at increasing depths do not show a decreasing trend in the major and minor axis lengths. While this limitation may be kept in mind, it is noteworthy that the developed approach within this study opens up the possibility of developing a generic method of taking complex shapes and observing similarities with known geometrical shapes such as ellipsoids, hyperboloids, cones, etc. and using the basic characterizing dimensions in developing suitable geometrical parameters that will then be used to develop accuracy response surfaces using MARS.
The validation work on the applicability of the mixed MARS model was limited to one substantial case study of the cranial implant. However, the individual models for ruled and freeform features have been tested to work well in predicting feature accuracy within reasonable limits for the maximum and minimum deviations in a number of test cases. These include the cases shown in Table 3. It was not necessary to include material properties in the current study, as this study only focused on grade 1 titanium sheet parts for cranial implant applications. However, the MARS models capture the effect of the mechanical properties such as tensile strength, which affect the spring back behavior of sheet materials during incremental forming.
In addition, it may be noted that the developed models within this paper were for sheets of thickness 0.5 mm and maximum dimensions limited to 110×90×35 mm. For parts outside these dimensions, simple extrapolation of model predictions may not yield accurate results. It was not an objective of this current study to explore the effects of material, size, and sheet thickness effects but to develop a methodology for freeform titanium sheet parts that would be applicable mostly to cranial implant applications and parts with regions of intermediate principal curvatures between two feature types.

Conclusions
In this study, a method to predict the accuracy behavior of incrementally formed titanium sheet parts was developed. MARS were used to develop models for ruled and freeform features that predict the accuracy at individual points in the STL file of the part to be formed. The models for freeform features are based on four training sets formed as ellipsoids while the models for ruled features were developed by varying the wall angles and major and minor axis dimensions at the top of the part. The major and minor axes of the parts at the depth (z-axis) cross-section of individual points in the part are used as variables in the model. In addition, both maximum and minimum curvatures are used as predictors for the freeform features model.
The results show that the generated model for freeform features is reliable in predicting ellipsoids from the training set with mean deviations for the prediction accuracy varying between −0.23 and 0.04 mm, while the prediction accuracy deteriorates for a new part such as a cranial implant. The model for ruled features predicts the accuracy of a new ruled feature with mean deviations varying between -0.77 and 0.29 mm, for the specific validation test cases performed. A new mixed model based on curvatures in the part was proposed, and it was verified to improve the average prediction for the cranial implant by 0.33 mm.
It may be however noted that the models generated using the ellipsoids for freeform surfaces should not be directly applied to any generic shape, unless an approximation in terms of major and minor axis is immediately apparent from the model geometry. Furthermore, the ellipsoids are positive curvature features, while a freeform part such as a human face model may consist of negative curvature regions where the models provided in this study will not be valid. However, the modeling strategy provided in this paper for parts with mixed curvatures using an index as shown in Eq. (4) has been shown to be a promising strategy.
Using the model for freeform features, part compensation was carried out for a grade 1 titanium implant part and the accuracy of the formed part was seen to improve vis-à-vis the uncompensated part with the best results achieved with a compensation factor of 0.7. Compensation using the model reduced the maximum deviations from 1.02 to 0.53 mm. Using the mixed model and the same compensation factor of 0.7, improvement in part accuracy was also realized in terms of both reduced maximum and mean deviation compared to the uncompensated part.
Further studies following this research can include improving the model with better predictors or a new prediction technique such as GAM. It would also be useful to study the effects of material properties and sheet thickness on the accuracy response functions, as material properties and sheet thickness can vary from one batch to another even for the same material and this will affect the plastic deformation and spring back resulting from forming the part. Systematic prediction of the accuracy of freeform titanium part surfaces using numerical methods such as finite element analysis is another potential area of investigation, where the results from the studies carried out within this research can be compared to numerical predictions.