A fully parametric toolbox for the simulation of single point incremental sheet forming process: Numerical feasibility and experimental validation

Single point incremental forming (SPIF) is a sheet metal forming process that allows man- ufacturing components without the development of complex tools in comparison with stamping process. This work is dedicated to the development of SPIF for microparts (shape or detail) and for thin metal sheets (less than 1 mm). This paper focuses on the deﬁnition of a numerical tool- box to simulate this process at these scales. Forming of a pyramidal shape of 4 mm in height from a copper sheet with an initial thickness of 0.21 mm is suggested. The complete meth- odology is proposed and numerical results are presented in terms of global geometry (shape and section proﬁles), thickness evolution and forming forces. Equivalent experimental tests are carried out to validate the numerical approach. Different ways of evolution are provided at the end of this work.

Single point incremental forming (SPIF) is a sheet metal forming process that allows manufacturing components without the development of complex tools in comparison with stamping process.
This work is dedicated to the development of SPIF for microparts (shape or detail) and for thin metal sheets (less than 1 mm). This paper focuses on the definition of a numerical toolbox to simulate this process at these scales. Forming of a pyramidal shape of 4 mm in height from a copper sheet with an initial thickness of 0.21 mm is suggested. The complete methodology is proposed and numerical results are presented in terms of global geometry (shape and section profiles), thickness evolution and forming forces. Equivalent experimental tests are carried out to validate the numerical approach. Different ways of evolution are provided at the end of this work.

Introduction
Single point incremental forming (SPIF) process is interesting both industrially and scientifically. In the first case, sheet metal components can be manufactured without specific tools using a CNC milling machine. This kind of process can be produced complex parts in small batch or for single part. In particular, Jeswiet and Hagan [1] used SPIF as a rapid manufacturing process to custom-made parts. However, this advantage is limited by the important thinning of the sheet, the occurrence of defects and an important working time as mentioned in the complete review offered by Jeswiet et al. [2] on asymmetric SPIF.
For scientists, incremental sheet forming (ISF) exhibits local effects and needs then new characterization methods. Filice et al. [3] define a stretching test to investigate material formability for straining conditions in the range between uni-axial and bi-axial stretching. From these tests, they suggested, a forming limit diagram for an aluminum alloy different from the conventional one. Ham and Jeswiet [4], have built forming limit diagrams by using a Box-Behnken design of experiments and response surfaces method based on the two following criteria: maximum forming angle and effective strains.
In Jeswiet et al. [5] and Duflou et al. [6], the formability is determined by considering the evolution of forming forces during the process. This approach is based on force measurements during the production of conical part. Typical curves are reported and influences of some process parameters are revealed. Ambrogio et al. [7] used this approach to evaluate a ''spy variable'' for defining a correction strategy to prevent failure during process.
Numerical simulations are very useful to develop manufacturing processes (feasibility, optimization). Therefore, numerical simulations based on the finite element method have been carried out for developing ISF process. Henrard et al. [8][9][10] have performed studies to propose the best ways for numerical modeling to predict the process correctly. Malhotra et al. [11] have investigated the use of several material models and element formulations to simulate SPIF. From these studies, it has been shown that element formulations (solid element), integration algorithms (transient dynamic explicit algorithm), material models and contact algorithms are the most influent parameters.
A second way to improve SPIF by virtual forming is based on toolpath strategy. Malhotra et al. [12] have proposed an automatic 3D helical path generator for the simulation of SPIF. Yamashita et al. [13] have conducted numerical simulation with different strategies and parameters on thin sheet.
Time simulation is often much longer than in the case of the other forming processes. This observation is also true for the simulation time, and Robert et al. [14] have proposed numerical techniques to reduce it by modifying the classical resolution of material behavior. In the study of Yamashita et al. [13], time simulation is also reduced by modifying the value of density which is equivalent to the use of a mass scaling algorithm.
Recent studies focused on the development of incremental sheet forming for microparts or for thin sheets [15][16][17][18]. The first study was due to Saotome and Okamato [19], and consisted in the development of a specific SPIF device for realizing 3D microparts in a scanning electron microscope. The main goal of the present paper is the development of a parametric toolbox for the simulation of micro-SPIF. In first, the considered material and the experiments used to validate the numerical approach are presented. Secondly, the fully parametric toolbox is defined and the numerical algorithm and methods used to simulate micro-SPIF are described. Finally, comparisons between numerical results and experimental data are discussed.
This paper shows the ability of the proposed numerical approach to predict SPIF process for manufacturing defect-free microparts.

