Selective element fission approach for fast FEM simulation of incremental sheet forming based on dual-mesh system

Incremental sheet forming (ISF) is a highly versatile and flexible process for producing low batches of sheet metal parts. Although finite element (FE) method is a key approach in the study of material deformation in metal forming processes, the application of FE in ISF process is still limited. This is due to the enormous simulation time required for ISF processes. Focusing on this problem, this paper presents a new selective element fission (SEF) approach for simulation of ISF process based on LS-DYNA. In the approach, the computational efficiency is improved by reducing the unnecessary calculations of elements outside the localized deformation zone in the ISF process. The introduction of a background mesh for simulation data storage and a separate simulation mesh with varied mesh density for simulation ensures both the efficiency of computation and accuracy of results. To verify the proposed SEF method, two ISF case problems including a hyperbolic cone part and a pyramid part are studied by comparing to the conventional FE approach and the H-adaptive approach in terms of the CPU time, the forming load, the final part profile, and thickness distribution. The influence of two key factors: element size ratio and toolpath segment length is also studied. The results suggested that the developed SEF approach can save the CPU time by up to 74 % with satisfactory accuracy as compared to conventional FE method, which demonstrates both the effectiveness and robustness of the developed SEF approach.


Introduction
Incremental sheet forming (ISF) is a flexible sheet forming method, which is commonly used for rapid prototyping of sheet metal parts in small batch manufacturing and for customized products. During the ISF process, a ball-nose forming tool is driven to move along a pre-defined toolpath according to the part design. Along the tool movement, the sheet is gradually deformed to the designed shape without the use of a forming die. During the ISF process, localized deformation of sheet occurs around the area where the tool is in contact with the blank. Due to the incremental deformation nature, ISF process has obvious advantages over conventional forming processes including greater formability, reduced forming load and lower cost of tooling and forming equipment [1]. These advantages have attracted increased attention in the metal forming research community, and this process was considered to have the potential to revolutionize sheet metal forming industry [2]. In recent years, progress has been made in many aspects of ISF technology, including in-depth understanding in deformation mechanism [3], development of new toolpath generation strategies [4], improved the geometric accuracy [5], new method for thickness prediction [6], and new material processing methods [7]. These advances have helped to make the ISF technology closer for practical applications in different industrial sectors.
In the study of the metal forming process, finite element (FE) method is a key approach in understanding the material deformation of complex processes [8], forming load using different tool paths [9], effects of process uncertainty [10], and dimensional accuracy of forming parts [11]. This is also true for using FE method to simulate conventional sheet metal forming and ISF processes. Concerning the ISF simulation research, Qin et al. [12] presented a numerical investigation on the influence of forming process parameters by modeling the forming process. Eyckens et al. [13] applied finite element method using sub-modeling technique to model the small plastic zone in the incremental forming process with an adequate fine mesh at a computationally acceptable cost. Ayed et al. [14] proposed "simplified analysis modeling" for numerical simulation of incremental sheet forming process. By considering the constitutive strains/ stresses equations and the tool-sheet contact conditions, this model was developed by using shell element with satisfactory results. Henrard et al. [15] studied the influence of material law, FE code, and yield criterion on the accuracy of ISF FE simulation especially in predicting the tool force. These existing works indicate a good potential of using FE approach to analyze the ISF process.
Although FE is an effective mean in ISF research, it also has some considerable limitations for ISF processes as compared to other metal forming processes, such as forging and stamping. Excessive simulation time required for ISF processes is a significant limitation. For example, Bambach et al. [16] reported that 73 CPU hours was required to implement an FE approach to simulate a simple case of forming a cone part. In a similar case, Henrard et al. [17] also pointed out that the simulation of the whole cone was time-consuming even with a very coarse mesh. Two key reasons of the long computing time can be considered as follows [18]: (i) The sheet blank is deformed by tool along a sequence of small increments, which cause hundreds of thousands of increments to be calculated when toolpath is long; (ii) As the contact area between tool and sheet are very small, very fine mesh is required in the contact area to ensure accurate contact and deformations calculation, which requires a large number of elements to be involved in the simulation. Due to the expensive computation time according to the above factors, the application of FE simulation in ISF process is still limited.
Concerning the reduction of computational scale, many efforts have been made. Muresan et al. [19] suggested a decoupling method to study the deformation of localized area. Sebastiani et al. [20] presented a decoupling algorithm to decrease the computing time in ISF simulation. This algorithm focuses on separating the FE model into an elastic zone and an elastic-plastic deformation zone. During the simulation, timeconsuming iterations focus on the plastic zone and only a few key elements are involved. Hadoush et al. [18] proposed a direct sub-structuring method to reduce the computation time. In their work, the assumption was made that plastic deformation zone exists only in local contact area. In this way, the substructures can also be divided into two categories: the plastic nonlinear substructures and the elastic pseudo-linear substructures. Bambach et al. [21] developed a method that introduces a fine-size mesh template which moves along the forming tool instead of adaptive remeshing method to speed up explicit FE simulations. Using this approach, a reduction of up to 80 % in CPU time can be achieved under specific conditions of the simulation cases.
Another approach for reducing the computational time is based on the simplification of contact model. Robert et al. [22] presented a simplified model for the tool-sheet contact conditions for SPIF process simulation. In their approach, large incremental steps can be used, and approximately 65 % of CPU time can be saved. Brunssen et al. [23] introduced a primal-dual active set formulation for material involving elastic-plastic behaviors. Comparing to the conventional penalty function contact algorithm, primal-dual active method allows for less calculus problem, and possible deterioration of stiffness matrix in penalty function method can be avoided. Henrard et al. [24] develop a new model dealing with the tool-sheet contact. This method enables moving contact modeling anywhere on the element surface without the use of penalty method. This method could save considerable simulation time for the ISF process. In spite of the above efforts made, computing time for FE simulation of ISF processes is still considerably long as the time step restricts the calculation efficiency due to the continuous changing of contact conditions [25].
In this work, a selective element fission method was developed to improve the computational efficiency in FE simulation of ISF processes. In the proposed approach, the computational efficiency can be increased by reducing the number of elements that involves in the FE model. This is achieved by assigning different mesh density on the blank sheet according to pre-defined split toolpath segments and the prediction of potential deformation regions. In this way, the total number of elements involved in the calculation can be reduced without losing the calculation accuracy especially at the plastic deformation area. In the following section, the methodology of the selective mesh fission approach is described. To validate the efficiency and accuracy of the developed approach, two cases of ISF processes are studied by comparing to the conventional FE method with the developed selective element fission (SEF) approach in terms of computing time and accuracy together with detailed evaluation of forming load and thickness distribution. This is followed by discussion and conclusions that may be drawn of the developed method.

