Numerical Methods for the Modelling of Chip Formation

The modeling of metal cutting has proved to be particularly complex due to the diversity of physical phenomena involved, including thermo-mechanical coupling, contact/friction and material failure. During the last few decades, there has been significant progress in the development of numerical methods for modeling machining operations. Furthermore, the most relevant techniques have been implemented in the relevant commercial codes creating tools for the engineers working in the design of processes and cutting devices. This paper presents a review on the numerical modeling methods and techniques used for the simulation of machining processes. The main purpose is to identify the strengths and weaknesses of each method and strategy developed up-to-now. Moreover the review covers the classical Finite Element Method covering mesh-less methods, particle-based methods and different possibilities of Eulerian and Lagrangian approaches.


Motivation
Machining is used for manufacturing components in various industries including automotive, aerospace, defense, electronics, construction, electrics, shipbuilding and others. It is the most prevalent manufacturing operation in terms of volume and expenditure. Machined components are used in almost every type of manufactured product. However, it has been estimated that machining expenditure contributes to approximately 5% of the GDP in developed countries, while in the U.S. alone it translates into approximately $250B per year [47]. Common machining processes includes sawing, broaching, grinding, drilling, milling, turning and shaping (see Fig. 1). It is classified as a subtractive manufacturing method and it provides various advantages to the finished products such as closer dimensional accuracy, surface texture or finish, complex shaping, and required sizes.
Machining presents a challenging and exciting intellectual problem that has fascinated researchers and practitioners for decades. However, commercial codes employed for the numerical evaluation of cutting conditions to analyze the cutting forces, torques, and chip loads are using mainly leveraging methods that have not evolved during last decade. Without successful models, experimental testing will continue being the reference for the process development in practice. Although experimental testing is expensive, right now it is the reference source for the models used in the industry. The above considerations constitute the compelling reasons that pushed researchers to pursue different lines of research in the field of numerical simulations of metal cutting processes.

Objectives
During last years new revolutionary techniques appeared extending the capabilities of the classical numerical methods. The research in most of these techniques became stuck presenting a lot of shortcomings and a limited range of application. With the aim of finding a way to go further in the development of a numerical method applied to the modelling of chip formation we present here a review of the numerical methods that have been used or created for that purpose. The common denominator for the most of them is the Finite Element Method (FEM) from which other techniques has evolved. In this review, the methods under analysis that in this review are classified as methods supported by a mesh, including the FEM in a Lagrangian, Eulerian and Arbitrary Lagrangian Eulerian formulations (ALE), the Multi-material Initially all characteristics of the physical problem will be presented, giving a description of the problem, the theoretical approaches and the magnitudes of interest. Then, the current experimental methods to analyze the problem will be mentioned as well as the most common numerical method to solve it. Numerical challenges arising from the theoretical approach and solutions developed to resolve them, are going to be extensively explained in the frame of the FEM. They include the thermo-mechanical coupling, the material behavior, the treatment of the incompressibility, the thermomechanical contact and friction, adaptive discretizations, breakage and separation criteria for the chip. Finally, an evaluation of the extended characteristics of the different numerical techniques (finite element methods, particle methods and mesh free methods) will be presented telling how good they are resolving the identified numerical challenges.

Contents
The paper is organized as follows. After this introduction, Sect. 2 presents the chip formation fundamentals based on classical models for the different machining techniques in a descriptive manner. Furthermore, in Sect. 3, all the numerical ingredients that are necessary to the numerical simulation of machining processes using the Finite Element Method are explained in detail. In Sect. 4, improved Eulerian Formulations capable to reproduce chip formation in orthogonal cutting conditions are explained. In Sect. 5, most important meshless methods that have been applied in the numerical simulation are presented. Finally, in Sect. 6, some conclusions are presented and a discussion provided on the need for improvements of the available numerical techniques.

Chip Formation Fundamentals
The purpose of this section is to explain the usual types of machining processes and parameters for their description and to bring together observations on the shape of chips. The role of mechanics in this context is more to aid the description than for the chip formation prediction. However, the chip formation in all machining processes (turning, milling, drilling and so on) can be described in a common way, a tip of a cutting tool removing a layer of material and creating a chip (see Fig. 2). Therefore the physics and numerical predictions explained in this paper relate to any process.
For modeling of material removal of material caused by a machining process (chip formation), some important issues are: how the thicknesses of chips vary with tool geometry, how the thicknesses of chips vary with the friction between the chip and the tool, which influence has the work hardening behavior of the machined material, how the forces on a tool during cutting may be related to the observed chip shape, how the friction between the chip and the tool and the plastic flow stress of the work material will vary. In order to answer these questions we must start describing the physical process of chip formation and its characteristic parameters.

Physical Description
In a cutting process, the tool cutting edge penetrates into the workpiece material, which is thus plastically deformed and slides off along the rake face of the cutting edge. This is the fundamentals of a chip formation process. Fig. 1 The three most common types of machining processes. a Turning, b milling and c drilling [36] The main terms are the cutting speed v cut , the chip flow v chip , and the relative directions between them determine if the formation of the chip is orthogonal or non-orthogonal, see Fig. 2. Considering the geometrical description, the general used terms are: rake angle (positive by convention), the clearance angle , inclination angle , major cutting edge angle , feed f and depth of the cut d.
Moreover, the definition of some terms changes a little bit between machining processes, for example the definition of the feed f. In turning and drilling, the feed f means the distance moved by a cutting edge in one revolution of the work; in milling it means the distance that the cutter advances across the workpiece in the time taken for each cutting edge to move to the position previously occupied by its neighbour. Feed f and depth of cut d always refer to displacements from the point of view of machine tool movements. They are described in Fig. 2.
Usually the most appropriate terms, from the point of view of the chip formation process, are the uncut chip thickness f ′ and the cutting edge engagement length d ′ . These terms closely related to feed and depth of the cut and they coincide in some particular cases of turning, in orthogonal chip formation with a 90 • major cutting edge angle.
Usually the processes in chip formation can be examined within the orthogonal plane (Fig. 2a), because the essential part of the material flow take place within this plane. This means that one can assume that the deformation is two-dimensional. The two-dimensional deformation is only disturbed at the edges of the cross section of the undeformed chip at the free surface, and in front of the cutting edge corner. This happens because there is material flow at an angle towards the orthogonal plane. However it is not a significant behavior for the characterization of the cut.
Depending on the deformation behavior of the workpiece material, there are different mechanisms of chip formation with either continuous or discontinuous chip flow. This mechanisms are presented in the next section.
Before presenting the different type of chips there are some important mechanical parameters that describe the process in terms of forces. The tool-workpiece interaction is usually characterized by the cutting and thrust forces. However, the resultant of the contact force between the tool and the workpiece can be decomposed in the cutting force and the thrust force. Moreover, the cutting force F c is the component of the resultant in the direction of the cutting speed (which is the relative speed between the work and the tool). On the other hand, the thrust force F t is the component perpendicular to this direction, in Fig. 3 is decomposed in two forces, the feed force F f and the passive force F p .

Fig. 2
Orthogonal and non-orthogonal chip formation [54] Fig. 3 Cutting force F c , feed force F f , and passive force F p Usually the direction of the resultant force can be predicted if the direction of the workpiece plastic flow is known. There exist different models of deformation zones, one of them is presented in next section, where a primary shear plane of the chip formation is defined and provides one of the main characteristics of the cutting conditions.
Other characteristics of the tool-chip interaction are the chip thickness, the chip-tool contact lengths and the resultant chip-tool contact pressures. Furthermore, the chip-tool contact pressures will determine the resultant force between the tool and the workpiece. Finally, the collection of all these relationships induce to the formation of different type of chips.

