Parameters for the FEA simulations of single point incremental forming

ABSTRACT Single point incremental forming has shown immense potential in the manufacture of prototypes. The sheet metal deformation in this technique is carried out by a small hemispherical or flat tool that moves from the sheet periphery to the sheet center while also pushing the sheet down. Simulations of single point incremental forming can be a tedious task since the area required to be formed is large and the two sheet contact area is fairly small. The effect of tool however, affects the entire geometry by causing deflections and radial strains. Hence, it becomes essential to understand the effect of certain parameters that influence the simulation results for this forming process. The paper aims to study factors influencing the simulation of single point incremental forming process. It gives a description of techniques employed to simulate non-linear behavior of sheet metal forming along with the implicit and explicit methodology. Information on the selection of element type including Reduced Enhanced Solid-Shell element is presented. Techniques employed to avoid locking are also discussed.


Introduction
Deep drawing has been used over the years to form metal sheets by using dies and punches to obtain required geometrical shape. This proved to be a cost effective process for mass production of components. However, the requirement of a special tooling/die for every shape makes an expensive and time-consuming process (Kopac & Kampus, 2005). Recent changes in the consumer behavior have led to an increase in the demand for customized components (Farias et al., 2014). This has resulted in high manufacturing cost for the small batch production especially when the required batch sizes are fairly small (Hirt, Ames, Bambach, & Kopp, 2004).
Several new processes and techniques have been developed and proposed overtime for flexible and cost effective manufacturing in sheet metal forming industry (Fratini, Ambrogio, Di Lorenzo, Filice, & Micari, 2004). Asymmetric Single Point Incremental Forming (SPIF) has gained interest over past few years with continuous studies being carried for further optimization (Jeswiet et al., 2008). Several studies are being carried in both numerical and experimental field to understand SPIF to optimize the process for better industrial implementation (Farias et al., 2014). With the increase in possible FEA have also been used to understand the complex forming phenomenon that is an important facet of SPIF. A geometrical model was developed to predict the membrane strains and sheet thickness to provide better approximations than the sine law. It was validated by comparing the results obtained by sine law, experiment and the FEA (Bambach, 2010). Implicit approach was used to simulate the sheet discretized with three layers of brick elements to explain the twist phenomenon by strain measurements (Duflou et al., 2010). The effects of various process strategies were investigated to improve the process precision. It aimed toward reducing sheet bending near the backing plate and springback. The process was simulated by discretizing the sheet with two layers of linear solid elements and using implicit method to solve the equations. FEA have been utilized to study the effect of various process strategies to improve the process precision (Essa & Hartley, 2010). In another study, the experimental results were compared to the numerical simulation results obtained by explicit method to investigate the sheet thinning. FEA were successfully implemented to analyze stress and strain during the sheet deformation (Shanmuganatan & Senthil Kumar, 2012). Various studies have investigated the parameters like wall angle, tool diameter, incremental step size and material thickness critical in SPIF by the aid of numerical models (Arfa, Bahloul, & BelHadjSalah, 2012;Bahloul, Arfa, & Belhadjsalah, 2013;Shanmuganatan & Senthil Kumar, 2013). FEA are now widely employed for better understanding of the defects occurring in the formed component. Numerical model was used to examine defects exhibited in SPIF namely corner fold, pillow and wall by using shell elements (Hussain, Al-Ghamdi, Khalatbari, Iqbal, & Hashemipour, 2014). Another study examined the pillow effect by the aid of experimental and numerical simulation results .
FEA have shown promising results in their ability to predict the SPIF outcomes in close proximity with the experimental results. However, successful prediction depends upon the type of calculations selected, element type and computationally efficient strategies employed. Understanding of these parameters is vital for proper execution of the finite element simulation. Some of the most important parameters for the SPIF process will be discussed in this paper.