Material
For this study, the selected material is a FPG copper sheet with an initial thickness of 210 lm and an initial hardness of 124 HV (known as half hard condition state). Its composition has been studied by Gréban et al. [20] and is given in Table 1.
The initial grain size is 10 lm.
The FPG material is used for the realization of lead frames (metal supports of electronics components). This kind of low alloy copper guarantees a high electrical conductivity required for the best performances of lead frames.

Experimental device
To validate the numerical simulations, experiments have been conducted. For this purpose, a dedicated apparatus has been designed and realized. The representation of this tooling is done in Fig. 1a. It is composed of a fixed die support, a modular die, a fixed blank holder clamped with the die by screws and the forming tool (end ball tool).
The modular die allows to define different shapes depending on the geometry of the part to be produced. In particular, it limits the non-desired bending obtained on the base of the part by specifying the nearest contour of the final part.
The forming process is performed by the use of a CNC machine tool (3-axis milling machine). The modular device is linked to a 4-axis dynamometer to get the forming forces during the process: three forces in CNC-axis and the torque along the z-axis. This set is clamped on the machine table as shown in Fig. 1b by bolting. The dynamometer acquisition frequency is set at 1 kHz. This approach has already been used in the case of conventional incremental sheet forming by Ambrogio et al. [7] and Duflou et al. [6].

The shape of the part
A pyramidal shape is proposed to investigate the single point incremental sheet forming of thin sheet and microparts. The geometry definition is illustrated in Fig. 2.
This shape has been chosen because there is alternation between the x and the y directions during the forming process. Therefore it will be possible to determine the influence of the tool position on the forming forces.
The geometrical parameters used in this study are given in Table 2. The tool radius is equal to the corner radius R.
In the following, it will be shown that the draft angle a is a discriminant parameter to choose a forming strategy.

Forming strategies
Different strategies are possible to produce the part. For the pyramidal shape, two approaches are considered: the constant Z-level and helical paths (Fig. 3).
The first strategy follows the internal contour of the part resulting from the intersection of the shape by perpendicular planes to the spindle axis (so parallel to XY plane) located at different z-position (Fig. 3a). Each outline is offset of a gap d 2 in z-direction and of a gap d 1 in x and y direction. The larger profile is done at the beginning of the process that ends by the smaller one. The parameters d 1 and d 2 are linked together by the relations The second approach consists in performing the shape by a helical path (Fig. 3b). By starting from one corner of the pyramid, the tool moves on the theoretical shape by coupling the 3-axis displacements according to the following relations   Table 2 Geometrical parameters of the pyramidal shape.
To ensure the spindle integrity, the tool rotates with a constant speed rate X and moves with a constant feed rate f.
To compare the two strategies, helix angle u and d 2 -increment are related with Eq. (3) by considering the length L i equals to the basis length l 1 .
The different parameters used in this study are given in Table 3.
The experimental results for these two strategies are illustrated in Fig. 4.
As shown in Fig. 4, these geometrical and forming parameters, with an initial grain size of 10 lm, discriminate the forming strategies. In the case of the helical path (Fig. 4a), a safe part in term of geometry is obtained. For the equivalent constant Z-level strategy, crack occurred before the desired part is not obtained. From experimental and numerical points of view, the influence of various material, geometrical and process parameters as the initial grain size (size effects), the forming strategies and tools geometry on the forming forces, defaults predictions, shape accuracy, will be investigated from these two tests.

Automatic mesh for parametric studies
One of the most important characteristics of the ISF process is associated with the use of non-specific tools (die and forming tool). From the numerical point of view, this specificity enables: -to build a parametric representation of the process; For that purpose, a dedicated toolbox, programmed in MATLAB Ò language, has been developed to create the parametric model (mesh, boundary, load and initial conditions, material behavior),  Table 3 Forming parameters of the pyramidal shape.

4
-to run the simulations (forming and springback) with LS-DYNA Ò software [21], -to post-process the results (deformed mesh display, springback and thickness measurements, forming forces control).
The complete loop is described in Fig. 5. This parametric toolbox also defines the experimental path in CNC language machine by the way of parametric paths (mathematical definition) or CAM files (APT language). The advantage of creating a unique file for the toolpath definition is the insurance to perform experiments and simulations with exactly the same trajectories.
Experiments and simulations are then compared by the use of a specific post-processing toolbox also developed in MATLAB Ò language. This last toolbox permits the calculation of specific results as well as conducting sensitivity analysis, identification and optimization procedures.