Mechanisms of Chip Formation
Depending on the workpiece material and the cutting conditions, the following mechanisms of chip formation can be distinguished ( Fig. 4): continuous, lamellar, segmented (or serrated) and discontinuous chip formations.
The Continuous chip formation is produced when the chip slides off along the rake face at a constant speed in a stationary flow. Continuous chip formation is promoted by a uniform, fine-grained structure and high ductility of the workpiece material which favor low undeformed chip thickness. The Lamellar chip formation is a continuous, periodic chip formation process similar to pure continuous chip formation. However, there are variations in the deformation process that cause more or less significant cleavages or even concentrated shear bands. The Segmented (or Serrated) chip formation is the discontinuous formation of a chip with still more or less connected elements. Furthermore, the chip use to have a higher thickness and present significant variations in the degree of deformation along the flow path. Finally, the Discontinuous chip formation occurs if the plastic ductility of the workpiece material is very low or if predefined slide paths are formed due to high inhomogeneities. Parts of the workpiece material are ripped out of the compound material without significant deformation. The workpiece surface is then rather produced by the ripping out process in the chip formation than by the tool traces.
The most common chip formation process under analysis is the continuous chip formation. This process can be described using a model of five deformation zones (Fig. 5). Furthermore, the main part of the plastic deformation takes place in the primary shear zone where a plastic band is formed due to shear deformation. In the secondary shear zones, one in front of the rake face, another one on the tip and a third one on the flank face, the workpiece material is additionally deformed under the influence of high friction forces. A stagnant zone (zone with high pressure from all sides) develops in front of the cutting edge (on the tool tip). In this zone built-up edges can occur. Build-up edges are formed by particles of the workpiece material, which adhere to the rake face and to the cutting edge. This is also the location where the actual separation of workpiece material takes place. Furthermore, minor plastic deformations occur in the preliminary deformation zone, see Fig. 5. This zone has an essential influence on the penetration depth of plastic deformations in the workpiece and so on the external workpiece zone.

Chip Formation Analysis
The most important experimental methods for the analysis of the chip formation process are the interrupted cut and the micro-cinematography. The basic principle of the interrupted cut is an abrupt separation of the tool from the workpiece. The current status in the deformation process is thus "frozen", can be subject to metallographic preparation and analyzed under the microscope [94]. In contrast to the interrupted cut method, micro-cinematography allows  [93] for examining the chip formation while the process is still running. For this purpose, a polished and etched specimen (black-and-white structure) is pressed against a fused quartz glass plate and cut by orthogonal turning. The process can be observed through the fused quartz glass plate and magnified using a microscope [102]. The analyses provide information on the mechanisms of chip formation, the plastic deformations in the chip formation zone and the position of the shear plane. They provide the basis for a calculation of the kinematic, mechanical and thermal conditions in the chip formation zone.
The Finite Element Method (FEM) is one of the most important methods to visualize the deformation process in front of the cutting edge, and to analyze the chip formation process, and the material behavior in the working zone. The numerical modelling of chip formation is a challenge. Except the physical phenomena explained above two more challenges need to be addressed. The first one is to provide accurate data to the model. The second is to actually choose the correct approaches and formulations for the treatment of thermo-mechanical coupling, the material behavior and incompressibility, and the treatment of contact and friction. Furthermore different discretization of the space, adaptive methods and transfer operators, separation criteria for segmentation and breakage of the chips, are used for approximating a solution. The combinations that have already been tried by researchers are numerous. These approaches are discussed in the following sections. In general, from the definition of the model and the cutting conditions (cutting angles, feed and depth of the cut, and velocities), numerical simulations try to reproduce the different cutting mechanisms described in Sect. 2.2. Parameters characterizing the chip formation, like cutting forces F c , chip flow, chip thickness and chip continuity will be obtained as results of these analysis.

The Finite Element Method
In this section we will start with the description of the classical finite element method FEM, for solid mechanics and large deformations. The method is commonly based on the Lagrangian description from continuum mechanics, but also Eulerian approaches can be used. After a general review of the well-known method, we will put emphasis on element technology and the wide range of techniques used for the treatment of the incompressibility. After that, time integration schemes and contact algorithms will be addressed. Finally, the adaptive remeshing, transfer operators, as well as the criteria for chip separation will be presented.

Problem Formulation
Finite Element Lagrangian formulations embeds a computational mesh in the material domain and solves for the position of the mesh at discrete points in time. As a consequence the Finite Element Lagrangian formulation is related to the problem of mesh distortion. The chip formation process occurs due to large deformation of the material where the mesh distortion becomes a critical factor. If the distortion is too large the calculation process can even be impossible to continue because Jacobian determinants become negative at  [102] some integrations points. At times, Finite Element Lagrangian formulations use a criterion to separate the chip from the workpiece. Those criteria included element deleting based on a geometrical distance of the tool tip to closest workpiece element, plastic strain and strain energy density. Other times, Finite Element Lagrangian formulations are used with mesh adaptivity and automatic remeshing, and as a consequence this strategy does not require a chip separation criterion. Using remeshing means that the fields of state variables have to be mapped from the old mesh to the new mesh. This mapping is not a straight forward task and introduces some numerical diffusion to the state variables. This technique has been successfully applied in simulations of continuous and serrated chip formation. A lagrangian description of motion and adaptive remeshing was used to simulate orthogonal cutting in [58,67,86] (see Fig. 6). Instead, a Lagrangian formulation including node separation criterion to predict chip piece separation was used in [51,88,90] (see Fig. 7).
Finite Element Eulerian formulations have been used by many authors to simulate continuous chip formation at steady state. Finite Element Eulerian formulation avoid the problem of mesh distortion but needs a predefined chip shape to develop the numerical simulation, while Finite Element Lagrangian formulation is able to predict chip formation from incipient to steady state. A proper assumption of Numerical modeling of metal cutting processes using a Lagrangian FEM and a node separation technique [107]. When the separation criteria are fulfilled elements are disconnected. The mesh and the material domain are the same Fig. 7 Numerical modeling of metal cutting processes using a Lagrangian FEM and adaptive remeshing [58] the chip shape is very difficult to obtain since it depends on many factors related especially to the Lagrangian description of the particles movement. On other hand, Finite Element Eulerian formulation the material flows through the fixed mesh that facilitates large deformation calculation. The main disadvantage of Eulerian formulations is that is no easily adaptable for modeling the unconstrained flow of the material as the chip evolves during the process. As a consequence Finite Element Eulerian formulations cannot simulate serrated and discontinuous chip formation. An example of a Finite Element Eulerian formulation applied to the numerical simulation of metal cutting is presented in [76].
In order to avoid the disadvantages of Finite Element Lagrangian and Eulerian formulations, other computational techniques have been investigated. One of them is the Finite Element Arbitrary Lagrangian Eulerian (ALE) formulation in conjunction with adaptive mesh techniques [32,33,63,75]. The ALE formulation combines the best features of pure Lagrangian analysis (in which the mesh follows the material) and Eulerian analysis (in which the mesh is fixed and the material flows though the mesh). In the ALE framework mesh motion is independent of material motion, for that reason high quality finite element meshes are preserved during the numerical simulation of machining process. Another advantage is that the ALE formulation does not need a criterion to separate the chip and the workpiece. Generally, ALE formulation is computationally cheaper than a Lagrangian formulation because remeshing is not needed only node relocation. ALE formulation can be used to simulate chip formation from incipient to steady state, but the problem is to define a mesh motion scheme in order to preserve a high quality finite element mesh during the simulation (the improvement of mesh quality is not guaranteed using only a node relocation technique). Numerical simulations in 3D using ALE are difficult to carry out, because the tracking of the mesh motion in needed to preserve a high quality mesh is more difficult in 3D than in 2D. A detailed information on use of ALE formulations in modeling metal machining are presented by [63,75,32] (see Fig. 8).