Methodology
In the conventional ISF simulation process, elements with smaller size have to be employed as the area of contact between tool and sheet is small. These small-sized elements require a large number of elements and significantly increase the simulation time. To overcome this problem, commercial software such as LS-DYNA has developed strategies such as H-adaptive mesh, in which the mesh density increases gradually at the critical area of tool-sheet contact [26]. However, in this H-adaptive mesh approach, the calculation becomes slower with the increase of element number due to movement of contact zone: the mesh density at the previous contact zone cannot be reduced when the contact area is moved forward. The accumulation of high-density mesh increases the intensity of computation. However, if the mesh density at the previous contact zone was reduced, which requires the replacement of the previous smaller elements by larger-sized elements and the accuracy of previous calculation will be lost. To overcome this problem, a new strategy is proposed by employing a predefined adaptive mesh strategy and an approach of data storage on a denser background mesh.

General description of the selective element fission approach
During the incremental sheet forming process, plastic deformation takes place locally at limited area [27]. This area normally covers the tool-sheet contact region and its vicinity. In the simulation process, the use of small elements that are far away from those areas would not affect on the computation accuracy; instead, it significantly reduce the computation efficiency. To overcome this problem, different mesh densities may be assigned to different regions: in the deformation zone, finer mesh may be employed to increase the computational accuracy; in the non-deformation zone, coarser mesh may be used to improve the computational efficiency. Due to the localized deformation in ISF processes, the deformation zone always changes with the tool movement. However, the accumulated total deformation zones for ISF process cover most of the sheet area. If the denser mesh is assigned to the whole deformation area, it would generate large number of elements. Instead of simulating the whole ISF toolpath at a time, the toolpath is divided into a number of segments and the FE simulation is implemented for each toolpath segment in sequence. In this way, denser mesh is only required for a particularly small deformation area in that toolpath segment other than the whole area so the number of elements can be reduced.
To facilitate the proposed strategy for ISF simulation, the ISF toolpath is split into a number of segments, and each toolpath segment will be simulated in sequence. For the simulation of each toolpath segment, the coarser mesh is employed and some of the elements in the coarser mesh are divided into denser meshes if it is established that the tool passes through a specific toolpath segment. With this varied mesh density, the element number can be reduced comparing to a full denser mesh. After each FE simulation, the result data are transferred into a fine background mesh for storage. Concerning the simulation of another segment, based on the updated geometrical shape and state variables from background mesh, new simulation mesh can be generated and FE simulations can be carried on segment by segment until finish. Figure 1 described this simulation principle of the proposed SEF approach. The key steps of this SEF approach include (1) toolpath division and simulation planning, (2) selective element fission, (3) data transaction and interpolation, and (4) generation of new mesh. It is noted that the number of element required in the actual simulation mainly depends on the length of a toolpath segment (or the number of segment division) as well as the size ratio. This size ratio (SR) may be defined as the area of an element before fission comparing to the area of element after fission while the segment number (SN) suggested the amount of segments or in other words the number of calculation steps. Thus the element SR and the SN are the key parameters that affect the simulation accuracy and efficiency.