Geometrical discretization
The finite element method is chosen to simulate the process. Due to large shear deformations, the blank is meshed with solid elements. Fully integrated eight nodes solid elements are used to get as much information as possible in the thickness. Three elements in thickness are considered. For the present study, 120 elements are imposed in the length and the width to discretize the blank with a total of 43,200 solid elements.
Each tool (forming tool, die and blank holder) is meshed with rigid shell elements (quad elements). The element size of each tool is considering with the smallest value of the die radius and a number of five elements in this radius is imposed to ensure a good geometrical representation of the surfaces.
In order to decrease the simulation time, an initial study was carried out to limit the blank dimensions. As the deformations are located at the contact between tool and sheet, the elastic field does not propagate in the entire blank, and so the dimensions can be reduced, by half of their initial values at least. As a consequence, the initial length and width are assigned to 17 mm.
In this model, the boundaries of the blank are considered to be clamped. It results that all the degrees of freedom of each nodes on the blank boundary are set to be null during simulation.
The complete mesh for this study obtained automatically from the toolbox is presented in Fig. 6.

Numerical algorithms and clustering
To perform robust and predictive simulations, it is crucial to choose the numerical algorithms and the associated numerical parameters. The explicit algorithm is chosen for time integration according to its robustness provided that the stability condition is respected. The classical mass scaling algorithm commonly used to increase time step is not applied in this study.

5
The choice of the contact algorithm is most often critical in the simulation of this type of process. For this reason, the classical penalty method proposed in LS-DYNA Ò by Hallquist [22] is utilized. To ensure the stability condition and a good repartition of nodal forces, the rigid parts are meshed elements of size equivalent to those of the blank.

Virtual simulation time
The stability condition of the dynamic explicit algorithm limits the time step size. In principle, this type of algorithm is dedicated to short transient simulations. Nevertheless, it can be carried out to simulate highly nonlinear problems, including problems of multiple deformable bodies in contact. It is this last feature that leads to implement an explicit approach to simulate the incremental sheet forming.
However physical time cannot be used for the simulation time because the actual run time of the process is significant. In the case of helical path, the experimental execution duration is 25 s which is too important to run with explicit approach. So it is necessary to find a virtual simulation time that will leads to numerical results similar to the experimental ones. This adaptation must not generate kinetic energy that must stay negligible compared to internal energy. But this condition doe not overcome the dynamic effects.
As the material behavior is not time dependent, the simulation time can be chosen as small as possible, but dynamic effects may occur which have no physical meaning. Therefore several simulations were conducted with different duration to find the smallest one which does not introduce numerical dynamic effects. Finally a virtual simulation time of 0.2 s, equivalent to a tool feed rate of 1300 mm/s, has been selected.
Due to the important size of the model and the stability condition, the simulation time is important. The massively parallel processing (MPP) release of LS-DYNA Ò is utilized with 20 processors for reducing the simulation time.

Behavior law
In relation with the studies of Touache et al. [23], the mechanical behavior of the copper alloy is considered through an elastic-plastic law with isotropic hardening. This copper alloy exhibits an isotropic elastic law and no strain rate sensitivity. According to these considerations, the following behavior law is proposed where _ k and u are respectively the plastic multiplier and the potential function. The plastic multiplier must satisfy the following Kuhn-Tucker's conditions _ k P 0 and _ k/ ¼ 0 ð12Þ The internal variables R and p are respectively the isotropic hardening variable and the equivalent plastic strain. The fourth order tensors C and P are the classical elastic operator and the deviatoric projector. The r y , K and G parameters are respectively the yield stress, the bulk modulus and the shear modulus. The set of material parameters was obtained by tensile tests and by ultrasonic tests. This set of identified parameters is summarized in Table 4. The density of this material is equal to 8500 kg m À3 . The friction law chosen to simulate the tribological behavior at the interfaces between tools and blank is the Coulomb's friction law. Lubricant (oil and water) is used during experiments and Table 4 Identified material parameters for FPG alloy. we considered only a static friction coefficient f s = 0.2 in the simulations to take into account steel/copper contact. This choice derives from studies on the influence of friction on forming forces level and will be explained further.

Results and discussions
For the validation of the numerical model, comparisons between experiments and simulations were only helical strategy was considered due to premature failures with the constant Z-level strategy (Fig. 4b). To perform these comparisons the temporal evolutions of the experimental and numerical results are synchronized with reference to the percentage of the forming cycle rather than actual or virtual time.