Thermo-Mechanical Coupling and Material Modelling
One of the most important aspects in the definition of the problem is the description of the material behavior. The constitutive modeling of the material is something Fig. 8 Numerical modeling of metal cutting processes using an ALE formulation [107]. Note that the mesh follows the material but keeps the same structure and connectivities inherent in the numerical treatment of the problem. Several material models have been developed for the study of chip formation, see ( [4,56,80,100]). As was mentioned in the chip formation analysis Sect. 2.3 it is known that plasticity zones appear in the material workpieces inducing the characteristics of the chip. The plastic deformation generates heat and it must be considered for an accurate reproduction of the problem conditions. Furthermore, the cutting interaction is also a heat source of the problem. Accurate numerical simulations must account for the thermo-mechanical coupling and for the heat transfer between tool and workpiece. Most common approaches for the thermo-mechanical analysis can be found in the literature, see ([26, 58, 86, 87, 100]). Under last mentioned conditions, numerical techniques must also solve the issues arising from the modeling of the plasticity. One of the major problems is the numerical treatment of the incompressibility. Material turns to incompressible when plastic flow occurs in deformation zones. This numerical problem has been solved using different techniques which at the same time, have conformed the characteristics of some methods of finite elements themselves. Section 3.3 presents a extended review of the different numerical techniques developed for the treatment of the incompressibility constraint due to plasticity. Von Mises plasticity model is the mostly used material model for reproducing metal-type behavior. This is a well known type of associative plasticity where the yield surface is defined by a function coming from the definition of the deviatoric part of the strain energy (see Eq. 1).
where �dev( )‖ is norm of the deviatoric part of the Kirchhoff stress, ( y + )(ē p ) is the mathematical definition of the stress flow which depends on the plastic deformation ē p .
From this basic definition of the elastic ( Φ < 0 ) and plastic zones ( Φ = 0 ), other phenomena are added to the model: temperature dependency ( ), mathematical definitions of the stress flows producing material hardening or softening ( y + ) or models for the material failure. The description of these effects are usually added to the definition of the yield surface, so Φ , ( y + )( ,ē p ) .
Some of the most used hardening models for thermoplasticity ( y + )( ,ē p ) are the Simo model (for steel) [87] the Johnson-Cook models(for titanium) [7,8,49,55], the micromechanical models [38,59] or the crystal plasticity models (dislocation density based models) [80], considering material anisotropy (for aluminium and other alloys), and others. Some of them consider also a relevant aspect, the dependency on the evolution on the plastic strain rates ̇ē p , so ( y + ) ,ē p ,̇ē p . This is a crucial aspect for the correct modeling of the shape of the chip during the machining process. The speed of the cut is a variable that is going to change the chip morphology basically because the material behaves in a different manner if the strain rate changes. Those models will be used to define the chip separation criteria explained in Sect. 3.7.
There are also models that consider other global or local effects, some of them try to take into account the dislocations produced in the material microstructures or nanostructure. Using different mathematical formulations this can be described without the need of using a multi-scale analysis. One of the more important effects to take into account is the generation of cracks in the stressed regions, coupling plasticity with damage models allows for an easier determination of the shear bands and the sliding surfaces that will generate the breakage of chips during the cut. A brief explanation of the most common breakage criteria is explained in 3.8.

Numerical Treatment of the Incompressibility Constraint Due to Plasticity
The most common finite elements used in the numerical simulations of metal cutting are the following: a plane strain quadrilateral isoparametric finite element used in [90,107], a 6 noded isoparametric triangular elements used in [58,86], an enhanced four node quadrilateral with 1-point quadrature used in [67] and a 3 noded linear triangle plus the Average Nodal Pressure formulation to deal with the incompressibility constraint used in [26]. The 6 noded isoparametric triangular element presented in [58], is one of the elements used in the commercial software AdvantEdge [112] and the 4 noded quadrilateral with reduced integration is used in Deform [111]. Those softwares are the most common numerical tools used in industry in the numerical simulation of metal cutting processes. In Table 1 a summary of the main characteristics or these commercial tools is presented. Also, a number of different finite elements have been developed to improve the poor performance of linear triangles and tetrahedral under incompressible and nearly incompressibility conditions. These finite elements can be classified in four groups mainly: (1) Mixed Enhanced Element, (2) Pressure stabilized finite elements, (3) Composite pressure fields and (4) Average Nodal Pressure/Deformation Gradient. The following lines present a summary about the advantages and the disadvantages of each of the improved linear triangle and tetrahedron, and in case it is available, a reference where apply the improved element in the numerical simulations of metal cutting processes.

Mixed Enhanced Elements
Enhanced Strain Technique, essentially consists in augmenting the space of discrete strains with local functions, which may not derive from admissible displacements. A suitable choice of these additional modes can improve the numerical performance of low-order elements and, more importantly, it can greatly alleviate the well-known volumetric locking phenomenon in the nearly incompressible regime.
Most of the schemes taking advantage of the Enhanced Strain technique have been designed in connection with quadrilateral elements, because linear displacement finite elements enriched with an enhanced strain locks. Instead, the mixed linear displacement linear pressure finite element enriched with enhanced strains is able to deal with the incompressibility constraint. Literature review [6], remarks that the straightforward extension of the Enhanced Strain approach to large deformation problems generally leads to unstable methods, representing a disadvantage of Enhanced Strain Technique. To the authors knowledge, Mixed enhanced finite elements have not been applied in the numerical simulations of metal cutting.

Pressure Stabilization
The pressure field in mixed linear-displacement/linearpressure finite elements presents unphysical oscillations when used in the numerical simulations of incompressible or nearly incompressible materials. Mathematically, it means that equal order interpolation for displacement and pressure does not satisfy Babuska-Brezzy condition. In order to remove these undesirable oscillations, a literature overview shows four different strategies mainly: characteristic based split (CBS) [24], Finite Calculus (FIC) [64], Orthogonal Subgrid Scales (OSS) [22,23,85] and the Polynomial Pressure Projection (PPP) [12,28].

The Characteristic Based Split
The characteristic based split (CBS) was originally developed in the field of fluid mechanics [24]. This method is based on the introduction of an artificial compressibility into the mass conservation equation, in such a way that the final results do not depend on the artificial compressibility. The other main ingredient of CBS is the fractional step method used in the time integration of momentum balance. This fractional step proposes a split of momentum equation in two equations such that its sum is equal to the balance of momentum equation. The equation split is equivalent to split the velocity update in a time step into deviatoric and hydrostatic components. In summary, CBS algorithm uses four main steps: (1) Compute the velocity update using an explicit time integration scheme of the equation of balance of momentum; in this time integration schemes hydrostatic forces are not taken into account; (2) Using the balance on mass, calculates the nodal pressure and (3) Using the velocity obtained in (1) and the gradient of the pressure field obtained at (2) update the velocity field. (4) Given the updated velocity using an explicit integration scheme to obtain the value of nodal displacements and update nodal positions. After this four steps, the next time step start. One advantage of CBS algorithm is the possibility to evaluate the pressure in a complete explicit way, but at the same time CBS allows to solve the pressure using an implicit scheme, in case it is needed. As far as authors know, there is no a complete CBS implicit scheme for velocity and pressure, so does not matter if the pressure is integrated implicitly CBS algorithm is conditionally stable, and as a consequence there is a restriction in the maximum allowed time steps.

The Finite Calculus
The Finite Calculus (FIC) method [64] is the based on the satisfaction of the balance of momentum and mass conservation in a domain of finite size. However FIC retaining higher order terms in the Taylor expansions used to express the different terms of the differential equations over the balance domain. The modified differential equations contain additional terms with the main function to suppress the unphysical oscillations of the pressure field. The mixed linear displacement/linear pressure tetrahedral could be used with implicit, explicit and semi-implicit time integration schemes. This represents an advantage in comparison with other modified tetrahedral finite elements such as the CBS [24]. Also, FIC method needs 5 degrees of freedom per node (2 displacement, 1 pressure, 2 projected pressure gradients) in case of linear triangle and 7 (3 displacement, 1 pressure, 3 projected pressure gradients) in case of tetrahedral. However, the number of degrees of freedom per node represents a disadvantage of FIC, in terms of computing time, required memory and an extra nodal variable to transfer between remeshings. Another disadvantage of FIC is that, the term added to the mass conservation equation depends on the shear modulus, the mesh size and a constant that is problem dependent.

