Modeling of the fracture energy on the finite element simulation in Ti6Al4V alloy machining

One of the main problems that exists when working with Finite Element Methods (FEM) applied to machining processes is the lack of adequate experimental data for simulating the material properties. Moreover, for damage models based on fracture energy, the correct selection of the energy value is critical for the chip formation process. It is usually difficult to obtain the fracture energy values and requires complex tests. In this work, an analysis of the influence of this fracture energy on the cutting force and the chip generation process has been carried out for different sets of cutting parameters. The aim is to present an empirical relationship, that allows selecting the fracture energy based on the cutting force and cutting parameters. The work is based on a FEM model of an orthogonal turning process for Ti6Al4V alloy using Abaqus/Explicit and the fracture energy empirical relation. This work shows that it is necessary to adjust the fracture energy for each combination of cutting conditions, to be able to fit the experimental results. The cutting force and the chip geometry are analyzed, showing how the developed model adapts to the experimental results. It shows that as the cutting speed and the feed increase, the fracture energy value that best adapts to the model decreases. The evolution shows a more pronounced decrease related to the feed increment and high cutting speed.

for the friction coefficient at the tool-chip-workpiece interface, verifying the results with experimental analysis. Childs et al. 13 combined FEM simulation with experimental tests to develop a constitutive equation. This combination better describes the flow stress,failure behavior and also examines which parameters are critical to obtain good results. Sadeghifar et al. 14 , implemented the Johnson-Mehl-Avrami-Kolmogorov (JMAK) recrystallization model. This model is used to study the effect of grain size in Ti6Al4V machining using different cutting parameters. They studied how the surface integrity is linked to the fatigue life and corrosion resistance of machined parts. Bai et al. 10 have also proposed another analytical model for chip formation prediction in orthogonal cutting of Ti6Al4V.
Even though there are several studies on Ti6Al4V, the material parameters are, in most cases (mainly the flow and failure data), provided by different sources which are assumed to be valid. However, it is uncertain that these obtained parameters match the material behavior for the specific machining process being simulated 13 . It is well known that these material parameters can differ from one article to another, sometimes under the same machining conditions, because there is no specific FEM model in each case or a general model. This might occur when, for example, the JC equation 15 and the fracture energy (G f ) are used to simulate the material fracture behavior. G f is associated with the material in order to characterize fracture during the cutting process (shear failure criterion 16 ). This value can be obtained from previous research, although it does not have a stable value, or through experimental methodology. However, this is time consuming and might not be economically and technically feasible.
So, this study aims to establish a relationship between the cutting parameters and the fracture energy to avoid being dependent on complicated experimental methods or values published in previous studies. This adapts the fracture energy value to the cutting parameters of the selected process. The cutting parameters implemented in machining processes highly influence the material behavior and so, the fracture energy should not be maintained as fixed but adapted. Especially in the machining processes, where extreme conditions are presented throughout the entirety of the cutting process and the material behavior is affected by it. Also, the simulation process can be better adjusted to the material behavior adapting the fracture energy, depending, in this case study, on the cutting parameters. One of the main objectives of the present work is to study the influence of the fracture energy on the machining process for different cutting parameters using FEM and to compare it with experimental results. The parameter G f has been selected since the JC fracture equation only defines the initiation of the crack. However, G f governs the damage evolution thus providing results that might fit better to the real behavior of the material. It should also be emphasized that the fracture energy is usually mode dependent and can therefore depend on the loading. The JC-model and JC damage criterion considers plasticity and damage initiation. This takes into account the strain, strain rate and temperature. The JC damage criterion has been commonly employed in simulation of segmented chip formation when machining titanium alloys like the one considered in this study. The main JC parameters that can be found implemented in the literature do not, in general, present a significant variation 11,17,18 but G f can deviate considerably depending on the fracture mode considered (from 11.53 MJ/mm 217 to 33.67 MJ/ mm 219 ). This material parameter is obtained from complex experimental tests 20 , taking the direction of the load application among others into account. However, the forces developed during the machining process are complex and so, its direct implementation in the numerical model can be complicated 21 . G f is considered in general as a constant of the material even though it is affected by the process conditions 22,23 .
It is well-known that an inherent property of this type of failure modeling produces mesh dependent results 24,25 . In the case of JC dynamic failure criterion, where elements are continuously degraded and deleted when the failure criterion is reached, the released energy depends on the element size. When a material point softens in continuum damage models, the softening becomes localized, promoting further softening. The mitigation of damage localization and mesh dependency is outside the scope of the current paper, being focused on the relation of the cutting parameters and the fracture energy. Therefore, the element mesh is kept fixed during simulations and the reported fracture energies are only valid for the specific mesh used in this paper.
The main intention of this study is to analyze the influence of the fracture energy mentioned above using the Finite Element Analysis (FEA) and to establish an empirical relationship between the G f , the cutting parameters and the force at fracture in the machining of the Ti6Al4V titanium alloy. The aim is to avoid being dependent on complicated experimental methods or values published in previous studies, presenting the fracture energy values that better adapt to the conditions simulated and fit the experiments conducted.