Toolpath division and simulation planning
Concerning the toolpath generation in the ISF process, the designed part may be sliced in the vertical direction thus a series of contours can be obtained. Based on these contours, spiral toolpath may be generated by interpolating the contours or these contours may be directly employed as layered toolpath. Other types of toolpath may also be employed for ISF process, such as the feature-based toolpaths generated from the equal-potential lines of the designed part or the Zig-Zag type toolpaths [4].
No matter what specific toolpaths are used, these paths are generally formed by a series of discrete points on the trajectory of tool movement. In this work, the ISF toolpath was divided into a series to toolpath segments with equal length. Each segment is simulated in the order to form an integrated ISF simulation. In this way, the computational efficiency can be improved.
Concerning the toolpath division, if there are too many segments, more data transformation operations between the simulation mesh and the background mesh for data storage are needed. The large number of data transformation operations may potentially reduce the simulation accuracy due to the loss of data accuracy in the interpolation process. On the other hand, if the toolpath segment is too long, the potential deformation area is increased which requires a finer mesh with a larger number of elements. Because of this, a proper segment length plays a key role in balancing the computation efficiency and simulation accuracy. With a given segment number SN, the length of toolpath segment L S can be determined by: In the toolpath division, as the generated toolpath is represented by a series of points, the division of toolpath segment can be obtained by: where k is the start point number of a toolpath segment while n is the number of points in a toolpath segment. In practice, with a given target toolpath segment length L S and a starting point k, the number of points n that involves in this toolpath segment can be calculated by summing up the length of toolpath points until the accumulated length is approximately equal to the target length L S .

Selective element fission algorithm
As mentioned in the previous section, the deformation of sheet metal occurs in a localized tool-sheet contact region. With properly defined toolpath segment, the potential deformation area during the simulation process can be predicted. In this way, fine mesh can be assigned to the potential deformation area while coarser mesh can be assigned to the other area.
In ISF simulation, structured mesh with shell elements is usually employed. For an initial large shell element with square shape, this element can be divided into smaller square-shape elements by equally splitting the element edge into two or more edges and adding new nodes for the newly generated elements. Figure 2 shows the fission of coarser element into smaller element. It is worthwhile to point out that the element fission method described above would result to an incompatibility at the element edge with different element sizes. The middle node on a neighboring edge shared by elements with different sizes is calculated by interpolating from the two vertex nodes on the neighboring edge [28]. In this work, this calculation is managed by the LS-DYNA during computation, which ensures the mesh deformation compatibility [26]. Although the use of unstructured mesh could also produce FE mesh with different density, this kind of mesh may cause complex data transformation operation at later stage. Thus, structured mesh seems a better choice for the SEF solution.
Another issue in the element fission method is to determine which elements are needed to be split. In this work, this is achieved by looking for any elements that are close to the toolpath segment. As presented in Eq. 3, for any point P j on the toolpath, if the distance d N-P of an element node N i to the point is within a critical distance d 0 , this element is considered to be in deformation affected zone and is thus be split.
The critical distance d 0 directly determine the number of elements to be split and affect the simulation efficiency. An appropriate value of d 0 would enhance the simulation efficiency and at the same time maintain satisfied simulation accuracy. In the incremental forming process, the critical plastic deformation zone is limited to the tool-sheet contact and the surrounding area [27]. The value of d 0 has to be large enough to include this area so the plastic zone can be covered by the refined mesh in the simulation. In this study, d 0 is twice the radius of the forming tool. Using Eq. 3, all the nodes near the toolpath can be found and the elements that include these nodes are identified and divided into smaller elements. In this way, meshes with different density can be generated for a particular toolpath segment. This process is illustrated in Fig. 3.