Iterative methodology
High material and geometric nonlinearity is observed for incremental forming process. The software and FEA technique employed for the simulation should be capable of taking nonlinearity into account. This includes consideration of material properties exhibiting accurate predictions of strain, shape and thickness distribution. This nonlinearity is a result of material yield, creep, large deflections or shape changes within the system of applied forces. The flexibility of the system changes as deformation takes place in the system (Whirley & Engelmann, 1993). In order to accommodate nonlinearity, change in stiffness matrix of the system is required (Getting Started with Abaqus: Interactive Edition, 2008).
FEA offers a variety of procedures to solve dynamic problems. Direct integration technique includes implicit and explicit methods to calculate the linear and nonlinear problems. Application of incremental loads that lead to geometrical and material changes may also cause material yielding in the region away from the forming region. This transient excitation in the structural systems is solved by selecting a solver algorithm from direct integration technique (Liu & Quek, 2003). The final outcome is established on the results obtained after every time step for better accuracy. The accuracy also depends on the selection of the suitable algorithm for a specific problem. The selection of the suitable technique depends on factors such as loading, type of material and mesh selected for the problem to be analyzed (Yeung & Welch, 1977). The general matrix equation for the structural dynamics including the inertia forces and damping forces is written as (Cook, 1995): where M is the mass matrix, C is the viscous damping matrix, K is the global stiffness matrix, R is the vector of equivalent nodal forces applied externally, D is the displacement, D˙is the vector of velocity components and D¨is the acceleration. The damping coefficients of C, cK and cM are derived from the experiments and can be further written as follows:

Implicit method
Implicit method is numerically stable and usually more efficient for a slow phenomenon. It calculates the quantities in the current step using the quantities obtained in the previous time step. At each time step, the matrix system is solved to obtain the solution and the process continues till the final time step hence, making the implicit method a time consuming process. The stiffness matrix K varies with time in the dynamic problems and therefore is computationally determined and inversed after being assembled. The constructed global stiffness matrix is extremely large for three-dimensional geometries and storing the matrix becomes a demanding task. K is inverted after NB 2 computations where, N is the number of degrees of freedom and B is the bandwidth of the stiffness matrix. The inverted value of matrix K is further utilized to calculate D. Implicit method makes automatic time increments during each process for calculations and is not governed by a stability criteria. This method is based on the accuracy and convergence considerations. The time increments are determined by monitoring the convergence of solution obtained for a time step. The step is repeated again by reducing the time step if non-converging solutions are obtained. Difficulties are usually associated with the solution convergence in problems with large deformations, large nonlinearity and multiple surfaces in contact. It is possible to increase the process speed (even a few steps can be computationally intensive) by using large time steps when the external activity in the model is of a slow time variation. This process has a high cost per step as a set of nonlinear algebraic equations must be solved at each step and the storage requirements increase with the meshed area. Newton-Raphson method and Newmarks method are generally employed as equation solvers (Harewood & McHugh, 2007;Mashayekhi, n.d.;Whirley & Engelmann, 1993;Yeung & Welch 1977).

Explicit method
Explicit methods are able to predict large deformations in three dimensional contact constraints. They have comparatively easy implementation for simulation than the implicit method. They can advance using known variables from the previous steps by avoiding the use of inverse of stiffness matrix. Use of diagonal element mass matrix M enables efficient inversion of the matrix to compute acceleration D¨. Finite difference equations are employed to calculate the values for next time step until the desired time step. Mean velocities are generally defined with a zero value or any specific value to begin the process. The time marching in explicit technique is extremely fast and suitable for processes with large deformations and highly nonlinear events.
The explicit methods however, become unstable and diverge rapidly over large time increments. This rapid divergence is a result of nonavailability of convergence checks. Hence, the mesh size and the time step in this case are 100-1000 times smaller than that of the implicit method. The smallest time taken by a wave to propagate across the smallest element in the mesh is selected as the critical time step size to get convergence (Harewood & McHugh, 2007); Sun, Lee, & Lee 2000;Theory Manual, 2011;Whirley & Engelmann, 1993;Yeung & Welch, 1977).
A comparison between the explicit and the implicit method has exhibited that the computation time required for simulation of SPIF process by implicit method is double the time taken for simulation by using explicit method even for small parts. However, prediction of accurate geometry remains a problem in the simulations done by the explicit method (Bambach, 2014;Bambach, Cannamela, Azaouzi, Hirt, and Batoz (2007)).