Numerical forming
As mentioned above, simulations were performed with the MMP version of LS-DYNA Ò using 20 processors. The total computational time took 4 h. The mesh evolution in function of forming cycle is presented in Fig. 7.
The ability of such finite element code to simulate the incremental forming parts of small dimensions or low thickness is proved. The study of the energies evolutions (total energy, kinetic energy, internal energy and contact energy) shows that the initial assumptions are valid: the kinetic energy is well negligible compared to the internal energy. The values are positive and demonstrate that no numerical energy is introduced. The sum of energies gives the total energy. From these considerations the simulations can be compared with experiments.
From the final mesh corresponding to 100% of the process cycle, the maximum value of the equivalent plastic strain reaches over 240%.

Springback simulation
The geometrical comparisons may be carried out provided that the numerical and experimental stress states are the same. Springback prediction is done in a second simulation step by the use of an implicit algorithm (LS-DYNA in implicit mode). This step is associated to the elastic stress relaxation in the part when the tools are withdrawn. Springback result is illustrated in Fig. 8.
Its influence on the final shape of the part shows a global displacement, in the forming direction, of the pyramidal shape. This displacement of 0.06 mm average can be considered as a geometry defect. However it can provide an advantage: a switch in micro-electronics could be designed by using the geometrical defect as a constrained displacement for the equivalent spring.

Qualification of the numerical final shape
One of the most important criteria to validate the formed part is associated to the validation of the overall geometry. This inspection n is realized by comparing the numerical and experimental shapes. The real part, presented in Fig. 4a was digitized by non-contact 3D laser scanning system (Steintek Mobilescan 3D) in high resolution mode. This optical method provides a high density of measured points (more than 100,000) with a resolution of less than 0.01 mm and an accuracy of 5 lm. The comparison with the numerical shape requires the extraction of the outer surface mesh. From the developed postprocessing toolbox, this operation is realized automatically by extracting and converting the outer surface quad mesh in triangular mesh without loss of geometrical definition. The representation of the different meshes (solid mesh, extracted mesh and subdivided mesh) is done in Fig. 9.
This operation is necessary to create a compatible mesh (STL format) for the inspection software (Geomagic Qualify 2012).
The experimental data (measured points) are then compared with this triangular mesh. The initial operation consists in fitting the two geometries in two steps. Firstly, a manual positioning of the experimental scatter with the triangular mesh is done thanks to the advantages of the pyramidal shape: theoretically this shape possesses two planes of symmetry that permit an easier positioning of the two parts. This is more difficult in the case of axisymmetric parts. Secondly, a precise positioning is ensured by the way of a best fit method. Once these operations are performed, the geometric comparisons can be conducted. The different inspection results are shown in Fig. 10.
The 3D comparison presented in Fig. 10a, shows an excellent matching of the experimental and virtual parts (an average difference of -0.031 mm for a standard deviation of 0.056 mm) excepted at the edges of the workpiece (±0.2 mm) and at the base of the pyramid (0.137 mm). These differences can be easily explained: -the experimental part is quite flexible due to its thickness and so it is easy to deform the part away from the forming zone during handling, -the value of the die radius and its position strongly influence the final shape, -the actual die radius is not strictly identical to that used for the simulations.
In addition this area is highly dependent on springback effect which is again a difficult parameter to control. The second comparison (Fig. 10b) is based on the local profiles of the pyramid obtained by cutting the 3D geometries with different planes parallel to a symmetry plane of the part (plan no. 7). All the cutting planes are separated by a distance of 0.4 mm. All the resulting profiles are plotted in their respective cutting plane in Fig. 10d: profiles in black are associated to the experimental measurements as the gray ones with circle markers correspond to the numerical results). These profiles are also shown in Fig 10c. These results confirm the excellent correlation between geometries obtained by experiments and simulations, except for the base of the pyramid. No pillow effect, as suggested by Ambrogio et al. [24], is observed for the small dimensions considered in the present study.
Finally, some observations are revealing the influence of the helical strategy. The numerical results presented in Fig. 11 show a mesh which is twisted about the symmetry axis parallel to the z-axis and denoted z c .
A significant rotation of the mesh is observed in the forming area as presented in Fig. 11a where one half of the mesh after forming simulation is represented. A twisting angle parameter denoted c is introduced to quantify this distortion. This angle can be evaluated by considering a node A in the mesh in its original and deformed position. We can construct the isovalues of the torsion angle shown in Fig. 11b, where its maximum value is 25.7°. The places of maximum twist angle are in areas where the tool undergoes a change in direction, and this area increases as the forming zone decreases.