Data transformation and interpolation
To store the simulation data, a background mesh is employed in the approach. To avoid the loss of accuracy during data transformation, the element size of background mesh is defined as the same size of the smaller element in the simulation mesh. In doing so, the nodes on the background mesh correspond to those on the simulation mesh. The nodal value such  as updated coordinates and thickness can be directly assigned to the corresponding nodes on the background mesh. Similarly, elemental values such as stress and strain components can also be assigned to the corresponding element in the background mesh as well. For those nodes and elements that do not have direct correspondences due to the different element size between the background and simulation meshes, data interpolation method are employed by using the shape function approach: where N i represents the shape function of node i and X i are the state variables to be transformed. In this way, after the completion of FE simulation for each toolpath segment, the background mesh is updated according to the newly calculated results from simulation mesh as shown in Fig. 4.

Generation of new mesh
Before simulation of each toolpath segment, a new simulation mesh has to be generated according to the current state of background mesh as well as the toolpath to be simulated. Based on the geometrical shape of present background mesh, new mesh is formed and finer mesh is generated along the toolpath according to the SEF approach proposed in Section 2.3. After the formation of meshes with different density, state variable results are transferred from the background mesh to the newly generated simulation mesh. Similar to the data transformation technique described in Section 2.4, if the elements are kept the same, the state variable results are directly assigned to the corresponding elements. If the element in the simulation mesh is larger than the corresponding element in the background mesh, the average values is taken from the corresponding background elements as presented in Eq. 4. In this way, a new simulation mesh can be generated as shown in Fig. 5.
where X i are the state variables to be transformed and n is the number of elements in the background mesh which correspond to a larger element in the simulation mesh.  In this work, the implementation of SEF approach is achieved via an in-house developed C# program and the commercial software LS-DYNA. The framework of this program is illustrated in Fig. 6. As shown in the figure, for any given toolpaths, they are divided into several segments according to a target segment length. For each toolpath segment, ISF FE simulation is implemented in sequence: for a given toolpath, elements of a coarse mesh are split into smaller ones according to the position of the toolpath and ISF simulation is carried out on the refined mesh. After that, the simulated result data including the geometrical information, stress, and strain results are updated corresponding to the background mesh. Based on the stored data in the background mesh and new toolpath information, a new FE input file can be generated and submitted to LS-DYNA solver for simulation. The newly simulated results of the background mesh are updated again. Therefore, the whole toolpath can be simulated on a segment basis.

Case studies
In order to evaluate the accuracy and efficiency of the developed approach for ISF simulation, two case study problems including a hyperbolic cone and a pyramid are performed for verification purposes. In each case, the FE simulation results are compared by using conventional FE and the H-adaptive method [26] within LS-DYNA as well as the developed SEF approach. The computing time is also measured and compared so the efficiency of the developed approach can be evaluated.

Experimental setup
The experiments of the two cases were performed on a 3-axial CNC machine as show in Fig. 7. In the ISF process, a ballnose tool with radius of 5 mm was employed. The feed rate of the ISF process for both cases is 1200 mm/min. Spiral toolpath was generated for the ISF process. In case 1 of the hyperbolic cone, a constant step down value (the movement of tool in the vertical direction in each toolpath contour [27]) of 0.5 mm was employed while for case 2 of the pyramid, step down value was defined to be 0.3 mm. Aluminum alloy AA5052 with an initial thickness of 1 mm was used for validation. The stressstrain curve from the uniaxial tensile test for the AA5052 sheet was illustrated in Fig. 8 and the material properties were calibrated as presented in Table 1. Using these data, the power law material model in LS-DYNA was employed for simulation. The initial size of the sheet is 180×180 mm. In order to reduce the friction between tool and sheet, lubricating grease was employed between the tool and the blank. During the ISF process, a backing plate was placed under the blank sheet to ensure geometric accuracy. Clamps were applied to the edge of the blank. In this way, the blank can be fixed and inaccurate  bending of the sheet can be avoided. During ISF experiment, a JR3 multi axis load cell was employed to measure the forming force. After the forming process, a laser displacement sensor was used to measure the cross-section shape of the finished part.

