Implicit or explicit time integration schemes in the PFEM modeling of metal cutting processes

This work presents the development of an explicit/implicit particle finite element method (PFEM) for the 2D modeling of metal cutting processes. The purpose is to study the efficiency of implicit and explicit time integration schemes in terms of precision, accuracy and computing time. The formulation for implicit and explicit time marching schemes is developed, and a detailed study on the explicit solution steps is presented. The PFEM remeshing procedures for insertion and removal of particles have been improved to model the multiple scales of time and/or space of the solution. The detection and treatment of the rigid tool contact are presented for both, implicit and explicit schemes. The performance of explicit/implicit integration is studied with a set of different two-dimensional orthogonal cutting tests of AISI 4340 steel at cutting speeds ranging from 1 m/s up to 30 m/s. It was shown that if the correct selection of the time integration scheme is made, the computing time can decrease up to 40 times. It allows us to affirm that the computing time of the PFEM simulations can be excessive due to the used time marching scheme independently of the meshing process. As a practical result, a set of recommendations to select the time integration schemes for a given cutting speed are given. This is intended to minimize one of the negative constraints pointed out by the industry when using metal cutting simulators.


Introduction
The description of a metal cutting process can be defined by the interaction of two physical elements. A tool cutting edge penetrates into the a metal workpiece which material is thus plastically deformed and slides of along the rake face of the cutting edge. Large plastic strains, mainly located within the primary and the secondary shear zones, produce the detachment the material and create metal chips. Machining analysis focuses on the study of chip formation, which is associated with large deformations, high strain rates and high temperatures. Even if new manufacturing techniques, such as additive manufacturing [1], have recently been developed for metal processing, machining continues to be the most used manufacturing operation in terms of volume and expenditure. This process contributes to around 5% of the GDP of industrialized countries, representing a significant part of the wealth of these economies [2,3]. The use of cutting edge technology such as numerical simulation can help to find the optimal cutting parameters and to visualize variables that are difficult to measure (e.g., the temperature and the strain rate field) with-out the need for costly experiments, thereby increasing the financial profitability of manufacturing companies.
Most of the time, the numerical modeling of metal cutting processes is based on continuum mechanics (CM) [4][5][6][7][8], where the linear momentum, energy equations and constitutive model are expressed as a set of nonlinear partial differential equations (PDE). This set of equations depends on space and time; thus, it is necessary to discretize these equations on space and time to obtain an approximate solution of a variable of interest at a time t.
The most used method for spatial discretization of the PDE for the problem of interest is based on finite element method (FEM) [9], including the FEM in their different descriptions: Lagrangian [4,6], Eulerian and arbitrary Lagrangian-Eulerian formulations [5]. Other techniques have been developed like the multi-material Eulerian method [10], the volume of solid [11], the material point method or point in cell [12], which extend the FEM by applying other spatial mesh-based discretizations, showing some advantages in comparison with the FEM. The new methods based on a fixed mesh (background mesh) eliminate the problem of mesh distortion, using elements that are filled and unoccupied with material as it moves along the fixed mesh. However, they present limited possibilities for the modeling of some cutting conditions that needs for and accurate boundaries representation (e.g., segmented chip formation). The modeling is difficult if not impossible with these techniques when a fixed mesh is used [11].
Techniques not based on CM have been developed; among them are: the discrete element method (DEM) [20,21] and molecular dynamics (MD) [22,23] with a focus on the nanometric scale. In DEM, the workpiece is represented as a bulk of identical spheres arranged at random [20] or in a facecentered cubic lattice organization [20]. The motion of DEM particles is modeled using Newton's second law, and the interparticle contact forces are modeled using springs and a damper approach. In contrast with FEM and 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. The main challenge of DEM is to find appropriate force laws and parameters in order to synthesize the solid with correct physical properties [24]. MD is a method for analyzing the physical movements of atoms and molecules. In the most common version, the trajectories of atoms and molecules are determined by numerically solving second Newton's law of motion for a system of interacting atoms/molecules, where forces between the particles and their potential energies are often calculated using interatomic potentials or molecular mechanics force fields [22,23]. Though molecular dynamics is very closely related to DEM, however, DEM is generally distinguished by considering rotational degrees of freedom for particles as well as using complex geometries for the geometrical description of particles (e.g., polyhedra) [25]. MD has been used to model the cut of pieces at nanoscale [22,23] while DEM models the cut of pieces with a size of the order of a few millimeters [24]. The main disadvantage of DEM is that the parameters of the model need to be adjusted for each case depending on the number, arrangement and size of particles [24] and for MD is the high number of particles needed to carried out a realistic numerical simulation [22,23].
Since particle-based methods and meshless methods do not depend on a mesh, they do not have a mesh distortion problem. Therefore, they have an inherent advantage when modeling problems with large deformations such as machining [9]. However, the main disadvantage of these methods comes from the needed spatial searches, for instance, the search of particles neighbors and patches (updating the data base of neighbor particles takes usually a long time in comparison with other calculations needed during each time step [9]). This can be the most time-consuming computational operation for some of these techniques when it is used with explicit time integration schemes. It means that computational time used by the integration of the governing equations changes the relevance of the auxiliary techniques used in particle methods. More information about the advantages and disadvantages of the particles and meshless methods for the numerical modeling of metal cutting processes can be found in [9].
The time integration of the differential equations has been carried out with implicit [5,6,8,17,18] or explicit [4,7,10,13,14,19] integration schemes. The selection of the scheme is usually based on its availability and recommendation of the computer software used and not based on the criterion of minimizing the computing time. Most commercial software is developed and specialized in the use of one type of time marching scheme, although some others offer different alternatives. Some top commercial software packages used in industry for the numerical simulation of metal cutting processes are: Abaqus (Implicit/Explicit) [9], Comsol (Implicit) [9], ANSYS (Implicit/Explicit) [9], AdvantEdge (Explicit) [26], Deform (Implicit) [27] and MSC. Marc (Implicit) [9]. Most of them use FEM to discretize in space the balance equations governing the physics of the problem. Abaqus also includes a SPH formulation able to model machining processes.
In the literature or in the manuals of commercial software, the following suggestions are made regarding the use of an implicit or explicit scheme: Explicit scheme must be used for applications of short duration, where high frequencies are relevant, such impact problems and explosions [28], Implicit schemes are appropriate for low velocity dynamics where the phenomena is dominated by low frequencies, such as for earthquake engineering with significant damping or for metal forming processes [28].
Following this recommendations, we can identify examples of cutting mechanics (processes including high speed cutting) where explicit schemes can be more efficient than implicit schemes. These examples define the term of High-Speed Machining (HSM), which conditions are complex and difficult fulfill by the available machine tools. It demands high precision ranges of cutting feeds and speeds, suitable tools and materials [29]. Furthermore, there is also no speed that defines the limits between low-and high-speed machining, arising the following question: At what cutting speed, are explicit schemes more efficient than implicit schemes? In the literature, there is no comparison between explicit and implicit time integration schemes of the numerical modeling of cutting process using the same spatial discretization and numerical method. Therefore, the conditions under which one scheme is better than the other have not been identified.
The main reason to develop numerical methods to simulate metal cutting processes is to observe and predict variables that are difficult to measure, with the goal to increase knowledge and understanding of the cutting process. The previous developments of the PFEM have shown its capacity to predict forces, strain, stress and to predict the transition from continuous to serrated chips with the increase in cutting speed [9]. The second reason to simulate cutting is to determine what processes parameters give the best cutting conditions, i.e., cutting speed, cutting depth and tool material. This analysis needs for expensive computing resources and large computing times, thus preventing the massive use of simulation by the industry.
The main goal of this work is to contribute to the second reason and develop criteria that allow the user of metal cutting simulations to select the most efficient time integration scheme according to their cutting conditions based in the spatial discretization. The modeling and the presented developments will be based on the Particle Finite Element Method (PFEM). The selection of time integration algorithm should not be based concerning only on computing time; it should also take into account the robustness and the convergence of the method. For example, when chip shapes or contact conditions are complex, the implicit algorithms use to present some convergence difficulties.
In order to analyze the robustness and the accuracy of the results, we will also show how to avoid the locking phenomena of linear triangle finite elements for implicit and explicit time integration by using stabilized mixed linear displacement-linear pressure finite elements. This work also shows that the excessive computing time of numerical sim-ulations of metal cutting processes is not only due to the remeshing scheme typical of the PFEM, but it is due the incorrect selection of the time integration scheme of the balance equations.
The article is organized in eight sections starting with Sect. 2, where the strong and the weak form of the balance equations is presented and the spatial and time discretization of these equations explained. Section 3 develops the transient solution of the discretized equations used in the PFEM for solid mechanics, which is briefly explained in Sect. 4. In Sect. 5, a regularized Johnson-Cook constitutive model that represents the material behavior at different temperatures, strains and strain rates is proposed. This model allows to face a large range of cutting speeds and conditions. In Sect. 6, we describe the contact phenomena at the workpiece-tool interface with a rigid cutting tool for both implicit and explicit schemes. The numerical modeling of AISI 4340 steel with implicit and explicit schemes at different cutting speeds is presented in Sect. 7. The set of results coming from a comparative analysis will allow us to determine how to select the time integration in order to minimize the computing time. The discussion about it and some recommendations are presented in Sect. 8.