Materials and methods
As mentioned, the present work considers the study of the Ti6Al4V alloy machining process using finite element analysis and experimental tests. An empirical relationship based on the cutting parameters for the simulation model is established and fitted to the results. The empirical equation is also validated by comparing the simulations with the experiments performed for different ranges of cutting parameters. Experimental tests. The composition of the selected Ti6Al4V alloy is presented in Table 1 and is obtained by arc atomic emission spectroscopy (AES).  (Table 2) have been chosen based on the industrial requirements. A total of 180 tests were performed, with 15 repetitions for each cutting parameter combination to ensure repeatability. The experiments are the same as in previous studies made on this alloy by the authors [28][29][30] . This facilitates a quick evaluation of the experimental results, reducing experimentation times.
The tests were performed in a parallel lathe (ECLIPSE model, ALECOP), equipped with FAGOR 8055 T Numerical Control, in dry conditions. To be able to minimize the geometric variable influence, all tests were carried out in an orthogonal configuration. The specimens were designed with a tailored geometry to maintain orthogonal conditions throughout the tests. Different grooving operations were carried out on a billet (L = 170 mm, D = 105 mm) to achieve a tubular geometry. Then, two crowns were formed, corresponding to the two diameters machined previously. Each crown was machined with a specific thickness equal to cutting depth. Additionally, a relief zone was established, eliminating a sector of the crowns, to ensure that the spindle reached a permanent regime. Therefore, the aim was to obtain results close to a two-dimensional behavior with this geometry and to allow comparisons with the FEM studies carried out, as they are usually implemented to minimize computational time 18,29,31 . A new coated WC-Co insert was used for each test to maintain the same initial conditions. For the monitoring of the cutting forces (Fc), a piezoelectric sensor dynamometer, 9121 Kistler with an amplifier and dynamic signal analyzer and a Pulse Labshop data acquisition system were used.
The generated chips were collected, codified, photographed and stored after each test using metallographic techniques for analysis. For the study of the chips, an inverted metallurgical microscope (EPIHOT 280 NIKON, Tokyo, Japan) and a CF Optical System (1.5 × to 400 ×) for the SOM images were used. The corresponding measurements of the samples were obtained using a digital-image-processing software (Omnimet BUEHLER, Lake Bluff, IL, USA). A complete study and representation of the chip geometry is presented in 28 .
The experimental results indicate that a relationship can be obtained between Fc, f, v c , the height of the peaks (h p ) and the height of the valleys (h v ). For the validation of the model and the fracture energy equation, the results are compared with the results obtained from the FEM models.
Finite element analysis. For the Finite Element Analysis, a 2D FEM Lagrangian formulation model is implemented using Abaqus/Explicit. The general JC flow stress model is used for the description of the viscothermo-plastic behavior of the material [Eq. (1)]. This law considers in a multiplicative way the work hardening, strain rate hardening and thermal softening of the workpiece material at high cutting speeds. The JC law is usually applied in FEM studies of machining including new parameters 32,33 or not 11,12 . Equation (1) shows the JC law, where σ is the flow stress, ε is the equivalent plastic strain (representing the work hardening based on the constants A, B and n), ε the strain rate, ε 0 the reference strain rate (representing the strain rate based on the constant C), T r the ambient reference temperature and T m the melting temperature of Ti6Al4V (representing the thermal effect based on the constant m). The parameter A is the initial yield strength of the material at quasistatic strain rate. The parameter B and n represents the flow stress on strain hardening behavior at quasi-static strain rate. And the parameter C represents strain rate effect, and m represents thermal softening effect. Table 3 shows the JC model parameters implemented for this study 11 . The reference strain is defined as 0.7 s −117 .
Also, for the damage initiation, the JC fracture model or damage criterion 12,34 , given in Eq. (2), is used to define the initial failure strain. Table 4 shows the model parameters, where ε −pl oi is the equivalent plastic strain at the onset of damage (assumed to be dependent on ε −pl /ε 0 a non-dimensional plastic strain rate, p/q a dimensionless pressure-deviatoric stress ratio (being p the pressure stress and q the von Mises stress), T * the non-dimensional temperature defined previously in Eq. (1) (T − T r /T m − T r and d i are the failure parameters 11 . where D is the damage, σ the effective stress and u p the equivalent plastic displacement. When the damage progresses in a specific element and reaches D = 1, the element is considered fully damaged and is consequently removed from the analysis. The element deletion is activated, so the elements that reach the distortion marked for the simulation (G f ) are eliminated from the simulation.
For the contact between the chip and the cutting tool, a surface-to-surface contact definition is applied. A penalty friction formulation based on Zorev's model [Eq. (4)] is used to define the tangential behavior 35 where τ the frictional stress, µ the apparent friction coefficient at the cutting tool/chip interface, µN f the friction along the contact length and τ limit the material shear stress, remaining constant µ and τ limit (µ = 0.3 and τ limit = 331 MPa). The normal behavior is defined as a hard contact. However, the friction variable has been set based on previous analysis and it is not part of this analysis 36,37 . The heat transfer coefficient at the interface was modelled using a high thermal conductance in order to reach the steady state quickly. Figure 1 shows the mechanical boundary conditions for the model. As illustrated, the displacement is set to zero in the y-direction at the bottom of the workpiece. In addition, the cutting speed is applied in the x-direction at the bottom of the workpiece. The cutting tool is constrained to a Reference Point (RP) and is set to zero. With regard to the thermal boundary conditions, a starting temperature of 298 K is applied to the workpiece and the tool in the initial step.
According to the consulted literature, the mechanical and thermodynamic parameters for the workpiece and the tool are as presented in Table 5 32,38 .
For decreasing the computational time, the mesh of the workpiece is divided into two different domains (Fig. 1). The upper part, which is subjected to severe deformation due to the machining process, is meshed with smaller elements. The rest of the workpiece is meshed with larger elements further away from the cutting area. The tool is divided into two areas; the first one around the contact zone with the chip and the second far away from the cutting area. The element size and shape, as well as the meshing control for each area, are shown in Table 6. As for the boundary conditions, the bottom of the workpiece has been fixed and no movement is allowed in the x-direction or y-direction. As for the cutting tool, a reference point (RP, Fig. 1) has been placed on the upper part to be able to include the cutting speed in the x direction.
The general element type applied for the workpiece and cutting tool is CPE4RT (4-node plane strain thermally coupled quadrilateral, bilinear displacement and temperature, reduced integration, hourglass control). For a better transition between the areas with different elements size, the element type CPE3T (3-node plane strain thermally coupled triangle, linear displacement and temperature) is used. The model has a total of 6514 number of elements and 6639 number of nodes.
µN f , τ < τ limit (sliding) τ limit , τ ≥ τ limit (sticking) www.nature.com/scientificreports/ Also, in order to complete the simulations and consume as little computational time as possible, it is important to avoid severe element distortion as much as possible. Therefore, it is essential to create a suitable mesh. As mentioned above, it is also well-known that the results are mesh dependent for materials which exhibit softening 19,39,40 . The mesh dependency is unavoidable and therefore it is referred to as model data instead of true material data. However, using the same mesh configuration, size, elements, etc. with no remeshing, it is assumed that the results are valid for this model. An optimization study has been conducted to set the final mesh used in all the simulations.

Results and discussion
To obtain the empirical relation, a set of simulations are carried out using a set of G f values. In the literature, a fracture energy of 33.67 mJ/mm 2 is commonly used for this alloy 17 . However, it is not possible to maintain this value for the set of cutting parameters, because of the differences between the simulation and the experimental results, but also due to the simulation problems encountered with mesh distortion etc.
Having this parameter fixed for the wide range of cutting parameters simulated is also inadequate. This is mainly due to its dependency on the changes in the fracture mechanics of the cutting process. A G f adaptation for a low or a high speed, as well as for the different feeds, provides more suitable results, as can be seen in this section.
Several simulations are performed using different fracture energies to obtain the dependence of the cutting forces on the feed, cutting speed and G f . Due to the computational time requirements, the simulations are limited to two different feeds, as seen in Fig. 2. Also, in Table 7, the simulated chips are shown. Depending on the applied G f different behavior can be observed. All the parameters have been fixed except for G f and the feed.
Based on the cutting forces provided from the experimental results in 28 and presented in Table 8, the fracture energy is obtained from Fig. 2 using the corresponding cutting speed and feed. The empirical relation [Eq. (5)] is then obtained by applying a multivariable linear regression to determine the unknown parameters, presenting a R 2 of 0.91. By studying the results obtained from the experiments and the FEM model, the empirical relation for G f , the cutting parameters and the cutting force has been obtained [Eq. (5)]. This can be used to determine G f from experiments by measuring the cutting force from the experiments, using the corresponding feed and cutting speed. The final equation is given by Two more simulations are made with the predicted G f for a cutting speed of 65 m/min and a feed of 0.10 and 0.20 mm/r in order to confirm the results obtained from Eq. (5). The resulting cutting forces from each simulation give a good fit to the fracture energy given by Eq. (5), as it can be seen in Fig. 2 (v65f01 check and v65f02 check).   www.nature.com/scientificreports/ Once the equation is assumed to be valid, a set of simulations are performed using the cutting parameters in Table 2. From the obtained results, the cutting forces and the chip geometry are compared with the ones acquired from the experimental tests. Figure 3 shows the comparison of the cutting forces for the implemented range of cutting parameters. It can be seen that the empirical relation works good for the implemented feed range. Table 9 shows the calculated G f and the simulated and experimental cutting forces and the corresponding error for each couple of cutting parameters. It can be seen that the error is in general below 10%. Only for a 125 m/min and 0.05 mm/r, the simulation result error is over 10%, which can be explained due to the mechanical behavior of the Ti6Al4V at a higher cutting speed. This type of alloy presents a low thermal conductivity which leads to a rapid increase in temperature in the primary deformation zone, having an adiabatic process. So, the material softens and it is easily deformed 29,41 . Also, as v c increases, friction is reduced and so the cutting forces are lower 42 . These kinds of mechanisms are   www.nature.com/scientificreports/ not included in such an empirical relation as the one presented in Eq. (5). However, the main objective of this analysis is to work with a simple equation that can be used to determine the fracture energy from cutting forces and cutting data. Knowing its limitations, it is presented as a first approximation. In Fig. 4, G f is presented as a function of feed and cutting speed. It can be clearly observed that the fracture energy tends to converge to a specific value at a higher feed. This tendency also applies for higher speeds. This demonstrates that less fracture energy is needed to fracture the material at higher values of feed and cutting speed. So, Eq. (5) shows how it is necessary to adapt Gf depending on the cutting data choice. This implies the importance of having a mixed-mode Gf model due to the complex interaction of loading, geometry etc. as the fracture mode varies in the machining process depending on these parameters 43 .
In order to illustrate the obtained results more clearly, a geometry comparison of the chip formation process and the experimental results 28 has been carried out. The results of the study are presented in Table 10. It can be www.nature.com/scientificreports/ seen that the cutting forces are in good agreement with the experimental results. The model presents similar chip formation for higher range of feeds (0.30-0.20 mm/r). However, it is not possible to capture the serrated geometry for the lower range of feed (0.10-0.05 mm/r). This might be explained due to the changes in the fracture mechanism. For higher f, the cutting process develops higher shear forces that can be well represented with the JC formulation implemented. However, for lower f, the portion of material that is submitted under the cutting tool thrust is less (thinner). Because of this, the normal force is bigger than the shear force and the model is not able to present the characteristic serrated geometry. It can also be observed that the mesh consists of fewer elements in the chip zone resulting in a stiffer behavior. This can be improved with a thinner mesh but, as explained in previous paragraphs, a thinner mesh, brings more elements and that more time is needed for the simulation. So, for this first approximation, the mesh is maintained. JC material model in simulations can lead to continuous chip geometry for low v c and f 17,44 but is one of the criteria generally used to simulate the crack initiation and propagation into the primary shear zone during the chip formation. The literature shows several numerical models developed but none of which are adequately adapted to predict the chip formation profile which highlights the difficulty behind the replication of this mechanism and the existing gap 4 . Also, the chip formation mechanism of Ti6al4V is still not well understood due to the complexity of the material behavior itself 42 . It can be seen in the literature, that some authors consider that chip formation in this kind of alloy is due to a thermoplastic instability, while others consider the initiation and propagation of cracks inside the primary shear zone of the workpiece material. Moreover, the microstructural state of the alloy strongly influences the chip formation. This alloy presents a transitional chip at lower v c , showing a more segmented chip at higher v c 41 , as it shows in the experimental tests. As stated in previous research 28 , over the experimental tests performed the serrated chip morphology remains continuous across the v c and f range studied, due to the high plasticity levels of this alloy and its low thermal conductivity, resulting in thermal softening, making the chip more difficult to break. This kind of behavior and the characteristics explained in previous paragraphs are what make it difficult to translate the alloy behavior to the simulation experiments.  www.nature.com/scientificreports/ Taking into account this brief explanation, it can be understood that the simplicity of the model is not taking into consideration all the parameters that influence the chip formation for this alloy. However, this first approximation of the model fits well to high feeds but not to low feeds. Therefore, the numerical model can be improved in this aspect. It can be identified in the simulations that, for the lower f of 0.05 mm/r, the chip tends to have a serrated shape with an increase in v c . For 125 m/min and 0.05 mm/r, the chip has a more pronounced wavy shape than for lower v c , but does not create a serrated shape. The same happens, less obvious, for f = 0.10 mm/r. In this case, the wavy effect disappears. This effect depends, in this case, on the mesh implemented. For f = 0.10 mm/r it seems that there are not enough elements to appropriately represent the chip shape. For f = 0.05 mm/r, the distortion of the mesh is such that some spaces and jumps between elements appear. This behavior, although is not well represented, follows the behavior presented by the alloy during the experimental tests. In the lowest feed rate range (0.05-0.10 mm/r), the chip was continuously helical and showed a tendency to create chip nests for 0.05 mm/r. This could be understood with the less serrated morphology of the chip obtained during simulation.
The geometrical data obtained from the experimental tests are shown in Table 11 28 and the comparison between the average h p and h v is presented in Fig. 5 and Table 12. This is only shown for 0.20 and 0.30 mm/r because, as shown in Table 10, for the other f, the chip geometry does not present a serrated shape. The average results of h p and h v present a larger error compared with the ones obtained from the experimental tests (Table 12). Considering the spread in the experimental measurements (maximum and minimum experimental values, as shown in Fig. 5 as "min" and "max") and not only the average value ("exp"), but it can also be observed that the results obtained from the simulation are in general within the range of the experimental results. The maximum error is 43.5% for v c = 30 m/min and f = 0.20 mm/r and the minimum is 8.5% for a v c = 60 m/min and a f = 0.30 mm/r. It can be estimated that with an increment of the v c and the f, the error % is reduced, except for www.nature.com/scientificreports/ 125 m/min, where the trend changes. This can be attributed to the same mechanism explained in relation to Table 10. As an average the model presents an error of 12% and for the valley values and 30% for the peak values. It can be concluded that the initial model provides reasonably good fit to the experimental results regarding the obtained cutting forces. For the chip geometry analysis, the model can be valid for a feed range between 0.20 and 0.30 mm/r. Thus, it is necessary to develop the model further to better adapt to the low feed range (0.05 and 0.10 mm/r).
As a model control, the energy balance, which must be identified, quantified and minimized to prevent erroneous findings, is analyzed, presenting a 5% of error, intended to be improved in future work. To control that the numerical model adhere to the basic physical laws (conservation of energy), the global energy of the model is checked to analyze that there are no major inconsistencies in the energy of the system, e.g. keeping the sum of internal, kinetic, sliding, hourglass, system damping, and rigid wall energies within an acceptable range. According to the literature consulted, this acceptable range should be around a 5% of the total global energy 45,46 .
Nevertheless, the developed model is presented as a preliminary study for the range of values studied. It would require a more in-depth study to analyse and incorporate the fracture mechanism mentioned, by which the chip generation can be closer to reality. A mixed-mode fracture criterion with different fracture energies in tension/ peel and shear should be implemented and tested to further evaluate the model.

Conclusions
In this work, a 2D FEM model for the analysis of the machining process of the Ti6Al4V alloy with adapted fracture energy values for each cutting parameter set has been developed. Additionally, an empirical equation to obtain the optimum fracture energy is presented. This involves maintaining the simulation parameters fixed and adapting G f depending on the cutting parameters. A validation of the FEM model with experimental results obtained in previous studies has been made in order to establish the methodology and select the optimal values. This shows the advantages and disadvantages of the analysis. For this comparison, the Fc, h v and h p are selected to analyze the model mechanical and geometrically.
Implementing the G f obtained from the empirical equation developed, the simulation model shows results in good agreement with the experimental tests. The main comparison is made with the cutting forces, presenting a value error lower than 10% for most of the cutting parameters studied. For a v c of 125 m/min the data evolution presents a higher error, but only for a f of 0.05 mm/r this percentage is over 10%. This is due to the mechanical behavior of the Ti6Al4V at higher v c (low thermal conductivity, high chemical reactivity and high hardening) and the simplicity of the equation.
As for the chip geometry, the model is able to represent the serrated characteristic chip for high f (0.30 and 0.20 mm/r). However, this shows an average error between peaks and valleys of 12% and 30% correspondingly. For lower f, the model is not able to obtain a serrated chip, due, again, to the simplicity of the initial model. The fracture mechanism changes depending on the f implemented and the model is not developed to implement these variables in the first stages of the study.
This work attempts to highlight the necessity of adapting the fracture energy to the machining process parameters during FEM simulations, due to the influence that this parameter has upon the mechanical behavior of the material being machined. The aim is not to obtain the value from previous literature works or complicated experimental methods, which hardly represent the mechanism behind the machining process analyzed.
It is necessary to point out that this work is only the first stage of the cutting forces and chip morphology analysis of the Ti6Al4V alloy by FEM. The study and results are established for the specified cutting conditions but with an appreciable adjustment to the experimental values presented. Other mechanical and geometrical parameters, such as the temperature, shear stresses, chip segmentation ratio or shear angle, will be addressed in further works. In addition, the generality of the empirical equation should be tested in a wider range of cutting speed and cutting depth. www.nature.com/scientificreports/