Case 1: hyperbolic cone
A hyperbolic cone shape, as a commonly used geometry for ISF process, is employed for the process evaluation as shown in Fig. 9a. The parts obtained from FE simulation and experimental testing is presented in Fig. 9b, c, respectively. In order to compare the computational efficiency using different methods, a conventional FE simulation and the newly developed approach are compared. For the newly developed SEF approach, the elements with different size ratio from 4 to 16 are employed. In all these cases, the size of initial small element was kept constant as 1 mm while the size of large element varies according to the size ratio. In addition, different toolpath segments are also compared, in which the toolpath segment number are set as 474, 242, 123, and 62. In all the simulations, the friction coefficient was set to be 0.05 for the contact between tool and sheet. The virtual speed is scaled up by ten times according to the best practice reference, in which the ratio of kinematic energy to total energy can be controlled within a limited value.
In order to evaluate the computational efficiency, the developed fast simulation approach under different element size ratio (SR) and toolpath segment number (SN) are compared with the conventional FE approach as well as the H-adaptive approach. In the H-adaptive simulation model, the refinement level is set to 3 and the original element size is 4 mm, i.e., the same as the model SR=4. As can be seen from Fig. 10, a maximum reduction to 26 % of original simulation time can be achieved from a size ratio of 16 and segment number of 474. As can be seen from the figure, by employing a larger size ratio and larger segment number, the computational efficiency can be improved. As given in Table 2, it can be found that larger size ratio and larger segment number would reduce the element number in the simulation.
Concerning the computational accuracy, forming load in the vertical direction is compared between simulation and   Fig. 11. As can be seen from the figure, all the forming loads from simulations are slightly higher than the load measured from experiments. The forming load predicted by the H-adaptive method is always higher than the conventional FE approach. Concerning the effect from the size ratio and segment number, it can be found that the prediction results become slightly worse when the size ratio increases. By increasing the segment number, the forming load becomes lower and closer to the experiment result.
For the geometric accuracy, the final part profiles are compared as shown in Fig. 12. It can be found that the predicted results from all cases of simulation are almost identical. However, a maximum discrepancy of about 0.7 mm can be observed between the simulated profile and the actual profile. This result suggests a good agreement between numerical results and experimental testing. It also indicates that the size ratio and segment number have limited impact on the final part profile.
The thickness distributions are shown in Fig. 13. As can be seen in the figure, the thickness distribution obtained from the conventional FE simulation is quite smooth. However, for other approaches including the H-adaptive approach and the developed SEF approach, although the minimum thickness reduction value is similar, the predicted distribution results are quite bumpy. This bumpy curve may be caused by the interpolation of the state variable data during the data  transformation process, which may be overcome by introducing improved data interpolation algorithms in the future. By examining this cone shape case problem, it can be found that simulation time can be reduced to about 26 % of its original time under the size ratio of 16 and the segment number of 474. Concerning the simulation accuracy, satisfactory results can be obtained especially for the final part profile and the thickness distribution. In this way, the developed approach may be employed for quick simulation in forming load and thickness reduction evaluation.

Case 2: pyramid
In this case problem, the dimensions of this pyramid are shown in Fig. 14a while the simulation and the experimentally made parts are illustrated in Fig. 14b, c. Similar to the previous case of the hyperbolic cone, both FE simulations and experiment were implemented so the efficiency and the accuracy of different simulation methods can be evaluated. In FE simulations, the conventional FE approach, H-adaptive approach, as well as the developed SEF approach with varying size ratios and segment numbers were employed. In all these cases, the size of small element was kept constant as 1 mm while the size of large element varies according to the size ratio.
The same as the previous case, conventional FE method example and three remesh method examples are conducted. The H-adaptive method example is set according to the remesh example in which size ratio is 16. Figure 15 shows different computing times required under specific conditions. Compared to the previous case, a more than 70 % reduction is also obtained by using the SEF approach. Note that the Hadaptive method example requires less simulation time than the corresponding remesh method example when SL = 200 mm. This discrepancy may be resulted from two sources: (1) compared to the previous case, the forming area is smaller. Element reduction is not as obvious as the previous case. (2) In case 1, the forming step is 0.5 mm while it is 0.3 mm in case 2. In H-adaptive method, once an element is divided, the small elements continue to be part of the model until the end of the simulation. But in the remesh method, small elements are changed back to large elements and in the next step, new small elements are generated at the same position and the results of nodes and element are transferred for a second time   Fig. 16 Forming load curve and hence more time is needed than the adaptive method. When smaller SL value is used, the SEF approach understandably is more computationally efficient than the H-adaptive method, as shown in Table 3. Compared to the conventional and H-adaptive methods, the same length toolpath requires less time. Figure 16 compares the vertical forming load history from the experiment and FE simulation results. As can be seen from the figure, the results from conventional FE approach is quite close to that obtained from ISF experiment. Concerning the Hadaptive method, large fluctuations of the forming load can be observed. For the other cases, the prediction values are more stable. For the developed approach, as the size ratio increases, the forming load becomes less accurate and the error reaches 20 % in the case of size ratio 16. By increasing the segment number, the forming load is slightly lower, similar to the results obtained in case 1. In consideration of the computational efficiency gains achieved from the proposed SEF approach, the prediction results are still satisfactory when a reasonable size ratio is chosen.
To examine the geometrical accuracy, the scanned part profile is compared with the numerical results as shown in Fig. 17. Similar to the previous case of the hyperbolic cone part, the predicted results from FE simulations are almost identical. The maximum geometrical deviation of the formed pyramid is 0.7 mm for the formed pyramid shape obtained from the ISF experimental testing and FE simulations. This result suggests a satisfactory accuracy of the final part geometry. In addition, it is also noted that the size ratio and segment number has limited effect on the geometrical prediction.
Thickness distribution is another measurement to evaluate the simulation accuracy of this case study. Figure 18 shows the comparison of the thickness distributions of this pyramid part from both numerical and experimental results. As can be seen from the figure, the trends of thickness variation are quite similar. The FE simulation results are quite close to each other, and these results are similar to the actual thickness distribution from experiment. Again, this result suggests satisfactory numerical prediction and the effect from size ratio and segment number may be ignored.
By examining this pyramid case problem, it can be found that simulation time can be reduced to about 26 % of its original time at a size ratio of 16 and segment number of 808. Satisfactory prediction can be obtained particularly in geometrical profile and thickness distribution. This case study again proves the robustness of the developed simulation algorithm.