Governing equations in strong form
Let us consider the body B that due to the applied loads undergoes large deformations, which occupies a region of the two-dimensional Euclidean space E with a regular boundary of its reference configuration, with material particles labeled X. A deformation of B is defined by a one-to-one mapping: that maps each particle X of the body B into a spatial point x: which represents the location of particle X of the deformed configuration of B. The region of E occupied by B in its deformed configuration is denoted as ϕ( ). The problem is governed by linear momentum and energy balance equations where ρ is the mass density, c p is the specific heat, a is the acceleration (also known as the material derivative of the velocity), θ is the temperature,θ is the time derivative of the temperature, σ is the symmetric Cauchy stress tensor, k is the thermal conductivity, b is the body force, and r is the heat source due to plasticity per unit of volume. Mathematically, the heat source due to plasticity can be expressed as where σ y +β is the equivalent stress,ε is the equivalent plastic strain rate, and η p is the fraction of the inelastic work. As usual for metal cutting simulations, the value η p is set to 0.9 [4,30].
To deal with quasi-incompressibility due to J 2 plasticity, the Cauchy stress tensor is decomposed into its deviatoric and volumetric components as where dev(σ σ σ ) is the deviatoric part of the Cauchy stress tensor, p is the pressure and I is the second order identity tensor. The pressure is assumed to be positive for a tension state. As in the present work, a Neo-Hookean material [31] with free energy given in Appendix A is used, the continuity equation can be expressed as where κ is the bulk modulus, J is the determinant of the deformation gradient F = dx/dX and the term κ ln(J ) J comes from the derivative of the volumetric component of the free energy function U (J ) of a Neo-Hookean material with respect to J . The reason to introduce the decomposition of the Cauchy stress in its deviatoric and hydrostatic components in Eq. (5) and the introduction of the continuity equation, Eq. (6), is because in this work we use a mixed formulation with linear triangle finite elements able to handle the incompressibility constraint due to J 2 plasticity. The details about the finite element discretization are presented in Sect. 2.3 .
The balance equations are solved numerically for a twodimensional region ⊆ R 2 , for the time range t ∈ [0, T ], given the following boundary conditions on the Dirichlet (ϕ( u ) and ϕ( θ )) and Neumann boundaries (ϕ( t ) and ϕ( q )), respectively u =û on u (7a) whereû,t,θ andq are the prescribed displacements, tractions, temperatures and heat fluxes, respectively, and n is the unit outward normal to the boundaries of the solid body (the workpiece). In Eq. (7), the boundary = u ∪ t = θ ∪ q and u ∩ t = θ ∩ q = ∅. Additionally, we assume that the following initial data is specified for the mechanical and thermal fields where u 0 (x),u 0 (x) and θ 0 (x) are the initial displacement, velocity and temperature of the workpiece, respectively. A constitutive equation for the evaluation of the stressstrain relation that includes the strain hardening, the strain rate hardening and the thermal softening is also needed to fully define the boundary value problem (BVP).