Elements type
Increased formability is still not clear in SPIF and various studies are underway to explain this phenomenon (Emmens & van Den Boogaard, 2009). Two dimensional or three dimensional elements can be utilized to analyze deformation in the sheet caused by various mechanisms. Hence, the elements that are employed in the simulation of plastic deformation of the sheet are membrane, shells and solid elements (Makinouchi & Kawka, 1995).

Membrane elements
The efficiency and less expensive nature of membrane elements give them an advantage over other types of elements. These elements have two translational degrees of freedom and thus are unable to retrieve the bending effects accurately. Springback is not predicted accurately by the membrane elements making their use fairly limited for SPIF (Makinouchi & Kawka, 1995).

Shell elements
Three dimensional geometries can be discretized by shell elements that can model thick sections. Conventional shells developed by the superimposition and coupling of plate and membrane elements. They are able to define the behavior of model in two dimensions and thickness direction (Bathe, 1996;Ellobody, Feng, & Young, 2013;Hoogenboom, 2013;Muthukumar & Kumar, 2015). S4R shell element, from the ABAQUS library is most commonly used shell element for the simulation of the SPIF process. They are based on the thick shell theory and allow user defined integration points through thickness. They are suitable for applications that involve forming problems of large rotations and deformations.
Shell elements provide good results for the prediction of geometry and thickness but are unable to represent full three-dimensional stress field especially when the local stresses under the tool are important. They have a complex formulation and additional efforts are required for rotational degrees of freedom in deformations involving large rotations. Plane stress assumptions in formulation of shell elements require special integration treatment for constitutive equations. Complicated shell theories are required to consider thickness and to make up for the developmental assumptions. In addition, the two-dimensional material law and the assumptions employed in the shells lead to inaccuracies in the description of through thickness variables (Arfa, Bahloul, and Belhadjsalah, 2014;Bahloul et al., 2013;Computational Finite Element Methods in Nanotechnology, 2012;Liu, Daniel, Li, Liu, & Meehan, 2014;Petek & Kuzman, 2010).

Solid elements
The solid elements simulate deformation process where full three-dimensional analysis are required to study the sheet deformation. They calculate the strain variations within the thickness of the sheet, avoid complex shell type kinematics, treat large rotations easily, consider three dimensional constitutive laws and consider contact conditions on both sides of the structure (Hauptmann & Schweizerhof, 1998;Jeswiet et al., 2005;Sena et al., 2011;Trinh, Abed-Meraim, & Combescure, 2011).
The most common and basic three-dimensional solid elements used in simulation are tetrahedral and hexahedral elements. Tetrahedral elements have a constant strain matrix and are also referred to as constant strain elements. These elements are only accurate when the strains are constant over the entire span. These results are poor capability to represent deformations involving bending and twisting (Cook, 1995).
Hexahedral elements are unable to represent the beam bending action accurately as they tend to shear lock during the deformation. However, it has been shown that linear hexahedral elements give much better results than constant strain tetrahedral elements by comparison of results obtained by simulation of a simple cantilever beam (Logan, 2007). The use of these elements to obtain good results can be made possible by use of remedies employed for locking (Cook, 1995;Liu & Quek 2003;Zienkiewicz & Taylor, 2000). A study has employed the linear brick elements with reduced integration as a remedy for locking to analyze the mechanics of fracture in SPIF (Malhotra, Xue, Belytschko, & Cao, 2012).