Orthogonal Subgrid Scales
Orthogonal Subgrid scales (OSS) was applied in the field of solid mechanics in [22,23,85]. Orthogonal Subgrid scales approach is based on two main ingredients: (1) a mixed equal order interpolation of the pressure and displacement fields and (2) a decomposition of the unknowns into resolvable and sub grid scales orthogonal to the finite element space. However, the idea behind Orthogonal Subgrid scales is to approximate the effect of the continuous solution which cannot be captured by the finite element solution which is the cause of volumetric locking. Moreover, the main purpose of Orthogonal Subgrid Scales is to define a strategy to overcome the requirements of Babuska-Bressi conditions and in consequence make possible the use of equal order continuous interpolation for displacement and pressure. As a consequence of adding a subgrid scale displacement, extra degrees of freedom are added to a node in comparison with a mixed formulation using simplicial elements with equal linear interpolation for displacements and pressure. Furthermore, a term that depends on mesh size, shear modulus, and a constant that is problem dependent is added to the mass conservation equation. Then, subgrid scales needs 5 degrees of freedom per node (2 displacement, 1 pressure, 2 projected pressure gradients) in case of linear triangle and 7 (3 displacement, 1 pressure, 3 projected pressure gradients) in case of a tetrahedron. In the framework of implicit dynamics, a staggered scheme in which the displacement and pressure are solved implicitly, while the pressure gradient is solved explicitly is usually used as proposed in [23]. As a consequence the pressure gradients degrees of freedom added represent a minor cost in terms of computing time.
It is important to remark that the terms added to the balance equations using FIC and subgrid scales to overcome Babuska-Brezzi conditions are exactly the same. Being the difference between FIC and subgrid scales the idea of how the terms are obtained in order to stabilize the finite element solution.

Polynomial Pressure Projection
Polynomial Pressure Projection (PPP) was initially formulated and applied to stabilize stokes equations in [12,28]. PPP is based on two ingredients. First, use a mixed equal order the pressure and velocity fields and secondly, a L 2 pressure projection. FIC and OSS introduce the projection of the pressure gradient onto the displacement space as a new dependent variable, and use the difference between these two fields to relax the continuity equation, while PPP uses a projection on a discontinuous space and as a consequence can be implemented in an elementary level. FIC and OSS use mesh dependent and problem dependent parameters while PPP does not need. PPP has been applied recently in the numerical simulations of metal cutting processes (see [77-80, 82, 83]).

Composite Pressure Fields
These finite elements enforce a constant pressure field over a group of triangles or tetrahedrals to reduce pressure constraints. The most representative finite elements are F-Bar [60,71] and Composite Triangles [19,37] . The F-bar is another strategy to deal with the incompressibility constrains using linear triangles. F-bar formulation was proposed in [60] in the framework of implicit dynamics and in [71] in the framework of explicit dynamics. It relies essentially on the relaxation of the excessive volumetric constraint typical of low order elements through the enforcement of the incompressibility constraint over a patch of simplex elements. An important aspect of the present method is that it preserves the displacement-based format of the corresponding finite element equations. At the same time, this method presents an unconventional stiffness format that stem from the fact that the internal force vector of a particular element depends on the nodal displacements of all elements of its patch, breaking the typical elementary assembly of the internal force vector and the stiffness matrix.
A triangular element in which a six-node triangle is constructed from four 3-node triangles with linear displacement fields in each subtriangle and a continuous linear strain field over the assemblage is presented in [19]. They have called this finite element a "composite triangle (CT)". This element presents some advantages in terms of contact search and imposition in comparison with 6-nodetriangle and furthermore this element is locking free in comparison with three node triangle. Furthermore, this element does not satisfy BB condition. An improved CT triangle in order to satisfy BB conditions using a constant volumetric strain and a linear deviatoric strain over the six-node finite element is presented in [37].

Average Nodal Strains/Average Nodal Stresses
It uses the computation of the average volumetric-strain/ volumetric-stress or strains/stresses at nodes based on surrounding triangles or tetrahedrals. Then, the elementary volumetric-strain/volumetric-stress or strains/stresses are equal to the average of the nodal values that belongs to the element. Average Nodal Pressure (ANP) was presented in [13,14] in the framework of explicit dynamics and by [3] in the framework of implicit dynamics. ANP is applied to a simple linear tetrahedron element that can be used in applications involving nearly incompressible materials or incompressible materials modeled using a penalty formulation. The element prevents volumetric locking by defining nodal volumes and evaluating average nodal pressures in terms of these volumes.
Average Nodal Pressure (ANP) defines the nodal volume as follows: where V (e) is the element volume, ndim is the dimension of the problem to solve and V a is the nodal volume. This definition of nodal volume, allows to define a nodal volume change ratio and as a result a nodal pressure. The average of nodal pressures is used as the modified elementary pressure. This definition of nodal volume reduces the volumetric locking tendency of linear triangles and tetrahedral and allows an accurate prediction of deformed shapes and forces.
However, ANP formulation is found to produce considerable checkerboard-type hydrostatic pressure fluctuations, which limits the range of applicability of ANP pressure formulations. It is important to remark that ANP work well for volumetric locking but present locking due to bending [72]. ANP has been applied in the numerical simulation of high speed cutting in [26]. The critical time step of ANP formulations imposed by stability is more or less 7 times greater in comparison with CBS formulations [14]. In terms of computing time it represents a great advantage of ANP formulations.
Two of these formulations are the Node Based Uniform Strain Elements (NBUSE) [28] and the Average Nodal Deformation Gradient (ANDG) [15] which present linear triangles and tetrahedra that are free of volumetric and shear locking. These average gradients are defined as follows: where V e is the element volume, ∇ (e) is the element displacement gradient, ∇ a is the nodal displacement gradient, (e) element deformation gradient and a is the nodal deformation gradient. Then using ∇ (e) ∕∇ a and the constitutive equation, nodal stresses are calculated and finally, modified elementary stresses are calculated as the average of nodal stresses. As reported in [15,28] these formulation can lead to the presence of non-physical low-energy modes. An improved nodal deformation gradient based on Streamline Upwind Petrov-Galerkin (SUPG) which remove the nonphysical energy modes was proposed in [15]. A stabilization strategy to remove the unphysical energy modes in the framework of implicit dynamics was proposed in [72], even so, the pressure field in some examples presents unphysical oscillations. The great advantage of average nodal formulations is that they offer locking free elements without increasing the number of degrees of freedom per node. NBUSE has been applied in the numerical simulation of metal cutting processes in [34].
Another relevant technique based on average quantities is the Mixed Discretization Technique (MD) [57] which is based on the following ingredients. (1) Mesh the solid body in quadrilateral or hexahedral zones; then divide each quadrilateral or hedrahedral into triangles or tetrahedrons. (2) The deviatoric behavior is defined on an elementary basis (triangle/tetrahedron), while the volumetric is averaged over a zone (quadrilateron/hexahedron). Then, an improvement of (MD) is presented in [27], the authors of that work call their formulation Nodal Mixed discretization (NMD). Moreover, the main advantage of (NMD) is that the average of the volumetric behavior is carried out in a nodal basis rather than in a zone basis (quadrilateral/ hexahedron), it implies that only a mesh of triangles or tetrahedra is needed. NMD uses a nodal volumetric strain rates, defined as weighted average of the elementary surrounding values. Then, a modified elementary volumetric strain rate is defined as the average of the nodal values. However, the difference between NMD and ANP is that in NMD the constitutive model is called on an element basis, while in ANP the constitutive model is called in a nodal basis for the volumetric behavior. Due to the similarities between ANP and NMD, is expected that NMD presents checkerboard-type hydrostatic pressure fluctuations, being a disadvantage of both formulations. In case the hydrostatic part of the stress tensor depends linearly on the determinant of the deformation gradient, ANP and NMD are exactly the same.