The weak form of the governing equations
Following the standard FEM procedure, the weak form of the momentum balance and energy equation is obtained in the following way: the L 2 inner product of Eqs. (3) is derived using arbitrary test functions w and ζ , such that where T is the space of virtual temperatures. By using the divergence theorem, the weak forms of momentum balance and energy equation can be obtained and expressed as where ∇ S w represents the symmetric part of the tensor ∇w, dv denotes the infinitesimal volume for the current configuration and da denotes the infinitesimal area for the current configuration. The weak form of the pressure constitutive equation can be obtained and expressed as where q is the test function, such that q = {q ∈ Q|q = 0 on ( u )}, where Q is the space of virtual pressures. For the treatment of the incompressibility constraint, the Polynomial Pressure Projection (PPP), introduced by Dohrmann and Bochev [32,33] is used. The PPP modifies the weak form of the pressure constitutive equation as where α is a stabilization parameter (α is constant and equal to 1 for all the numerical simulations carried out in this work), G the shear modulus,p the best approximation ofp in (Q 0 ) andq ∈ Q 0 an arbitrary test function, where Q 0 is the space of polynomial functions with zero degree in each coordinate direction.
The best approximationp of p of Q 0 is obtained through the following integral The previous integral represents the L 2 projection of the pressure with linear shape function in a discontinuous space (zero degree polynomial).
The main advantage of PPP compared to other stabilization techniques is that the pressure stabilization is accomplished without the use of the residual of the momentum equation; thus, the calculation of higher-order derivatives and the specification of a mesh size dependent stabilization parameter are avoided. Furthermore, the addition of the second term in Eq. (11) in a FEM context is performed at the elementary level.

Spatial discretization using FEM
Following a standard Galerkin approach, the displacements, the pressure, the temperature and the test functions are interpolated with linear shape functions N: where the upper symbol (·) denotes a vector containing the nodal variables. The same order interpolation functions are used for scalar and vector fields: N and N u .
The matrix form of the Galerkin expression of Eq. (9) is obtained introducing the previous approximations to the weak form, The explicit form of each matrix of the system above can be written as where B u is the small deformation strain displacement matrix, B θ is gradient temperature matrix, respectively, andσ σ σ corresponds to the Voigt notation of the Cauchy stress tensor. The discrete counterpart of the stabilized pressure constitutive equation is given by As N contain the set of polynomials of order 1,Ñ contain the set of polynomials of order 0. The above global force vectors are obtained as the assemblies of element vectors as it is standard of finite elements theory. Given a nodal point, each component of the global force associated with a particular global node is obtained as the sum of the corresponding contributions from the element force vectors of all elements that share the node. In this work, the element force vectors are evaluated using Gaussian quadratures.

Time-integration schemes for the mechanical problem
We present subsequently the main steps for solving nonlinear dynamics problem presented in Eq. (14) (a) by using Newmark time integration schemes. To that end, the time interval of interest [0, T ] is subdivided in the chosen number of time steps, which will specify the time instants where the selected time integration scheme should deliver the solution: the solution at time t n is defined as: The time integration algorithm implemented in this work belongs to the family of Newmark methods. The method seeks the solution of linear momentum at time t n+1 starting with the solution a time t n in a the step-by-step basis.
Equation (14) (a) is re-written as The problem consists of finding a displacement function that satisfies Eq. (27) and the corresponding initial conditions; the set of equations for the Newmark method defines the displacement and velocity fields at t n+1 as an interpolation of the displacements and velocities a t n , and acceleration field at t n and t n+1 , which for the displacement fields reads: and velocitȳ Hence, by choosing β = 0 and γ = 1/2 we obtain the central difference scheme (called explicit scheme in this work), whereas for β = 1/4 and γ = 1/2 the general expression will reduce to the trapezoidal rule approximation (called implicit scheme in this work). In Eqs. (28) and (29), h represents time increment between t n and t n+1 , defined as: The time increment is also know as the time step.

Trapezoidal rule with displacement formulation (Implicit)
From Eqs. (28) and (29), and using the trapezoidal rule, the velocityv n+1 and accelerationā n+1 at t n+1 reads as The last two equations are implicit in the sense that they depend upon the displacement value at t n+1 . As Newton-Raphson method is employed for the iterative solution of the set of nonlinear algebraic equations with displacements as main unknowns, the consistent linearization of momentum equation is necessary, which can be written as: where the subscript (i) refers to the iteration number of the Newton-Raphson method and K n+1 As is usual for nonlinear FEM, the tangent matrix is evaluated as the sum of the material stiffness matrix and the geometric stiffness matrix (see [31,34]).
Equation (32) is completed with the update of displacement according tō where ū n+1 (i) is the increment of displacement during one Newton-Raphson iteration.