Discussion
The developed fast simulation method is based on selective element fission approach, which enables the reduction of the number of elements involved in simulation by representing the localized deformation nature in the ISF process. The case studies of a hyperbolic cone and a pyramid shape suggest considerable simulation time reduction and hence demonstrate its feasibility and robustness of using the developed approach in ISF processes. Comparing to the conventional FE method, the developed approach consumes much less computing time but at the same time maintains satisfactory accuracy in terms  Comparing to the other approaches, it is relatively easy to implement as only limited programming on element fission and data transformation is necessary without requiring complex FE formulation and derivation, which has a good potential of widespread applications in ISF process. In addition, the principal idea developed in this approach including the split of tool paths for individual simulation and the employment of dual-mesh system for data storage may also be introduced to other incremental-typed metal forming processes with localized deformation such as metal spinning or rolling.
In the proposed method, the size ratio and the number of toolpath segments are two key parameters that affect the computational efficiency. The study of two case problems in this work suggests that large size ratio with increased segment number tend to reduce the simulation time. However, the study on forming load suggested that larger size ratio may produce less accurate results. This inaccuracy may be attributed to the unstable element thickness: the variation of sheet thickness may affect the plastic work and hence the load required for deformation. As can be observed in the two cases, although the trend of thickness distribution is quite similar to that obtained from experiment, the distribution from the fast simulation approach is more uneven than that from conventional FE simulation. In addition, similar observation of fluctuation may be also found for the load prediction. This may indicate the unsmooth contact condition or unsmooth data interpolations values, which may be due to the loss of accuracy during data interpolation process especially for shell thickness and part geometry. Effort is to be made in the future for further improvement of the thickness distribution and load prediction by employing more efficient data interpolation, surface approximation, and thickness calculation methods. In addition, this SEF approach may also be applied in ISF simulation with solid elements for future research.

Conclusions
In this study, an FE-based selective element fission (SEF) approach for incremental sheet forming simulation process is presented. By using this method, the FE computing time, geometric profiles and the thickness distribution of form parts, and the vertical forming load are compared to that of the conventional FE method. To validate the proposed method, two common ISF geometrical parts including a hyperbolic cone shape and a pyramid shape are modeled by using the developed SEF approach. The results of two case studies showed a good agreement between the developed SEF approach and the experimental study. Conclusions of this work may be summarized as follows: 1. By considering the local deformation nature of ISF process, selective element fission is an efficient approach for FE simulation of incremental sheet forming process. 2. Satisfactory simulation accuracy has been observed comparing to existing H-adaptive method, especially for thickness distribution and the final part geometry. 3. The two key parameters of element size ratio and toolpath segment number, considered in this work, have a notable influence on the computational efficiency as well as prediction of forming load; however, they have less effect on the geometric accuracy and thickness distribution.