Solid-shell elements
Solid-shells are intermediates between shells and conventional solid elements and combine the advantages of both solid and shell elements. These elements are based on the shell kinematics assumption but have nodes and degrees of freedom configuration as that of solid elements. They can be more time consuming than the conventional shell elements and similar to solid elements can exhibit locking phenomenon (Duchˆene et al., 2013;Quak, 2007;Sena et al., 2011).
Various developments have taken place in the solid-shell elements since their initial development with only one integration point through the thickness. Locking was found to be eliminated by introduction of plane stress condition or linear strains in the thickness direction. Various ways to prevent locking phenomenon were given in another study (Harnau & Schweizerhof, 2002). An efficient Reduced Enhanced Solid-Shell (RESS) element was formulated that constricted the locking phenomenon and have any number of integration points in the thickness direction. This element formulation combined the reduced integration technique and the EAS for better accuracy. The element was tested for both linear (Alves de  and nonlinear applications (Alves de Sousa et al., 2006). The RESS was employed in another study to analyze the large deformation elastoplastic applications. It was found suitable for simulation of sheet metal forming accurately and efficiently (Alves de Sousa, Yoon, Cardoso, Fontes Valente, & Gracio, 2007). It was further modified to give Modified RESS (M-RESS) by increasing the accuracy and stability of the element. It was done by employing EAS with one enhancing parameter .The RESS element was unable to consider the transverse shear components and could lead to hourglassing in several nonlinear applications. ANS was incorporated to overcome these issues (Cardoso et al., 2008). RESS element along with ANS is found suitable to accurately predict the final geometry. However, the increased flexibility of element has shown to result in overestimation of force (Duchˆene et al., 2013). It was found that selection of four sampling points for ANS results in the maximum reduction of error from the results (Duchene, Ben Bettaieb, & Habraken, 2011).
These formulations have been successful in simulation of complex sheet metal forming processes (Cardoso et al., 2008). RESS element was used in a study to simulate SPIF process. They gave better results than the solid elements with less computational cost (Sena et al., 2011;Sena, de Sousa, & Valente, 2010). A new element SSH3D based on RESS found in good agreement with the experimental results. RESS element successfully demonstrated the capability to predict shape of the cone by making use of adaptive mesh refinement technique (de Sena et al., 2013).

Locking
Locking is a phenomenon observed when elements exhibit overstiff behavior in response to the deformation. It can occur due to various reasons including poor formulation of governing equations. These modes are responsible for reduction in the rate of convergence for a given model.

Membrane locking
Membrane locking is the inability of a shell to represent the extension of the element during pure bending deformations. It occurs due to the inability of shells to bend without stretching. Shells have high membrane stiffness with a low bending stiffness. This abnormality is observed as underestimation of displacements in elements which are required to model bending without stretching. This phenomenon can be eliminated by the reduced integration techniques and assumed strain fields like ANS and EAS (Ambroziak, 2013;Belytschko, Liu, & Moran, 2000;Meek, 1991).

Poisson's thickness locking
Poissons thickness locking is observed in continuum shells in bending dominated problems that assume linear displacement in the thickness direction. However, constant thickness strain is observed in the elements based on this theory instead of linear variation. This effect introduces artificial stiffness in the element. EAS technique is suitable option for the elimination of Poissons thickness locking (Carrera, Brischetto, and Nali, 2011;Chandran, Udaykumar, & Reinhardt, 2010;Schlebusch, Matheas, & Zastrau, 2004;Schlebusch & Zastrau, 2005).

Curvature thickness locking
Curvature thickness locking is observed in shell and membrane elements. This form of locking is caused by transverse bending stress or strain modes. This phenomenon is observed when the vector from the lower edge to the upper node is not perpendicular to the mid surface. The elements tend to lock upon application of bending stress and induce artificial thickness in the element. It can be eliminated by using assumed strain field methods (Echter & Bischoff, 2008;Long, Cen, & Long, 2004;Long, Cen, & Long, 2009;Macneal, 1987;Sze, 2000;Tan & Vu-Quoc, 2003).

Shear locking
Shear locking is one of the most common abnormalities that take place in the shell and volume elements. It occurs largely in the linear elements like eight node brick element that are fully integrated. This phenomenon is observed when the displacement in the element upon the application of forces is less due to exhibition of increased stiffness of the system. Energy in the element is lost in the shearing leaving far less energy to model bending. It can be eliminated by reduced integration procedures and assumed strain fields (Bower, 2011;Dhondt, 2004;Ochoa & Reddy, 1992).

Volumetric locking
This phenomenon occurs when Poissons ratio of the material moves toward a value of 0.5. The displacements in such cases tend to be zero as the material becomes incompressible when the Poissons ratio attains a value of 0.5. Inaccurate approximations of strain field occur as a large amount of energy is lost to the virtual power at any integration point. It can be prevented by using reduced integration procedures and EAS approach for the elements (Ambroziak, 2013;Bower, 2011).

Reduced integration
Shear locking and membrane locking can be avoided by using reduced integration (Arfa et al., 2014;Bahloul et al., 2013;Liu et al., 2014;Malhotra et al., 2012;Mirnia, Mollaei Dariani, Vanhove, & Duflou, 2013;Petek & Kuzman, 2010). It is the easiest way to avoid locking. Under the bending action, full integration scheme is unable to model the bending action of the cantilever accurately. Reduced integration makes use of one order less accurate integration scheme to determine the element stiffness. This technique is able to capture accurate displacement results by neglecting the shear strain in the element (Bower, 2011;Dhondt, 2004;Liu & Trung, 2010;Zienkiewicz, Taylor, & Too, 1971).

Selective reduced integration
This method is suitable for incompressible or materials very close to incompressibility where full integration might cause locking. It is implemented by decomposing the stress tensor in hydrostatic and deviatoric components. The deviatoric part undergoes full quadrature whereas, the integration of term having the hydrostatic component takes place by reduced integration. Consider an example of clamped plate that is deformed through the center. The application of selective reduced integration will use full quadrature of 22 for bending stiffness and reduced integration of 11 for shear stiffness matrix. The advantage of using selective reduced integration is that the rank of the element stiffness and global stiffness matrix will remain unchanged (Bathe, 1996;Belytschko et al., 2000;Gosz, 2005;Hughes, 1987).

Hourglass
Reduced Integration in linear elements can sometimes give rise to zero energy modes that exhibit nonphysical deformation. This phenomenon is called as hourglassing. These modes occur when the integration points are unable to capture any stress or strain in the elements and display zero strain. The resulting stiffness offers no resistance to deformation mode resulting in instability. Some elements appear severely deformed whereas the overall mesh remains stable. It can be eliminated by providing a restraint which the element lacks after employment of reduced integration. The restraint must not have significant effect on the stiffness of element to other modes. The two main ways of hourglass control are the viscous damping and introduction of artificial stiffness. They are unable to control hourglassing in coarse mesh with large nodal forces (Belytschko et al., 2000;Belytschko, Ong, & Kennedy, 1984;Brodland, 1994;de Souza Neto, Peric, & Owen, 2011;Dhondt, 2004;Tejchman, 2013;Waszczyszyn, Cichon, & Radwanska, 2013).

Enhanced assumed strain (EAS)
EAS is a technique to avoid instability due to locking phenomenon during simulation. It is a suitable option for the elimination of volumetric and transverse shear locking. This method enhances the displacement related strain field by introducing additional degrees of freedom to the element Duchene et al., 2011;Duchˆene et al., 2013). An important feature of this method is that it does not increase the size of the stiffness matrix as no strain energy is added to the model but improves the flexibility in the element (Quak, 2007;Quy & Matzenmiller, 2007;Valente et al., 2012;Wisniewski, 2010).