Trapezoidal rule with mixed formulation for the mechanical problem (Implicit)
The mixed formulation is based on linear interpolation finite elements for both displacement and pressure. In case of a mixed formulation with trapezoidal rule Eqs. (14) (a) and (21) should be solved in a coupled way. The time discretization of the previous equations is written as whereā n+1 is given by Eq. (31).
The consistent linearization of Eqs. (35)(a) and (35) where m represent the second order identity tensor in Voigt notation. It is important to remark that mix K n+1 uu,(i) = K n+1 uu,(i) , the mathematical expressions for both matrices are given for Appendix B. Details about how to obtain the matrices are given in nonlinear finite element books like [31,34]. Like K n+1 uu,(i) , mix K n+1 uu,(i) is the sum of the material stiffness matrix and the geometric stiffness matrix. Furthermore, in the mixed formulation the Cauchy stress in Eq. (41) is evaluated with Eq. (5).
Equation (36) is completed with the update of displacement and pressure according tō where ū n+1 (i) and p n+1 (i) are the increment of displacement and pressure during one Newton-Raphson iteration, respectively.

Central difference scheme with mixed formulation (Explicit)
In the central difference scheme, the displacement and velocity are updated at t n+1 as The last two equations are explicit in the sense that they depend upon the displacement, velocity and acceleration at t n . In Eq. (45), the acceleration at t n+1 is obtained from the momentum equation evaluated in a special way as follows To decrease computational cost for the central difference scheme, M n+1 u is approximated by its diagonal form (lumped mass matrix).
The pressure at t n+1 is obtained with the implicit solution of the pressure from Eq. (21) in the following way The consecutive solution of Eqs. (44), (45), (46) and (47) represents a typical solution of a time step from time t n to time t n+1 using the central difference scheme with mixed formulation. One way to avoid the implicit calculation of the pressure field (Eq. (47)) is the use of the displacement/irreducible formulation but with the locking problems due to incompressibility as it is shown for some examples presented in Sect. 7.
A formulation similar to the one presented in this section and stabilized using finite calculus has been presented in [35]. The previous formulation has the advantage that the pressure calculation is done explicitly, as a result a decrease of the computing time is expected. Another fully explicit formulation that solves the problem of locking of linear triangle finite elements that does not need the calculation of a implicit pressure have been developed; among them are F-bar [36] and average nodal strain/stresses [37]. The exploration of explicit formulations for both displacement and pressure that uses linear triangles finite elements will allow to find new finite element formulations that obtain the same results but with a lower computing cost. The study of the previous formulations as well as their computational efficiency will be investigated in future works.

Time-integration schemes of the thermal problem
To solve the thermal problem, it is necessary to approach the time derivative of the temperatureθ . The forward Euler scheme approximates the time derivative of θ at time t n aṡ meanwhile, the backward Euler scheme approximates the time derivative of θ at time t n+1 aṡ Solving for θ n+1 in Eq. (48), we obtain θ n+1 = θ n +θ n h; (50) the last equation is an explicit update because it depends on the values at time t n . Solving for θ n+1 in Eq. (49), we find the last equation is an implicit update because it depends on the values at time t n+1 . The time discretization of heat equation (Eq. (14)) is given by: For the solution of Eq. (52) using Backward Euler scheme, the consistent linearization is necessary; the linearization is given by: where the subscript i refers to the iteration number of the Newton-Raphson method and where Equation (56) is completed with the update of temperature according tō For the Forward Euler time integration scheme, we evaluate the heat equation (Eq. (14)) at time t n as Solving forθ n , we obtaiṅ After gettingθ n , the temperature at t n+1 is calculated using Eq. (51). It is important to remark that for computational efficiency the lumped mass matrix approximation will be used for the calculations.

Transient solution of the discretized equations
The thermo-mechanical problem will be solved using a staggered approach with a mechanical phase with the temperature held constant, followed by a thermal phase at a fixed configuration. The isothermal split solves the mechanical problem with a predicted value of temperature equal to the temperature of the last converged time step t n and, then, solves the thermal problem using the configuration obtained as a solution of the mechanical problem at t n+1 . In Sects. 3.1 and 3.2, the transient solution of the thermomechanical problem using the trapezoidal rule with mixed formulation (implicit) and central difference scheme with mixed formulation (explicit) is presented, respectively.

Implicit schemes
Equations (35) and (52) are solved in time with two implicit Newton-Raphson type iterative schemes, one for the mechanical problem and one for the thermal problem [28,31,34]. The basic steps within a time increment [n, n + 1] are: Initialize variables: =b e,n ,ē =ē p,n whereb e, is volumetric preserving part of the elastic left Cauchy-Green tensor,ē p is the equivalent plastic strain and x n (i) is the position of the particles a time t n and iteration i.
• Iteration loop for the mechanical problem: i = 1, . . . , N iter,mech , where N iter,mech is the maximum number of iterations.
1. Compute the nodal displacement and pressure increments (It is important to remark that the term R u,(i) includes the contribution of the contact forces as it is shown in Sect. 6) The solution of the linear systems is carried out using a LU factorization. 2. Update the nodal displacements and pressures where γ is a parameter between 0 and 1, γ is given by the line search strategy. The line search strategy implemented in this work is presented in [31]. 3. Update the nodal coordinates, velocities and accelerations 4. Calculate the relative deformation gradient measured between the current configuration x n and the updated configuration x n+1 : F = dx n+1 /dx n , the Jacobian J = det(F) and the area of all elements of the mesh 5. Calculate the Cauchy Stress σ σ σ ,b e andē p at n+1 using an implicit stress update presented in Chapter 9 of [38]. Also, calculate the tangent constitutive matrix. 6. Calculate the residuals R u,(i) and R p,(i) 7. Check convergence of the mechanical problem where TOL u and TOL p are prescribed error norms.
In the set of examples presented in the paper we have set TOL u = TOL p = 10 −4 .
If conditions (63) and (64) are satisfied proceed to the thermal problem. Otherwise, make the iteration counter i ← i + 1 and repeat Steps 1-7. • Iteration loop for the thermal problem: i = 1, . . . , N iter,therm , where N iter,therm is the maximum number of iterations.
1. Compute the nodal temperature increments using a LU factorization of the linear system.
2. Update the nodal temperatures 3. Calculateb e ,ē p , the heat source due to plasticity r and the derivative of r with respect to temperature 4. Calculate the residual R n+1 θ,(i) 5. Check convergence of the thermal problem where TOL θ is the prescribed error norms of temperatures. In the set of examples presented in the paper, we have set TOL d θ = 10 −4 .
If conditions (67) and (68) are satisfied, then make n ← n + 1 and proceed to the next time step. Otherwise, make the iteration counter i ← i + 1 and repeat Steps 1-5.