Time Integration Schemes
The Finite Element Method allows different time discretization schemes. The most common are the implicit and explicit time integration schemes. Each of them has its advantages or disadvantages. Implicit schemes are unconditionally stable; it means that there is no restriction on the time step used in the numerical simulation. In implicit formulations, mechanical problem can be solved in a static o dynamic way. Furthermore, implicit formulations can be used with standard and mixed (displacement/ pressure) finite elements. However, implicit schemes needs the solution of a linear system certain number of times each time step. Usually, the solution of the linear system represent most of the computing time. Implementation of a new constitutive equation is a more difficult task in an implicit scheme, due to the requirement to implement the algorithmic constitutive tensor. Moreover, in some cases an implicit scheme does not converge, due to the high nonlinearities involved in the problem. Finally, contact conditions decrease the size of the time step used in the numerical simulation and as a consequence increase the computing time of an implicit scheme.
Explicit formulations solve the mechanical problem in a dynamical way. The solution of each time step in an explicit scheme is simple and computationally efficient, provided by the use a lumped mass matrix in the simulation. Explicit schemes do not need the solution of a linear system; this is an advantage if the numerical solution is done using parallel computing. Implementation of a new constitutive equation is an easier task; it allows implement simple or complex constitutive equation without a large effort. Explicit scheme are conditionally stable, it means that the time step used in the simulations should be less or equal than a given critical time step, the critical time step correspond to the time that take for a sound wave to travel through the small finite element of the mesh. In case of an elastic material, the critical time step depends on the mesh size Δx compressibility modulus , Poisson ratio , density of the material and a constant that depends on the finite element used.
The restriction imposed on the time step by the explicit schemes, allows to conclude that for numerical simulation which involves long period of computing time or low speeds, implicit schemes are more favorable in comparison with explicit schemes. On the contrary, when velocities are high and the contact conditions are complex, is necessary to decrease the time step used in implicit formulations. In case explicit formulations offers efficient calculations, with a really interesting computing time. One example in cutting mechanics in which explicit schemes are more efficient that implicit scheme is high speed cutting. Now, the question is: At what cutting speed, explicit schemes are more computationally efficient than implicit schemes? It is important to mention that in the literature there is no a comparison between explicit and implicit time integration schemes, which shows under what condition one scheme is better than the other. In the literature, implicit schemes have been used in chip formation simulation in [7,86,88,90] and explicit schemes in [26,58,67].
Also, there are some mixed schemes in which the hydrostatic part of the balance of momentum is integrated implicitly and the deviatoric part is integrated explicitly. Some examples of mixed time integration schemes,related to the treatment of the incompressibility, are: Characteristic Based Split [84] and Finite Calculus [64]. These strategies have not been applied in the numerical simulations of metal cutting yet. Some commercial software like AdvantEdge uses explicit time integration schemes, while other software like Deform uses implicit time integration schemes.

Contact Algorithms
Modeling the complex thermo-mechanical phenomena that takes place at the tool chip interface is of paramount importance, because numerical results like feed force and contact depends strongly on an accurate modeling of the thermomechanical contact at the tool-chip interface. For that reason, in the numerical simulation of metal cutting processes is necessary to use accurate, robust and computational efficient contact algorithms. Contact problems using Finite Elements implies two basic challenges: First, the way in which the contact constraint is imposed and second, the contact detection strategy. Contact constrains are mainly imposed using the penalty approach [39,42,68,74,104] or the Lagrangian multipliers [74,104]. A mixed penalty-lagrangian formulation which uses the advantages of penalty and lagrangian formulations is presented in [87]. Several contact formulations developed up to now, enforce the contact constraints at specific collocation points (contact detection strategy). However, the most common strategy is the node-to-segment approach developed by Hallquist et al. [39] Its main idea is that a node located on the slave surface must not penetrate the opposing master side segment. This approach can be applied in a single and two pass algorithm. In a single pass algorithm only nodes on the slave side are checked against penetration into the master segment, and the nodes on the master side are free to penetrate the slave segments, while in the two pass algorithm, the nodes on the slave surface are checked against penetration into the master segment and the nodes on the master surface are checked against penetration into the slave segment. Both searching strategies have disadvantages because one pass algorithm allows penetration of master nodes into slave segments and do not pass the patch test, and the two pass is prone to lock due to overconstraining of the displacements on the contact surface, but it pass the patch test. Also, the node to segment approach has been extended to thermomechanical contact by Wriggers, Zavarise, and Miehe in [105,106,108].
Recently, a contact algorithm strategy based on following ingredients: (1) The continuity constraints imposition along the entire coupling boundary in a weak integral sense and (2) the use of segment-to-segment discretization strategies based on the so-called Mortar Method was presented in [30,95]. In contrast to classical node-to-segment formulations, in the segment to segment discretization the contact constraints are not imposed pointwise at a finite number of slave nodes but are defined along the entire contact boundary and therefore a more complete coupling between the degrees of freedom of the contact surfaces is obtained. However, an extension of the mortar method to thermo-mechanical dynamic contact problems including frictional heating and thermal softening effects at the contact interface was present by Heber and Wohlmuth in [44].
Oliver et al. [40,62] propose the Contact Domain Method (CDM). In this method, contact constraints are enforced using a stabilized Lagrange multiplier formulation based on an interior penalty method (this allows the condensation of the introduced Lagrange multipliers and leads to a purely displacement driven problem) and contact detection strategy is done with a fictitious intermediate region connecting the potential contact surfaces of the deformable bodies (this fictitious intermediate region is built using Delaunay triangulation). Oliver identifies the following advantages of CDM in comparison with other contact algorithms: (1) The solution does not dependent on the choice of slave and master sides, as the contact pairing is uniquely defined via a constraint Delaunay triangulation (2) the performance of the contact domain method (CDM) is superior to classical node-to-segment formulations and comparable to mortar based contact algorithms. The Contact Domain Method (CDM) has not been applied in the simulations of chip-formation problems.
The summary presented above about contact imposition and contact search is focused on implicit dynamics. In case of explicit dynamics special procedures have been developed [9,39], such as the momentum related techniques in which modifications are made to acceleration, velocities and displacements. The main goal of these modifications is to avoid the penalizing effect of the time step introduced by the high stiffness, associated with penalty approaches. Precisely, this strategy has been used in the numerical simulation of metal cutting processes by Marusich and Ortiz [58]. Bruchon et al. [17], propose the use the metric properties of the distance (or gap) function between two bodies in the formulation of contact problems. In this formulation, the vectors normal to the contact surfaces are defined through the gradient of this distance function. This contact strategy can be applied in explicit and implicit frameworks. However, the contact strategy presented by Bruchon has been applied in the numerical simulation of high speed metal cutting processes by De Micheli and Mocelin in [26]. Also, Sekhon and Chenot [86] present a very simple contact strategy applied to the numerical simulation of metal cutting processes using a flow formulation, such that if a boundary node is found to lie inside the tool, the length of time step is decreased in such a way that the node would lie on the tool face. Then, if the node at the start of the time step is in a compressive stress state, then the node is restricted to move along the tool face.
An overview about the contact algorithms available in the literature, shows that contact modeling using finite elements is not a simple a task, more research is need in order to develop a robust, efficient and general contact algorithm. However, the best algorithm able to predict the contact phenomena at the tool chip interface can be chosen based on the following parameters: (1) computationally efficient, (2) exact satisfaction of contact constraints, (3) decrease matrix ill-conditioning and (4) no extra degrees of freedom due to contact constraints.