Assumed natural strain (ANS)
ANS is incorporated to treat locking especially where results are distorted by the shear locking. It interpolates the strain components in two steps than the conventional single step method. In the first step, the interpolation is carried out in a classical manner at sampling points where locking is less likely to affect the result and depends on the order of quadrature employed. In the second step, the strain components are calculated from the sampling points to the integration points. It can sometimes calculate lower strain values than the actual values in an element Caseiro et al., 2014;Duchˆene et al., 2013).

Computation efficiency enhancing techniques
Care has to be taken in the simulation of SPIF process while selecting element size because of the small area of contact between the sheet and tool of a comparatively small diameter (Sena et al., 2011). The elements are required to be small enough to give accurate results but large enough to be computationally efficient. Certain techniques are also employed to reduce computational time.

Adaptive remeshing method
This technique was used to reduce the computation times while performing simulations for SPIF. The implementation of this technique refines the mesh near the tool and moves along the tool. The elements near the tool are divided in a fixed number of small elements. This phenomenon thus eliminates the requirement of fine mesh before simulation resulting in reduced computational time. However, the state of the fine mesh does not change if the elements are highly distorted to avoid irregularities in the result. This technique can modify the results as the force and stiffness matrix change due to variation in number of degree of freedom involved after refinement. It gives good results for SPIF simulation discretized by shell elements. It was found to reduce overall simulation time with implicit integration scheme (Lequesne et al., 2008). It was successful in geometry and shape prediction when employed in conjunction with the RESS element (de Sena et al., 2013). It was also employed in an explicit integration scheme to reduce computational time. The study reported errors due to explicit integration scheme were responsible for deviation from desired results (Bambach, 2014).