Explicit schemes
Equations (35) and (52) are solved in time with two explicit schemes, one for the mechanical problem and one for the thermal problem [28,31,34]. The basic steps within a time increment [n, n + 1] are: Given: ū n ,θ n ,x n ,v n ,b e,n ,ē p,n • Loop for the mechanical problem.

Calculate the time step
where l is the minimum distance between two particles of the mesh, 0 < C ≤ 1 is a constant (0.8 used for all the numerical simulations), ρ is the density and E is the elastic modulus of the material. As explicit scheme are only conditionally stable, the maximum size of the time step allowed is given by the well known Courant-Friedrichs-Lewy (CFL) criterion (Eq. (69)).

Update the nodal displacements and positions
3. Calculate the relative deformation gradient measured between the current configuration x n and the updated configuration x n+1 : F = dx n+1 /dx n , the Jacobian J = det(F) and the area of all elements of the mesh 4. Calculate the Cauchy Stress σ σ σ ,b e andē p at n + 1 using a implicit stress update presented in Chapter 9 of [38].

Calculate nodal accelerations
where M n+1 u is the lumped mass matrix. 6. Update the particle/nodal velocities 7. Correct the nodal/particle positions, displacements, velocities and accelerations of the nodes that are in contact with the tool. 8. Calculate implicitly the nodal/particle pressures.
The linear system is solved using LU factorization. • Loop for the thermal problem.
1. Calculate the heat source due to plasticity r . 2. Calculate the nodal temperatures.
where M n+1 θ is the lumped mass matrix.
Make n ← n + 1 and proceed to the next time step. The step 8 of the loop for the mechanical problem is not necessary for the irreducible formulation.

The particle finite element method
The particle finite element method (PFEM) is a numerical method based on the followings ingredients: (1) an updated Lagrangian formulation, (2) the finite element method with displacement/velocity and pressures as main variables, (3) the Delaunay triangulation and (4) the alpha shape method [39,40]. The combination of the four ingredients has allowed to model with great success fully and quasi incompressible free-surface flows using a Lagrangian framework, solving the mesh distortion limitation of the use of Lagrangian FEM formulations. As PFEM uses a Lagrangian formulation no convective terms and their inherent numerical difficulties are present in the balance equations. Since its original development, the application of PFEM has spread to problems in solid mechanics [41], thermo-mechanics [9,42], granular materials [43][44][45][46][47], fluid-structure interaction [43,48], contact of deformable bodies [49] and additive manufacturing processes [50]. Both explicit and implicit time integration schemes can be used for PFEM, although most of the times an implicit time integration has been preferred [51]. The PFEM has been used with an explicit time integration schemes, for example, in [52][53][54] while in [55][56][57] implicit schemes were proposed. The boundary of the set of particles for the original PFEM is identified with the alpha shape method [58]; in this work the boundary is preserved using a constrained Delaunay triangulation to guarantee mass conservation and to avoid nonphysical increase in the flexural strength of the chip (see Fig. 1 where the purple elements represent the unphysical mass added by the alpha shape). In this work, the internal variables that describe the plastic phenomena are projected from old integration points to the new integration points using a minimum distance criteria to minimize the diffusion of this variables of the mapping from mesh to mesh (see Fig. 2). This is another difference with respect to the original PFEM. The previous ingredients improves the performance of the PFEM when it is applied to the numerical modeling of metal cutting processes as it is shown in [8].
The basics steps of the PFEM used in this work are: Given a set of particles with a boundary that represent the workpiece at the beginning of the simulation Set tool incremental displacement = 0.
1. Generate a constrained Delaunay triangulation using the particles as nodes. 2. Explicitly or implicitly solve the Lagrangian form of the mechanical problem. 3. Explicitly or implicitly solve the thermal problem. 4. Update particle position, velocities, accelerations. pressures and temperatures. 5. Set tool incremental displacement = tool incremental displacement + time step tool displacement. 6. If the tool incremental displacement is greater than a critical value (0.002 mm for all the numerical simulations) remesh the workpiece with the insertion of particles at the Gauss point of the mesh whose plastic power is greater than a given tolerance or remove particles that are closer than a critical distance, otherwise go to step 2. Set tool incremental displacement = 0.
The minimum distance between particles is based on the following heuristic criteria: where h min is the minimum distance between particles, d is the cutting depth, and r is the tool radius. The previous criteria was designed to avoid the use of very small mesh size when tool radius is close to 0 (sharp tool).