Adaptive Remeshing, Error Estimator and Transfer Operators
In the numerical simulations of metal cutting processes large deformation, material and geometrical nonlinearities are present. Due to this reason, mesh degeneration through the numerical simulation is present, the first approach to tackle this problem was using predistorted finite element meshes [51,88], some of these references mention the limitations of predistorted meshes in the numerical simulation of metal cutting processes. Another approach to deal with this problem is adaptive remeshing. In numerical metal cutting simulations the magnitude and distribution of error estimators evolve during the incremental solution, showing that is necessary to refine the mesh where high gradients are taking place or de-refine the mesh where errors estimators are small in order to preserve a bounded computational cost. Furthermore, remeshing is used to preserve an adequate element shape and to predict cracking in the numerical simulation of serrated and discontinuous chip formation. Moreover, the elements of the boundary in contact with the tool are prone to distort and interfere with the tool, such an interference can lead to losses of volume of the workpiece and the undesirable effect of enlarging the tool radius, for that reason contact boundary conditions is one important remesh indicator. Generally, three steps are needed to remesh the workpiece: (1) Calculate error estimators and distortion metrics. They will determine if a new mesh is needed and then (2) create a new finite elements mesh, and (3) transfer data from the old to new mesh. Over the years, several error estimators have been proposed for elasto-plastic problems based on mathematical foundations and physical considerations. Zienkiewicz-Zhu [109,110] error estimator, calculates for each node an improved stress and defines the error as the difference between this stress and the one calculated by the standard finite element procedure. Ortiz and Quigley [73] propose an adaptive strategy based on equi-distributing the variation of the velocity field over the elements of the mesh. Marusich and Ortiz [58] used an adaptive criterion based on the equi-distribution of plastic power, the plastic work criteria has been applied in [58] to the numerical simulation of metal cutting processes. Lee and Bathe [52] propose a point wise indicator for error in stresses and plastic strain increments, based on the difference between smoothed (stress/plastic strain) and the (stress/plastic strain) at gauss points. Peric et al. [69] present an error estimator based on the projection and smoothing of plastic power rate and the rate of fracture. Furthermore, this method is not only able to capture the progression of plastic deformation, but also provide refined meshes at regions of possible material failure. Moreover, this error estimator has been applied in the numerical simulation of high speed machining in [67,69]. Micheli and Mocelin [26] presented an adaptive remesher coupled with mesh boxes used to mesh very accurately the area with adiabatic shear bands.
Once error estimation and distortion metrics are calculated, the next step is mesh generation. Structured and unstructured meshes can be used in the discretization of the domain. Also, in the literature three different refinement approaches have been proposed and used: h-version, in which the density of the finite elements is increased using the same interpolation order as for the elements. Moreover, the p-version, in which finite elements are fixed and the interpolation order of the elements is increased. Finally, the hp-version, which is an hybrid approach of the two approaches. Marusich and Ortiz [58] present an h refinement scheme plus continuous Delaunay triangulation at the nodes in their spatial position. Sekhon and Chenot [86] create a completely new finite element mesh based on a Delaunay-Voronoi type algorithm each time one or more elements of the mesh have got overly distorted. Peric et al. [70] present a remeshing scheme using quadrilaterals elements (constructed using Delaunay triangulation). Furthermore, the remeshing schemes presented above were applied to the numerical simulation of metal cutting in [58,86,90].
After creating a new mesh, the transfer of nodal variables (i.e. displacement, temperatures, pressure,etc.) and historydependent variables (i.e. Gauss point variables) from the old mesh to the new mesh is required. Here, the main goal of a transfer operator is the minimization of the numerical diffusion of the state variables. Ortiz and Quigley [73], proposed a transfer operator based on the weak form of the equilibrium equations in conjunction with the interpolation of nodal variables and applied it in the context of strain localization problems. That transfer operator has been applied in the numerical simulation of metal cutting by Marusich and Ortiz [58]. Lee and Bathe [52] and Peric et al. [69,70] present an adaptive mesh strategy for large deformation problems, which uses a different transfer operator for nodal variables and gauss point variables. Moreover, the nodal variables use the inverse isoparametric mapping. Where as the gauss point variables are smoothed to the nodes of the old mesh followed by transfer to the nodes of the new mesh and, finally interpolated to the Gauss points. This error estimator has been applied in the numerical simulation of high speed machining in [67,69].

Workpiece Separation Criteria
Material separation is a complex phenomenon involving many physical processes occurring at the micromechanical level. Fracture begins at the micromechanical scale and eventually is observed at macroscopic level. As a consequence, the chip separation is one of the challenges in the numerical simulation of metal cutting process. In the literature, methodologies to model chip separation can be classified as follows: continuous chip separation along a predefined cutting plane or chip separation using element deleting by killing elements based on some element deletion indicator. The last option is based on large plastic deformations and continuous remeshing. There are basically two types of indicators: those based on geometrical considerations and those based on physical considerations. Some examples of the most common indicators used in the numerical simulations with FEM are listed below: (a) a chip separation criterion based on the distance between the tool tip and the nearest node along a predefined cutting direction (based on geometry), (b) constant equivalent strain criterion and maximum shear stress criterion (based on physics), (c) fracture models like Johnson-Cook criterion and Cockroft-Latham criterion (based on physics).
The equivalent strain criterion has been a popular failure criterion for metal cutting simulations [90]. In this approach fracture is assumed to occur when plastic strain calculated at the nearest node to the cutting edge, reaches a critical value. The drawback of this method is, if uncontrolled, node separation propagates faster than the cutting speed, as a result forming a large crack ahead the tool tip. Similarly a critical stress criterion has also been suggested where node separation is activated once the material reaches a critical stress value [89].
The Johnson-Cook failure criterion is another criterion based on the postulation that the critical equivalent fracture strain is a function of stress triaxiality, strain rate and temperature. The Johnson-Cook fracture model is semiempirical in nature and necessitates the determination of constants from tensile tests with high triaxiality, shear tests and Hopkinson bar torsion tests at varying temperatures and strain rates. Johnson-Cook fracture model has been used to model machining of titanium alloy (Ti-6Al-4V) in [21]. Another similar fracture model used in machining simulations is the Cockcroft-Latham fracture criterion [20,97].
Examples of adaptive remeshing used as a strategy to separate the chip frorm the workpiece can be found in [26,58,67,86].

Chip Segmentation and Breakage
The process of segmented and discontinuous chip formation involves the propagation of fractures through the deforming chip. Chip segmentation by shear localization and breakage are important characteristics which can be observed in a certain range of cutting speeds when machining materials like titanium. Experimental studies have shown that shear localization in machining titanium alloys is due to the occurrence of thermo-plastic instability and ductile fracture. However, for structural steels depending on cutting conditions both ductile or brittle fracture can occur.
The separation criteria used to determine the propagation of the fractures is mainly based in physical models. Marusich and Ortiz [58] present a brittle fracture criteria formulated in terms of the toughness, used in conjunction with a multi-fracturing algorithm and a ductile fracture expressed in terms of the fracture strain, derived from Rice and Tracey's void growth criterion. Brittle or ductile fracture criterions are used depending on the machining conditions. Owen and Vaz [67] present a fracture criterion based on the equivalent plastic strain and the uncoupled integration of Lemaitre's damage model. These fracture criteria are able to model the material failure (thermal softening/failure softening) in problems involving adiabatic strains localization, where high speed cutting is under analysis. Boroushaki et al. [16] and Umbrello [97,98] also included damage mechanics in the simulation of crack propagation based on Lemaitre's damage model. Umbrello et al. [99] employed Brozzo's criterion to predict the effect of hydrostatic stress on the chip segmentation during orthogonal cutting. Other examples are Ceretti et al. [20] that combined the Cockroft and Latham criterion with a criterion based on the effective stress in order to optimize the material fracture in cutting operations. Furthermore, Chen et al. [21] use the Johnson-Cook fracture model to simulate serrated chip formation in titanium alloy (Ti-6Al-4V) high speed machining.

Improved Eulerian Formulations
In the previous section a review of the most important aspects related to the numerical modelling of chip formation using the FEM was presented. In this section an introduction to recent numerical techniques that have been applied for the simulation of machining process. The characteristics of each method will briefly be exposed showing the different solution approaches for the solution of machining problem. Sometimes the solution is very similar to the ones developed in the standard FEM, but other methods apply alternative solutions that does not have any relationship. Improved Eulerian formulations will be explained first. The main characteristic is these methods is that like the FEM, they are supported by a mesh. Mesh free and Particle methods will be introduced later exposing the particularities of each developed technique.