Substructuring
Substructuring is used to reduce computation time for implicit simulations for SPIF. The mesh is divided into several nonoverlapped structures categorized into plastic nonlinear substructures and elastic pseudo linear substructures. The part of sheet in contact with tool is categorized under the plastic nonlinear substructures calculated by the nonlinear integration schemes. On the other hand, the elastic deformation on the rest of the sheet is treated by the elastic substructures by assuming linear behavior. In addition, the internal degrees of freedom are eliminating by condensation of stiffness matrices and the internal force vectors. The iterations in this technique are carried out by updating the internal force vector by the multiplication of internal force vector with the stiffness. The stiffness matrices obtained in the condensed form are kept constant due to the linearity of the elastic deformation, whereas, the plastic substructures are updated nonlinearly. A good agreement with the experimental results was obtained in simulations employing this technique (Hadoush & van Den Boogaard, 2009a, 2009b.

Mass scaling
The simulation time in models employing explicit integration scheme increases due to the small time step size necessary to maintain stability. The computation cost can be reduced by the implementing mass scaling technique. This is done by increase in structure mass by increasing the density. This rise results in the increase of time step size of the elements less than the absolute size of DT2MS. Errors might result due to the increase in mass resulting due to mass scaling. A mass scaling factor of 10 can be used in the SPIF without significant effects on the result (Surech & Regalla, 2014).

Time scaling
Time involved in the forming of sheet metal by SPIF is relatively large. This is due to extremely long tool path required for deformation. Explicit integration schemes are able to model the processes that occur in a very short duration of time. Time scaling makes use of this feature by artificially increasing the tool velocity in the process to decrease the computational time. However, time scaling might produce dynamic effects in the solution and reduce accuracy. The kinetic energy of the system should be maintained at less than 10% of the internal energy of the system to obtain accurate results. The velocity of the tool can be increased to 40 m/s for the simulation of SPIF process (Surech & Regalla, 2014).

Discussion and conclusion
Explicit method offers reasonable accuracy considering the computation time while using the implicit method. However, they are unable to predict the springback which occurs due to the dieless nature of the process. Even though implicit method provides better prediction of geometry and thickness, it is a time consuming process. Techniques like adaptive remeshing are able to provide better results by reducing the time required for computation. Conventional shells have successfully been used for simulation of SPIF process. They are suitable to model the thin metal sheet when only the thickness and force predictions are essential. They are unable to predict the through thickness parameters such as through thickness shear and material flow. In order to predict the through thickness gradients the time consuming solids can be replaced by the solidshell elements. They can have unlimited integration points defined in the thickness direction for better prediction of through thickness gradients in comparatively less computation time. This statement can be supported by the constant enhancement of the solid-shell element and development of RESS element. Along with the EAS and selective reduced integration, the RESS has been integrated with ANS to prevent locking. This gives the element enough flexibility to avoid any kind of locking that might occur in the linear element. However, the underestimation of force is a major setback of this element. A proper parameter selection without major approximations can lead to better solution. However, the designer or engineer will either have to make a compromise of solution result or spend more resources for a more accurate result Jeswiet, Adams, Doolan, McAnulty, & Gupta, 2015.