Johnson-Cook strength model
The constitutive model used for the modeling of metal cutting process should describe the behavior of the material at low and high strain rates. The Johnson-Cook (JC) material model describes the behavior of the metal in a wide ranges of temperatures, strains and strain rates. The JC strength model [59] is a phenomenological model that describes the flow of the material as the product of three mathematical terms with A being the initial yield strength of the material at room temperature and a reference strain rate.
Here, ε is the equivalent plastic strain,ε is the strain rate non-dimensionalized by the reference strain rate, σ y denotes the flow stress, β the isotropic nonlinear hardening modulus, and B, C, m and n are fitting constants. In Eq. (77), θ is given by where θ 0 is the room temperature, θ melt is the melt temperature and θ is the workpiece temperature. The nondimensionalized plastic strain rateε is given by: e p is the effective plastic strain rate,ε 0 is the reference strain rate.
The JC model assumes an independent effect of the strains, strain rates and the temperatures. The expression for the first bracket of Eq. (77) represents the nonlinear hardening law, the expressions of the second terms model the strain rate hardening, and third bracket represents thermal softening, respectively. Table 1 presents the JC parameters for the AISI 4340 steel used for the numerical simulations.
The non-dimensionalized strain rateε of the JC is modified to increase the robustness of the numerical method in the following way: with the previous modification the problem of the convergence of the Newton-Raphson algorithm due to the change of lnε from positive to negative whenε is close to 1 is avoided with a minimum effect on the stress-strain curve. Details about the implementation of an implicit backward Euler time integration scheme of the JC model are presented in [60,61] and also included as an Annex A for completeness. For the present study, only a minor modification to the time integration presented in [60,61] is necessary to include the new non-dimensionalized strain rate.

Contact modeling
In this work, the cutting tool moves along the horizontal direction from right to left with a specific cutting speed and cutting depth. The tool is rigid and represented by three curves: one line that represents the rake face, one line that represents the flank face and a circle that represents the tool tip. So the tool is defined by the rake angle θ r , the flank angle θ f , the tool radius R and the center of the tool tip (x c , y c ) as shown in Fig. 3.
To identify which surface a particle is in contact with, the tool is divided in 3 regions, as it is shown in Fig. 3. Mathematically, the three region can be represented by the following inequalities: and For the implicit time integration scheme, a penalty-based approach is used. That approach allows the penetration of workpiece particles inside the tool surface. When a penetration is detected, a force perpendicular to the closest tool surface is added to Eq. (35)(a) as an external force. Mathematically, the contact force using a penalty formulation is represented as where κ is the penalty parameter, g N is the normal gap and n is the normal to the contact surface. The penalty parameter used here is: where m is the mass of the node that is in contact and h is the time interval.
In the case of the explicit time integration scheme, a predictor-corrector strategy (PCS) is used. The PCS allows the penetration of the workpiece particles inside the tool. When the predicted position Eq. (70) of boundary particle penetrates the tool, an acceleration along the normal direction of the closest tool surface is imposed to that particle such that the particle stays over the tool surface for the beginning of the next time step [4,62], eliminating the unwanted penetration. The PCS was preferred over a penalty approach to avoid that the contact penalty factor puts a narrower restriction on the time step [63].
The normal acceleration correction can be expressed as: where j represents one of the nodes that penetrates the tool, h is the time interval and n j is the normal vector to the contact surface related to node j. The previous acceleration is added to the acceleration obtained in Eq. (72) for the Step 7 of Sect. 3.2.

Examples
The PFEM with plane strain assumption is applied to the modeling of continuous chip formation process of orthogonal cutting of a rectangular block with cutting depth of 0.1 mm with both implicit and explicit time integration schemes of the balance equations. The workpiece is represented with a rectangle of height 0.5 mm and length 2 mm (see Fig. 4). The strain hardening, strain rate hardening and thermal softening of the AISI-4340 steel are described with the Johnson-Cook flow stress model. The material properties of the flow model   Tables 1 and  2, respectively. The cutting tool is assumed to be rigid and isothermal with a rake angle of 0 • , a flank angle of 5 • and a tool radius of 0.025 mm. There is no friction for the contact zone between the workpiece and the chip. The melting temperature and the temperature at the beginning of the simulation are set to 1648 K and 293 K, respectively. The cutting speed varies between 1 and 30 m/s. The results obtained will be analyzed and compared in the following sections for terms of computing time, chip shapes, stresses, temperature and cutting forces. Also, a comparison with the displacement/irreducible formulation as a possible form to avoid the implicit calculation of pressures is presented at the end of this section.