Thickness evolution
To complete the validation of the model, local comparisons with experiments were carried out. The most important disadvantage of the incremental sheet forming process is an important thinning of the sheet. In this section, the thickness     is evaluated and locally compared with the experiments. The thickness is calculated by measuring the normal distance between the inner and outer surface of the sheet on the whole part. The thickness repartition is plotted in Fig. 12b in a projection view defined by Fig. 12a.
The thickness repartition is very similar to that of the twist angle presented in Fig. 11b. This result shows that the twisting phenomenon, and implicitly the forming strategy, influences the thickness distribution. The minimum value of thickness is 0.07 mm and corresponds to a thinning of 66%.
For validation, the experimental part has been cut in two symmetric parts according to the section plane defined in Fig. 12a. The cut was performed by wire electro-discharge machining process (WEDM) to avoid introducing any additional mechanical stresses in the specimen. Numerical and experimental results are given in Fig. 13.
In Fig. 13a, the numerical model section is compared with the experimental one. The global shape of the numerical simulation is closed to the experiments and thickness evolutions plotted in Fig. 13b confirm these observations. Experimental data are obtained by optical metrology and image processing. The 2D numerical profile of the cutting part is extracted by the way of a pedestrian alpha shape extractor algorithm [25]. The thicknesses evolution is presented in Fig. 13b and shows a good correlation between experiments and simulations. The numerical evolution follows the sine law reported by Reagan and Smith [26], Jeswiet et al. [2] and also discussed by Jackson and Allwood [27].

Forming forces
The last comparison is based on the validation of the mechanical behavior and the friction laws. It is possible to investigate their effects by observing the evolution of the forming forces during the process as proposed by Szekeres et al. [28]. With these considerations, the forming forces in the three directions (x, y and z) are given in Fig. 14.
The results show an excellent prediction of the forming forces. Variations and levels are close to the experiments. In the case of z-axial force, a slight difference is observed after 30% of the cycle time. The experimental force is lower than its numerical counterpart of about 15%. A possible hypothesis is related to the softening phenomenon. The forming process results in large deformation (with plastic strains of more than 243% locally) and softening may be due to damage for example   [29]. As a consequence damage modeling is a way of improvement for the description of material behavior and for predicting risks of cracks. Experimental results with the constant Z-level strategy can be utilized to perform the complete identification of softening model and fracture evolution.
Another assumption may be the influence of the material microstructure (grain size) on the material behavior. The influence of scale effects [30] can be described by the Hall and Petch relation [31,32]. This relation defines the evolution of the yield stress with the grain size. This may explain the observed changes in the evolution of the forming forces. Resulting forming efforts are nonetheless consistent with the experimental results in first approximation.
Monitoring efforts during incremental forming process could allow the identification of material and friction behavior as they directly influence their levels. This can also be a method to identify the occurrence of defects. Mention may be made of the use of SPIF for the identification of Gurson's damage model [33,29] and the definition of failure in order to predict the occurrence of defects during simulations [34][35][36].

Conclusions and perspectives
This research focuses on the development of single point incremental sheet forming process and its numerical modeling. A complete numerical toolbox has been developed to perform fully parametric simulations with finite element software. The toolbox includes a post-processing tool to analyze the results to conduct comparisons with experimental data. These comparisons lead to the following conclusions: 1. Comparisons made on the overall geometry of the workpiece reveal very good agreement between the part obtained by simulation and that obtained during experiments. Only the parts at the base of the base of the pyramid show significant differences. These deviations are sensitive to the position and value of the die radius. Despite these differences, comparisons conducted at the section level are very convincing. 2. The prediction of thickness distribution is close to that obtained on the real part. 3. Forming forces obtained by numerical simulation show good correlation with measured values. However, it was a slight underestimation of the axial forces during thinning. It is assumed the influence of grain size and softening on the material behavior.
The developed approach has demonstrated its ability to efficiently predict the forming of small parts with thin thickness by the ISF process. In terms of perspectives, the presented methodology can be implemented to introduce size effects into the material and friction laws by adjusting either the grain size on the dimensions of the test (specimen and tools). The helical strategy is well suited for identifying the friction and the material behavior laws. As the constant Z-level strategy leads to fracture, it seems appropriate to identify failure models. As a result the SPIF process may be a specific characterization test for material and friction under complex loading.