Multi Material Eulerian Method
Multi-material Eulerian (MMEM) method [10,11] is an improved Eulerian formulation based on finite elements. The Multi-material Eulerian strategy is able to deal with large deformations and free surface generation , which usually take place in the numerical simulation of machining processes. MMEM overcomes the main disadvantage of standard Eulerian Formulation . Typically, in MMEM the contact strategy between the chip and the work piece is based on mixtures theory, and this represents one of the disadvantages of the method, because mixtures theory does not predict well the friction phenomena at the tool-chip interface. Recently, Vitali et al. [101] improve MMEM using X-FEM at the interface where two materials come in contact, getting as a result a discontinuity in the velocity field at interfaces and as a consequence predicting well the friction phenomena. MMEM has been applied to the numerical simulation of AISI 4340 steel under orthogonal cutting conditions in [11], getting results that agree well with experimental results. Furthermore, the numerical simulation presented in that work used a rigid tool; friction is neglected and does not take into account heat transfer between the tool and the workpiece. One 2D example of the numerical simulation of metal cutting processes using the Multi Material Element Method is shown in Fig. 9.

Volume of Solid (VOS)
Al-Athel et al. [1] propose a new strategy to simulate orthogonal cutting based on a modified Finite Element Eulerian formulation and a Volume of Solid (VOS) approach. Finite Element Eulerian VOS formulation uses the advantages of an Eulerian formulation plus a numerical VOS scheme that can track the free surface of the material and model the unconstrained flow of the chip without being limited to only steady state scenarios with an assumed shape for the deformed chip. The VOS method allows the tracking of the free surface though (1) the calculation of the fractional volume value of the solid in each element and (2) the definition of a free surface directional vector to find the location and shape of the free surface (an interface tracking scheme).
The results from numerical simulations show good agreement in values and behavior with the ones obtained from ALE and experiments. Authors remarks that Eulerian VOS formulation is only able to predict continuous chip formation and cannot handle segmental chip formation. One advantage of the Eulerian VOS formulation is that no mesh motion scheme or remeshing strategy is required. A 2D example of the numerical simulation of metal cutting processes using the Volume of Solid method is shown in Fig. 10.

Material Point Method (MPM) or Point in Cell (PIC)
The Material Point Method (MPM) or the Point in Cell (PIC) was introduced initially in fluid dynamics by Harlow [41], extended to solid mechanics problems by Sulky [92] and applied to orthogonal cutting simulations by Wickowski [103] and Ambati et al. [2]. In the MPM, state variables are traced at a set of points (materials points) defined independently of an Eulerian mesh on which the equations of motion are formulated and solved. More in detail, the weak form is solved on a background mesh at each time step and the computed acceleration is used to update the particle positions. Later the updated particle data is used to reinstate the position and coordinates of the background mesh nodes.
As the mesh is defined in an arbitrary way, the problem of Fig. 9 Numerical modeling of metal cutting processes using Multi Material Eulerian Formulation [11] mesh distortion, which leads to difficulties in Lagrangian FEM, is avoided. Another advantage of MPM is the easy way to solve free surface problems. The main drawback of the material point method is related to the condition of stability for the procedure of time integration of dynamic equations. For a given inter particle distance equivalent to a certain mesh size, the critical time step used in MPM is many times smaller than in case of standard FEM. Furthermore, strain localization in MPM is very sensitive to the density of material points used in the calculation. Also, MPM needs more memory in comparison with standard FEM, because it is necessary to save information of the Eulerian mesh and Lagrangian particles. Ambati and co-workers [2] compared MPM and FEM in terms of the plastic strain field and temperature field, finding a good agreement between the numerical simulations. Also, the comparison shows that MPM provides a smoother chip formation and less strain localization than FEM. A 2D example of the numerical simulation of metal cutting processes using the Material Point Method is shown in Fig. 11.

Meshless and Particle Methods
The finite element method is the reference technique in the simulation of metal forming and provides excellent results with both Eulerian and Lagrangian implementations. The latter approach is more natural and direct but the large deformations involved in such processes require remeshingrezoning algorithms that increase the computational times and reduce the quality of the results (due to the mapping the state variables from the old mesh to the new mesh). Instead, the advantages of the mesh free and particle methods may be summarized as follows: 1. They can easily handle very large deformations, since the connectivity among nodes is generated as part of the computation and can change with time; 2. The method can easily handle damage of the components, such as fracture, which should prove very useful in modelings of material failure.
Due to the reasons presented above, the mesh free and particle methods overcomes the main disadvantages of FEM.

Smooth Particle Hydrodynamics (SPH)
The first element free Lagrangian technique applied to cutting process is SPH (Smoothed Particle Hydrodynamics). Limido et al. [53] apply SPH to high speed numerical cutting of an Al6061 Aluminum alloy using the commercial software LS-Dyna. Authors report that SPH results are in good agreement with the experimental observations and numerical simulations carried out with AdvantEdge. The main advantage of SPH is that it does not need a finite element mesh to calculate derivatives. Material properties and state variables are available at a set of points, called particles. This avoids severe problems associated with mesh tangling and distortion which usually occur in Finite Element Lagrangian formulations involving large deformation and strain rates. In SPH the value of a function, or its derivative can be estimated at any particle i based in the set of particles that are within a given distance h from i particle. One of the advantages of SPH is the natural workpiece-tool separation; the workpiece matter flows naturally around the tip tool. An additional advantage of SPH is contact handling, because contact is modeled as a particle interaction and friction parameters (like Coulomb friction parameter) do not have to be defined. The main disadvantage of SPH in comparison with FEM is the neighbours search, because updating the data base of neighbor particles takes usually a long time in comparison with other calculations needed during each time step. One example of the numerical simulation of metal cutting processes using SPH is shown in Fig. 12.

Finite Point Set Method (FPM)
Uhlmann et al. [96] apply the Finite Point Set Method (FPM) to model cutting of Inconel 718. FPM is an implicit scheme which is based on the differential form of the conservation laws of mass, momentum and energy. In detail, FPM is a generalized Finite Difference scheme based on Moving Least Square Method (MLS). The main advantages of FPM are: (1) Remeshing is avoided by the mesh free approach, (2) Numerical losses due to remeshing does not occur, (3) Because FPM is a lagrangian formulation allows for an easy way to represent free and dynamics boundaries. Some disadvantages of the FPM are: (1) Relatively high computational cost in comparison to FEM, due to the high number of neighbours that each node has in FPM compared to FEM. In the referenced work [96], the authors consider the tool as rigid, and no heat transfer between the tool and the piece. Moreover, authors say that FPM needs further development to simulate a more realistic chip formation process. There are some similarities between SPH and FPM, because both use a sphere of influence to study particles interaction, but SPH uses for field function approximation the kernel approximation while FPM uses Moving Least Squares. One example of the numerical simulation of metal cutting processes using Finite Point Set Method (FPM) is shown in Fig. 13.

Constrained Natural Element Method (CNEM)
The Natural Element Method (NEM) was presented and applied to problems in solid mechanics by Sukumar et al.  [53] in [91]. Illoul et al. [46] applied CNEM to 3D numerical simulation of orthogonal and oblique cutting of a Ti-6A-4V alloy. The CNEMs shape functions are built using the constrained Voronoi diagram [48]. CNEM involves three main steps: (1) Build the constrained Voronoi diagram. Thus, for each node, a Voronoi cell is geometrically defined. (2) Compute the CNEM shape functions, and (3) Calculate the internal and external forces, acceleration, velocities, and displacements using a variational formulation. In CNEM, state variables are available on particles, so there is no numerical diffusion due to an update of the Voronoi diagram. Furthermore, CNEM does not need chip-workpiece separation criteria. The main disadvantage of CNEM is the computer time needed to calculate shape functions. Another disadvantage is the necessity to remesh surfaces, where two surfaces folds, excessive elongations are produced or where nodes become too close. The numerical results presented by Illoul and Lorong have not been validated against experiments. Figure 14 shows an example of a numerical simulation of metal cutting using CNEM.