Computing time
Traditionally, most of the research of numerical modeling of metal cutting processes has been carried out with implicit or explicit commercial FEM codes. In some works, a comparison of the numerical simulations of cutting using different commercial simulation codes with implicit and explicit time integration schemes has been carried out in [64]. However, the comparison of the efficiency of the time integration schemes for the previous works is very difficult or has not been carried out, because the FEM formulation, the finite element used, the programming language of the simulation code and the computer where the simulations were run are very different. For example, some simulation programs have a flow formulation (negligible elastic strains) and others an elastoplastic formulation. The comparison of the two previous formulations does not make much sense as the flow formulation cannot predict residual stresses while the elastoplastic can predict residual stresses. All the numerical simulations presented in this work were run on a Dell computer with an Intel Xeon(R) CPU E5-1660 v3 @ 3.00GHz x 16 processor and 32 GB of RAM (Located in Luleå, Norrbotten. Sweden). The PFEM was programmed using vectorized MATLAB programming for both explicit and implicit schemes. The routines used for explicit and implicit codes used for the calculation of the stresses, mechanical internal forces, mechanical external forces, thermal internal forces, thermal external forces, the inertial forces and the remeshing are exactly the same (see Sect. 3). The differences arise for the update of nodal position and temperatures as shown in Sect. 3.
A total of 34 simulations were performed, 17 with the implicit scheme and the other 17 with the explicit scheme (see Tables 4 and 5 ). In all the simulations, the number of particles varies between 787 at the beginning of the simulations and around 5706 at the end of the simulation, corresponding to a minimum distance between particles of 0.004 mm.
The remeshing scheme ensures that the element size is never smaller that 0.004 mm, avoiding that a very small finite element can jeopardize the efficiency of the entire code for explicit simulation. All the implicit simulations need 2500 time steps for a tool advance of 1 mm, such that for a time step the tool advance a distance of 1/2500 mm. The implicit simulations were carried out using 2500 time steps because with that number of time steps corresponds a minimum computing time as it is shown in Table 3 and also corresponds to the minimum number of time steps with which the iterative Newton-Raphson scheme converges. Each mechanical implicit time step needs in average three Newton-Raphson iterations to reach a convergence of forces of 1 × 10 −4 , and each thermal implicit time step needs two iterations to reach the same tolerance. The remeshing frequency of the PFEM takes place each 5 implicit time steps where particles are inserted and removed; the frequency of remeshing corresponds to a tool displacement of 0.002 mm. The average computing time of the implicit simulations is 0.89 h, being more or less the same for all cutting speeds (see Table 4). In explicit simulations, the remeshings were carried out after the tool had advanced a distance of 0.002 mm, such that a given simulation implies the same number of remeshings for both implicit and explicit schemes (each of the simulations presented in Table 4 Table 4). Table 5 shows that for velocities greater than around 27 m/s the explicit scheme outperforms the implicit scheme. The results in Tables 4 and 5 show that the correct selection of the implicit or explicit time integration scheme can generate savings of computing of up to 40 times. One way to make the explicit scheme competitive at low speeds is the use of mass scaling, although it should be used with care because it can change the results completely due to the significant increase in the chip inertia. For that reason, mass scaling was not used in this work. Table 5 presents a comparison at higher cuttings speeds but using a tool displacement of 0.5 mm keeping the other parameters equal to the ones used for the simulations presented in Table 4. The results presented in Table 5 verify that for velocities greater than 27 m/s the explicit schemes outperform the implicit schemes. The reason why higher cuttings speeds were run with a tool displacement of 0.5 mm is to avoid the contact of the chip with the top surface of the piece (self-contact) that take place if the simulations at higher cuttings speeds are run using a tool displacement of 1 mm.
In the following section, we present a comparison of the predicted chip shapes at different cuttings speeds with both implicit and explicit schemes. Figure 5 shows the predicted chip shapes at different cutting speeds using the explicit and implicit scheme, respectively. The curvature of the chip is larger at a higher speed, resulting in the chip coming into contact with the workpiece at a smaller tool displacement. This is due to the fact that at higher cutting speeds the effect of heat diffusion is lower, generating locally higher temperatures and a significant thermal softening of the workpiece material.

Predicted chip shapes
The comparison in Fig. 5 shows that the tendency of predicted chip shape of both schemes seems to be the same. A detailed comparison of the chip shapes is presented in Fig. 6 for cutting speeds of 1 m/s and 10 m/s, respectively. Figure  6a shows that at 1 m/s there is a small difference between the predicted deformed chip thickness for close proximity to the primary shear zone. The other difference occurs at the highest point reached by the chip (see Fig. 6a). A comparison of the chip shape made with a smaller mesh size shows that these small differences become smaller. When comparing the shape of the chip predicted at 10 m/s (see Fig. 6b), differences are observed in the neighborhood of the highest point of both internal and external surfaces of the chip. In the same way as for the 1 m/s case, that small difference of the shape of the chip disappears with a small decrease in the size of the mesh (distance between particles). The distance between the particles (mesh size) for all the numerical simulations was the same as the distance used in Fig. 6 because, with that mesh size, the convergence of the cutting and feed forces is obtained as it is shown in Sect. 7.4.  Fig. 7) and at 10 m/s is around 1344 MPa (see Fig. 10). The difference of the maximum von Misses stress between the numerical schemes is 4 MPa at 1 m/s and 2 MPa at 10 m/s (the same scales and ranges are used in Figs. 7 and 10 to simplify the comparison). These results show the similarities of the prediction between numerical schemes and the thermal softening arises at higher cutting speeds (smaller value of maximum von Mises stress at high cutting speeds).

Predicted fields
In Figs. 8 and 11, the predicted contour field of temperature at 1 m/s and 10 m/s using implicit and explicit schemes is shown. At 1 m/s, the maximum value of the temperature is around 686 K and at 10 m/s is around 920 K. As expected, the highest temperature occurs at 10 m/s, due to the shorter time allowed for heat conduction take place, during the generated local heating. The difference between the maximum temperature predicted with implicit and explicit schemes at 1 m/s is 3 K and the difference at 10 m/s is around 30 K. The last measurement is not an indication that the two integration schemes give different results because with a small decrease in the mesh size the temperature differences between the schemes disappear. In addition, the similarity between the results obtained by both schemes can be verified by comparing global variables such as the force as it is presented in Sect. 7.4.
The strain rate field is shown in Figs. 9 and 12 for 1 m/s and 10 m/s, respectively. The strain rate reaches a value of around 5×10 4 and 7×10 5 1/s for cutting speeds of 1 m/s and 10 m/s, showing higher strain rates at higher cutting speeds. According to Figs. 9 and 12, the predictions of the strain rate field for the implicit and explicit schemes are very similar, only showing a small difference near the tool tip.

Predicted forces
In Figs. 13, 16, 14 and 17, the predicted cutting and feed forces are shown for cutting speeds of 1 and 10 m/s for both implicit and explicit schemes. In Figs. 13, 16, 14 and 17, both the implicit and explicit schemes tend to the same value for both cuttings speeds; the only difference between the schemes is the noise which is typical of explicit integration schemes (No damping was used to attenuate noise for explicit simulations). In Fig. 15, a comparison of the prediction of the cutting force at 10 m/s with both explicit and implicit scheme is presented; a low-pass filter was used for the explicit signal to remove the nonphysical noise as it is suggested in [65]. The previous figure shows clearly that both schemes predict exactly the same cutting force. Table 6 shows the predicted cutting and feed forces for the cutting speeds shown in Table 4 using implicit and explicit schemes. The value presented in Table 6 corresponds to the average value at the steady state. Table 6 shows that the cutting force decreases slightly with increasing speed and then remains constant. Meanwhile, the same table shows that the feed force increases with increasing cutting speed up to a speed of 13 m/s and then remains almost constant. The difference between the predictions with implicit and explicit simulations is always less than 3.5 N and in terms of percentage is less than 3% as it is shown in Table 6 . This comparison shows that the two schemes yields equivalent results. However, the advantage of implicit scheme over the explicit is clearly seen in terms of computational time at velocities smaller than 27 m/s as is shown in Sect. 7.1. At

Locking with the irreducible formulation
In the following section, results with the irreducible/ displacement based formulation are presented. The goal of an irreducible formulation is to avoid the implicit calculation of the pressure field of the explicit scheme developed in Sect. 2.4.3. and as a result decrease the computing time of the explicit scheme.
In Figs. 18 and 19, a comparison of the predicted cutting and feed forces using the mixed and the irreducible/displacement formulation is shown with an explicit time integration scheme. The irreducible formulation overpredicts the cutting force with around 10 N (6% higher); meanwhile, the feed forces are correctly predicted. The previous effect is the well-known volumetric locking problem due to plasticity of a finite element formulation with linear triangle finite elements. The previous results show the importance of the correct selection of the finite element when modeling metal cutting, since a wrong selection of the element can generate large prediction errors of important variables for cutting such as forces. If an implicit formulation with the irreducible formulation is used, an over-prediction in the cutting of around 6% similar to the ones obtained with the explicit one is obtained.   Fig. 21 show that the mixed formulation is free of pressure oscillations; meanwhile, in Fig. 20 the irreducible formulation shows the typical nonphysical spurious oscillations of the pressure field.
The motivation of this section is to quantity the error of the prediction of variables like cutting forces, pressure among others made when an irreducible formulation is used as some researchers [67][68][69][70] continue to use linear triangle finite elements to model metal cutting processes knowing that this finite element is prone to locking phenomena due to plasticity.

Discussion and conclusions
Metal cutting or machining is a process in which a thin layer of metal, the chip, is removed by a wedge-shaped tool from a large body. Cutting is a complex physical phenomena where adiabatic shear bands, excessive heating, large strains and high strain rates are present. Cutting speed plays an important role for chip morphology, cutting forces, energy consumption and tool wear. Due to the complex physical process of cutting, machining processes is one of the most challenging indus-  trial problems to be analyzed from the numerical simulation point of view. To optimize and predict industrial machining an accurate, robust and efficient numerical simulation tool is needed. Here implicit and explicit numerical integration schemes for machining simulations are investigated with the goal to give insight of the best choices for accuracy and efficiency. Two time integration schemes of the balance equations were implemented within the PFEM framework and applied to the numerical modeling of metal cutting processes. Each of the schemes was tested at 17 cutting speeds, such that a total of 34 numerical simulations were carried out. The results presented here show that both schemes predicted the same value of stresses, temperature, strain rates and forces. The difference between the numerical schemes is evident when a comparison of the computing times is carried out. The results show that a selection of the implicit time integration scheme saves 40 times regarding computing time at a cutting speed of 1 m/s (low speed machining). Meanwhile, at 30 m/s the computing time with the explicit scheme is 75% of the time needed for the implicit scheme. As such, the excessive computing time attributed to the numerical simulation of cutting is not only attributed to the remeshing, but also to an incorrect selection of the time integration scheme. This results can be used as a guide that allows people who simulate cutting to select the time integration scheme that results in a minimum computing time. For example, for cutting under similar conditions to the ones presented in this work, the implicit scheme is recommend for cutting speeds  At the same time, the results show how harmful it is to use a displacement/irreducible formulation with linear triangular finite elements for the numerical simulations of metal cutting in terms of the over-prediction of cutting forces and the nonphysical spurious oscillations of the pressure field. For example, the spurious oscillations of the pressure field of the irreducible formulation would not allow the use of a failure criteria that depends on the pressure-deviatoric stress ratio. The previous results allow to quantify the error that some researchers [67][68][69][70] are committing when they use linear triangle finite elements to model metal cutting pro-cesses knowing that this finite element is prone to locking phenomena due to plasticity. The problem is that for some regions of the chip the irreducible formulation will predict compression pressures while in reality they are in tension thus predicting no damage where damage should develop. The previous situation will invalidate the application of the irreducible formulation with linear triangle finite elements to the numerical modeling of segmented chip, where fracture criteria are of paramount importance. The stabilized mixed formulation with linear displacement and linear pressure for both the explicit and implicit time integration solves the problems of the irreducible formulation mentioned above. The previous results show the relevance of the selection of the  The results presented here show that the PFEM solves the problems of mesh distortion and diffusion, typical of standard FEA, placing the PFEM in a competitive position with the FEM. A literature overview allows us to conclude that, among the new numerical methods developed for orthogonal cutting simulations, the only one that has been tested to model cutting with implicit and explicit time integration schemes is the PFEM. The other numerical methods have only implement one of the time integration schemes [30,71]. The development of a code with both implicit and explicit integration schemes will allow the activation of an explicit scheme when an implicit scheme is not able to converge, increasing in this way the robustness of the numerical code.
In summary, the research carried out in this work allows to conclude that the implicit scheme is recommend to be used   at velocities lower than 27 m/s for metal cutting problems with similar mesh size, number of finite elements, contact parameters and material model to the ones used in this work because its computational efficiency is better compared to the explicit scheme; meanwhile, the explicit is the option to choose for velocities higher than 27 m/s.