Discrete Element Method (DEM)
Cundall [25] introduced discrete element method (DEM) also called a distinct element method for the analysis of rockmechanics problems. Since then, DEM has been applied in different fields has been applied to simulate and analyze flow behavior in a wide range of disciplines including pharmaceutical and process engineering, mechanical engineering, materials science, agricultural engineering and more (see [31,43,50]). Moreover, the discrete element algorithm is a numerical technique, which solves engineering problems that are modeled as a system of distinct, interacting, and general-shaped (deformable or rigid) bodies or particles, subject to motion and deformation. Furthermore for DEM calculation cycle a time-stepping algorithm that requires the repeated application of the law of motion to each particle, a force-displacement law to each contact, and a constant updating of wall locations is used. Contacts, which may exist between two balls or between a ball and a wall, are formed and broken automatically during the simulation.
Fleissner et al. [31] applied DEM to orthogonal cutting process simulations. The authors represented the workpiece as a bulk of identical spheres arranged in a face centered cubic lattice. Particles interact by a visco-elasto-plastic law neglecting thermal effects. The numerical model of metal cutting presented in that work was not validated against experimental results. In contrast to other meshless methods, which are mainly designed to solve partial differential equations that describe the physical phenomena, DEM accounts for the simulations of particle interactions. DEM developers recommend this technology to problems which involves breakage, rupture and large deformations, and together with contact of multiple bodies. Precisely, the numerical ingredients needed for simulation of orthogonal cutting. However, at the same time, DEM developers recognize that FEM is superior to DEM for problems where small elastic strain are of interest or for the investigation on mode shape of structural oscillations. In those cases DEM needs to be further developed to predict accurate chip formation. One example of the numerical simulation of metal cutting processes using Discrete Element Method (DEM) is shown in Fig. 15.
Eberhard and Gaugele [29] show the applicability of the DEM for the modeling of an orthogonal cutting process of steel and aluminum. In that work, particles interact by a visco-elasto-plastic law including thermal effects. The numerical results are examined by comparing the simulated forces acting on the tool with experimentally obtained forces. Results presented in [29] show that using a workpiece with a regular lattice the cutting force component can be reproduced very nicely whereas the passive force shows considerable deviation. Instead, modeling a cutting process with a random lattice workpiece fails to reproduce basic qualitative characteristics of metal cutting. Furthermore, Eberhard and Gaugele [29] identify that the main challenge of DEM is to find appropriate force laws and parameters in order to synthesize the solid with correct physical properties, whereas the difficulty with FEM is found with regard to separation of material and remeshing. Eberhard and Gaugele [29] suggest that some ingredients of SPH can be added to DEM to overcome its drawbacks.

Maximum Entropy Meshfree Method
The maximum entropy meshfree method was presented in 2006 by Arroyo and Ortiz in [5]. The maximum entropy meshfree method involves three main steps: (1) calculate the Delaunay triangulation of the cloud of nodes (only used for integration). (2) calculate the maximum entropy shape functions, and (3) Calculate the internal and external forces, acceleration, velocities, and displacements using a variational formulation.
Greco et al. [35] propose a simple model of visco-plasticity, where both the pressure and velocity fields are discretized with maximum entropy approximants. The inf-sup condition is circumvented with a numerically consistent stabilized formulation that involves the gradient of the pressure. The performance of the method is studied with a 2D orthogonal cutting. To improve the quality of the results Greco et al. [35] suggest that is necessary to implement rezoning algorithms, but with respect to finite element methods that require time consuming rezoning procedures, this operation takes a negligible computational time in the Maximum entropy mesh-free approximants strategy. One example of the numerical simulation of metal cutting processes using Maximum entropy mesh-free approximants is shown in Fig. 16. Greco et al. [35]

Particle Finite Element Method (PFEM)
The Particle Finite Element Method (PFEM) is a particlebased method supported by a Lagrangian mesh. The initial developments of the Particle Finite Element Method (PFEM) took place in the field of fluid mechanics [45], because PFEM facilitates tracking and modeling of free surfaces. Later on, the Particle Finite Element (PFEM) was applied in a variety of simulation problems: fluid structure interaction with rigid bodies, erosion processes, mixing processes, coupled thermo-viscous processes and thermal diffusion problems [65,66]. First applications of PFEM to solid mechanics were done in problems involving large strains and rotations, multi body contacts and creation of new surfaces (riveting, powder filling and machining) [61]. In the PFEM, the continuum medium, considered as an infinite pack of particles each of them of infinitesimal size, is represented by a finite set (or a cloud) of infinitesimal-sized particles. The particles of the cloud describe and contain the properties of the continuum medium (displacement, pressure, temperatures, strains, stresses, internal variables, etc.) at their instantaneous locations, and, when necessary, the properties of the other particles of the continuum are obtained by interpolation from those in the cloud. Moreover, the characterization of the continuum by the PFEM has two basic ingredients: (1) continuous remeshing via Delaunay tessellations and (2) boundary recognition techniques.
A Delaunay triangulation is generated at each time step connecting the particles of the cloud on the basis of their updated positions (see [18]): this triangulation is used as the Fig. 16 Numerical modeling of metal cutting processes using Maximum entropy mesh-free approximants [35]. Where P is the pressure and V y the vertical velocity finite element mesh during the time step, and, when necessary for computational purposes, the properties of the particles (e.g. at the nodes of the mesh) are interpolated to the Gauss points. Due to the optimal properties of Delaunay triangulations for minimizing angle distortions, this continuous re-meshing minimizes the element distortions making the method suitable and advantageous, in terms of reliability and robustness, in front of more classical finite element methods.
Specifically the method becomes suitable for simulation of those industrial problems displaying material flow (cutting processes, granular material flow, metal forming processes, etc.).
Boundary recognition techniques can be naturally used in the PFEM (e.g. alpha-shape techniques). This facilitates modeling those complex mechanical processes in which new boundaries, different from the initial ones, appear.  Figure 17 shows a numerical simulation of continuous and serrated chip formation using the PFEM. The PFEM was applied to numerical simulation of metal cutting processes in 2D in Rodriguez et al. [77-80, 82, 83] and in 3D in Rodriguez et al. [81] 6 Conclusions A large number of numerical techniques have been developed and applied in order to give a response to the simulation of chip formation. It is proven that modelling machining processes still is a challenge to existing algorithms and computational tools. Large plastic deformations localized in small areas, complex contact conditions and large geometrical boundary changes are some of the major challenges associated with this class of problems.
In this paper the chip formation problem is introduced and a summary of the state-of-the-art developments in its numerical modeling is presented. As seen, significant advances have been made in developing more advanced computational tools and for fundamental modeling at the process level using the traditional 2D and 3D methods. Numerous modeling methods and techniques have emerged in recent years (particle based methods and meshless methods), with interesting application potential. Major findings of this paper are blended with statements of future research directions as follows: There is an urgent need to move from 2D to 3D model development. Most fundamental variables and machining performance measures are predicted in 2D, while industrial needs are in 3D models for specific operations. Still, the current challenges is the computational efficiency of the numerical methods in 3D. The desired numerical technique must fulfill next characteristics: should allow for the treatment of large deformations, nearly incompressibility and heat conduction, workpiece-tool contact including friction effects as well as the full thermo-mechanical coupling for contact. Furthermore, the method should be computationally inexpensive without precision loss.
New particle-based techniques seem to be the most promising lines in order to reach the ideal for the numerical modeling of the chip formation, accomplishing most of the needed characteristics in a natural way. The main goal is to make these techniques more robust and computationally cheaper.

Compliance with Ethical Standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.