Coupled Eulerian–Lagrangian simulation and experimental investigation of indexable drilling

In many industries, indexable insert drills are used to cost effectively produce short holes. However, common problems such as chatter vibrations, premature tool wear and generation of long curly chips that cause poor chip evacuation make optimization of the drilling process challenging and time-consuming. Therefore, robust predictive models of indexable drilling processes are desirable to support the development towards improved tool designs, enhanced cutting processes and increased productivity. This paper presents 3D finite element simulations of indexable drilling of AISI4140 workpieces. The Coupled Eulerian–Lagrangian framework is employed, and the focus is to predict the drilling torque, thrust force, temperature distributions and chip geometries. To reduce the computational effort, the cutting process is modelled separately for the peripheral and the central inserts. The total thrust force and torque are predicted by superposing the predicted result for each insert. Experiments and simulations are conducted at a constant rotational velocity of 2400 rpm and feed rates of 0.13, 0.16 and 0.18 mm/rev. While the predicted torques are in excellent agreement, the thrust forces showed discrepancies of 12 - 20% to the experimental measured data. Effects of the friction modelling on the predicted torque and thrust force are outlined, and possible reasons for the thrust force discrepancies are discussed in the paper. Additionally, the simulations indicate that the tool, chip and the local workpiece temperature distributions are virtually unaffected by the feed rate.


Introduction
Conventional drilling is the major hole-making process used in industry today, and the drilling process can be divided into different types such as twist, step, spade and indexable drilling. Solid twist drill is the most common one, but indexable drilling is a cost-effective alternative to produce short holes. Indexable drills consist of multiple inserts, usually two, mounted on a tool holder. The inserts can be indexed and replaced when they are worn out and no regrinding is necessary. Since the inserts can be almost tailor-made for its purpose, they are also suitable for difficult-to-cut materials.
Indexable drills have, under certain conditions, shown to produce less oversize holes and with better surface finish compared to twist drills [1]. Also smaller thrust force and torque compared to twist drilling were observed which yielded improved tool life and reduced power consumption. In [2], the authors conducted an experimental study on drilling under MQL conditions on the steel AISI 1045 and three difficult-to-cut materials: stainless steel AISI 304, titanium Ti64 and Inconel I718. The experiments were repeated using both an indexable drill and a solid twist drill, and the measured thrust force, inner surface roughness, tool wear and tool temperature were compared. They found that for Ti64 and I718, the thrust force was smaller and the tool temperature was lower for the indexable drills. However, for AISI 1045 and AISI 304 the tool temperature was higher for the indexable drilling. The thrust force for indexable drilling was smaller for AISI 1045 but larger for AISI 304. For the difficult-to-cut materials, similar and relatively low surface roughness values were obtained with the two drills. However, for AISI 1045, indexable drilling resulted in significantly worse surface quality and it is argued that this is explained by the built-up edge (BUE). The BUE phenomenon during indexable drilling has also been studied experimentally in [3].
Due to the asymmetric tool geometry, and often also different insert geometries, unbalanced process forces will be present during the drilling process. In [4], the authors investigated three variations of indexable drilling by introducing splitting grooves, change in point angle, cutting edge angle and relief angle. They found that both the surface roughness and hole roundness were improved and that the vibrations could be reduced. The analysis performed in [5] shows that the radial forces lead to both static and dynamic deflections that is affecting the shape of the drilled hole. The forces can also cause regenerative chatter that leads to unpleasant noise levels, damage on workpiece and tool, poor surface quality and reduced tool life [6]. Furthermore the workpiece temperature, clamping loads and drill run-out can also significantly impact the hole quality [7,8].
Decision support regarding the tool design and optimal process parameters are usually based on time-consuming and costly experimental testing campaigns, and in situ measurements during the drilling process are obviously precluded by the inherent visual limitations. The wide range of process states such as forces, vibrations, surface quality, workpiece and tool temperatures, tool stresses and wear, clamping loads, residual stresses and their complex interactions that governs the quality of the cutting process are very difficult to predict. Therefore, modelling and simulations of metal cutting processes have, during the three last decades, been subjected to intensive research, and nowadays, the finite element method (FEM) is a common approach used for simulating metal cutting processes. However, a machining process is very challenging to simulate due to the extreme local conditions with high strains, strain rates and temperatures. A relatively recent overview of the current state and needed developments within the field is given in e.g. [9]. In most cases, the Lagrangian description in which nodes follow material points is utilized, but due to the extreme localized deformations and the associated numerical difficulties, alternative FE formulations are commonly employed such as updated Lagrangian, Arbitrary-Lagrangian-Eulerian (ALE) [10] and the Coupled-Eulerian-Lagrangian (CEL) [11] as well as meshless methods like the Smoothed-Particle Hydrodynamics (SPH) [12].
Numerous numerical investigations have been reported on twist drilling using FEM with various formulations. For example, thrust forces, torques and chip morphologies were predicted reasonably well for small-hole drilling using the Lagrangian formulation with explicit time integration in [13]. In [14], the tool temperature was predicted using the updated Lagrangian approach with the purpose of identifying machining parameters that enhances the tool life. Recently the drilling torque, thrust force and temperature generation were successfully predicted in [15]. In [16], the thrust force and the torque during drilling of carbon fibrereinforced plastic was predicted for various fibre orientations using the Lagrangian approach with explicit time integration and cohesive elements that model the fracture processes. For drilling in cortical bone, the torque and thrust force where predicted in [17] using the Lagrangian formulation with explicit time integration and in [18] the updated Lagrangian formulation was used to predict the bone temperature. Since FE-simulations of drilling process are computationally demanding, other FE-based techniques have been proposed to predict the final shape of drilled holes where the actual cutting process is substituted by equivalent forces and heat flows to the workpiece, e.g. [19,20].
In this paper, the CEL formulation is used to model an indexable drilling process. In the CEL formulation, the tool is modelled using the Lagrangian formulation, whereas the workpiece is modelled using the Eulerian formulation. Since the workpiece mesh is fixed in space and the workpiece material flows through the mesh, mesh distortion is completely avoided and the chip formation process is modelled through the Eulerian/Lagrangian contact without the use of a failure criterion. In recent years, the CEL formulation has been increasingly used to model various metal cutting processes such as orthogonal cutting of various workpiece materials, e.g. [11,[21][22][23], slot and shoulder milling of Al6061-T6 [24], intermittent longitudinal turning [25] and hard milling of AISI4140 [26]. The authors in [27] employed the CEL formulation to predict the drilling torque, burr formation and residual stresses in aerospace grade aluminium and titanium for various cutting parameters and the obtained results were benchmarked with an equivalent Lagrangian model. The CEL approach showed higher computational efficiency and better ability to predict chip and burr formation as well as the torque.
Although a lot of research has been conducted on numerical modelling of drilling processes, research studies on indexable drilling is scarce and, to the best of the authors knowledge, numerical modelling of an indexable drilling process has never been reported in the literature before. However, indexable drilling has well-known problems such as chatter vibrations, premature tool wear and poor chip evacuation that precludes optimization of the drilling process. Therefore, predictive modelling of indexable drilling processes is desirable for supporting the development towards improved tool designs, enhanced cutting processes and increased productivity.
In this paper, we use the CEL formulation in Abaqus/ Explicit 2018 to model indexable drilling of AISI4140 and the focus is to predict the torque, thrust force, chip geometry and temperature distributions. To reduce the model size, the cutting process of the peripheral and central insert are modelled separately and the thrust force and torque are predicted by superposing the predicted result for each insert. This reduces the computational effort associated with each simulation. Through small scaled parametric studies, the influences of the friction, mesh size and mass scaling parameters on the simulation results are also investigated.

Experimental work
Dry drilling experiments are carried out in a GROB G300 4-axis machining centre to measure the drilling forces and torque. The experimental series is conducted with a variable feed rate but with constant rotational speed. The rotational speed is 2400 rpm and four repeated tests are conducted for each feed rate of 0.13, 0.16 and 0.18 mm/rev. The selected rotational speed of 2400 rpm results in a cutting velocity of about 3 m/s at the peripheral corner edge which is the maximum cutting speed recommended by the tool manufacturer. After each experiment, the insert is indexed and after four tests the insert is replaced.
Cylindrical specimens of height 25 mm and diameter 32 mm are cut out from quenched and tempered AISI 4140 (42CrMo4) rods. The drilling tool KOMET XV62 KUB-DUATRON, shown in Fig. 1a, gives the orientation of the two identical inserts as visualized in Fig. 1b. The outer diameter of the tool holder is 22.7 mm and with inserts placed in the holder the drilled diameter is 23.7 mm. The rake angle is 0 • for both inserts if the chip breaker groove is ignored and the clearance angles are both 8 • . The central and peripheral inserts are slightly inclined from the rotational axis such that the point angle is roughly 172 • where the back taper angles are 5 • and 3 • for the central and the peripheral insert, respectively. The corner angles of the insert are 90 • which means that the back taper angles are equal to the back taper angle of each insert. Square shaped KOMET inserts, SOEX 07T308-1 are used that have the side length 7.94 mm and the nose radius 0.8 mm. The inserts have the material code BK8425, and they are coated with a multilayer PVD-TiAIN/ TiN coating.
The experimental setup is shown in Fig. 2a, and its components are for clarity schematically illustrated in Fig. 2b. The specimen is clamped between two steel plates by tightening four bolts. That is, the specimens are held in place by contact and frictional forces that act on the top and bottom surface of the specimens. Interference between the tool and the upper plate is avoided by preparing the upper plate with a 25 mm hole. A HBM MCS10 multicomponent transducer (MCS) is connected to the lower plate to measure three perpendicular force components and the drilling torque. Four HBM ClipX measuring amplifiers, a National Instruments USB 6210 data acquisition card and the data acquisition tool box in MATLAB were used to amplify, acquire and analyse the experimental data. The bottom of the transducer is connected to an additional plate that is fixed to a vice which in turn is mounted to the machining table. Besides the forces and torque, the temperature and in-process displacements are measured using four thermocouples and a laser scanner, but the results from these measurements are not reported in this paper. After each experiment, the chips produced by the central and peripheral inserts are collected and photographed to enable comparisons with the simulated counterparts.

Finite element models
Abaqus/Explicit 2018 and the CEL formulation is used to model the drilling process. The cutting process of the central and peripheral insert are modelled separately. Figure 3 shows the initial configurations of the FE-models each consisting of two parts; i) a rigid thermally coupled Lagrangian model of the cutting insert; ii) an Eulerian space that is partly filled with workpiece material. A reference node placed at the hole centre is controlling the motion of the cutting insert and it is subjected to a prescribed rotational velocity about the hole centre axis while the remaining degrees of freedom are suppressed. The Eulerian spaces, indicated with dashed lines, are constructed to extend well beyond the workpiece material boundaries to allow the chip to move and deform. The outer free surfaces of the Eulerian spaces are subjected to fixed boundary conditions and they are for simplicity modelled as adiabatically insulated and an initial temperature of 20 • C is defined for both the insert and the Eulerian space.
Only a 55 • and a 30 • rotation of the drilling tool is simulated for the central and peripheral insert, respectively, in order to achieve reasonable computational times while still obtaining steady-state cutting forces and workpiece temperatures near the cutting zones.

Euler space and workpiece material boundaries
The workpiece material boundaries shown in Fig. 3, are carefully constructed to model the nominal chip thickness along the cutting edges. Cross sections of the workpiece material boundaries are depicted in Fig 4. As shown in Fig. 4a, two regions along the central insert cutting edge can be identified. The region near the centre of the hole is exclusively cut by the central insert giving a constant chip thickness equal to the feed per revolution. Outside this region, the hole is cut by both inserts resulting in a gradually decreasing chip thickness. This mixed central/peripheral region is indicated with the dashed lines that shows the position of the peripheral insert a half rotation earlier. The corresponding analysis for the peripheral insert is depicted in Fig. 4b.  Figure 5 illustrates the nominal chip thickness distribution along the cutting edges for the two inserts.
The initial material boundaries are finally constructed by sweeping the grey regions in Fig. 4 about the rotational axis. The Lagrangian insert is at the start of the simulations placed such that the dark grey regions represent the uncut chip cross sections.
The Eulerian space for the central insert cutting simulations shown in Fig. 3a has a 10 mm radius and a 3 mm height. It is discretized with an unstructured mesh consisting of about 1.6×10 6 8-noded thermally coupled Eulerian linear brick elements (EC3D8RT) with the dimensions 60 -80 m in-plane and 30 m in the feed direction.
For the peripheral insert cutting simulations, the Eulerian space shown in Fig. 3b is a circular trapezoid of height 3 mm and with an inner and outer radius of 4 mm and 15 mm, respectively. It is discretized with a spider web mesh consisting of about 1.7×10 6 elements. The dimensions of the elements are about 80 m in the radial direction, 30 m in the feed direction and along the tangential direction, the element size increases along the radial direction from 35 m to about 120 m near the outer edge.

FE-model of the cutting insert
As shown in [28], the micro-geometry of the insert has a major influence on the simulated results in terms of, for example, cutting forces, chip form and temperature distributions in the insert. Therefore, special attention is given to model the precise geometry of the insert. First, the insert was scanned using the GOM ATOS III Triple Scan optical 3D scanner system to extract a STL-mesh representation of the insert geometry as shown in Fig. 6a. Subsequently, the solid FE-mesh shown in Fig. 6b was generated by using the commercial pre-processor

Material and modelling parameters
Plastic flow of the workpiece material is modelled using the Johnson-Cook (JC) constitutive model [29] where T, ̄p , ̇̄p and ̇0 represent the temperature, equivalent plastic strain, equivalent plastic strain rate and reference strain rate, respectively. The governing material parameters for AISI4140, shown in Table 1, are taken from the literature. As enabled by the Eulerian formulation, chip formation is simulated without the use of initiation and evolution of material damage.
A rigid body constraint is imposed on the cutting tool nodes to reduce the computational times while still allowing for heat transfer across the tool/workpiece contact surfaces. The thermo-mechanical parameters for the cutting tool and workpiece in Table 2 are adopted from [31] where the same alloy and a similar insert material was used.
The default Eulerian-Lagrangian general contact algorithm is employed to track the contact interfaces between the tool and the workpiece material. The frictional shear stress, f , is governed by a sticking-sliding friction model which is based on Zorev's friction model in [32] and it is given by Eq. (2), where , p and p denotes the friction coefficient, the contact pressure and the frictional shear stress limit. As in [33], the sliding friction is modelled with Coulomb friction with = 0.4 that is independent of both temperature and relative sliding velocity. In the sticking region, the frictional shear stress equals the shear stress limit, p , which is chosen using the approximation outlined in [10] for orthogonal cutting conditions. Assuming a state in the sticking region where T = 1000 • C, ̄p = 4 , ̇̄p = 10 5 gives a flow stress of ̄= 583 MPa using the JC-parameters in Table 1. With the known von Mises stress at these conditions, the corresponding shear flow stress is set to p =̄∕ √ 3 ≈ 337 MPa. It should, however, be noted that this choice of shear stress limit is somewhat arbitrary since both the temperature and strain rate differs between the inserts and for each insert, they also vary greatly over the tool/workpiece contact surfaces.
Both plastic work in the workpiece material and frictional work are assumed to be completely dissipated into heat, i.e. both heat fractions are set to unity. Moreover, the frictional heat is partitioned equally into the tool and workpiece contact surfaces. Heat exchange across tool-workpiece contact surfaces is governed by the pressure-dependent heat conduction given in [34], and the data are for clarity repeated here in Table 3.

Mass scaling
The simulations presented here are computationally expensive due to the large number of elements and small element sizes that results in very small critical time increments. In [35], the authors investigate two strategies for reducing the simulation time in thermo-mechanical coupled CEL analyses, i.e. time scaling and mass scaling. Since the presented FE-simulations involve a wide range of strain rates and a strain rate dependent flow stress, it is only feasible take advantage of the mass scaling technique. As the cutting insert is rigid, and therefore excluded in the critical time step calculation, only the density of the Eulerian elements is scaled.
As outlined in [35], the density is replaced by a fictitious density ⋆ = k m which increase the critical time increment with a factor of √ k m . This means that the cpu time is reduced by a factor of √ k m . The thermal parameters are then compensated for by transforming the specific heat capacity to c ⋆ p = k −1 m c p . This strategy is applied on the central and peripheral cutting simulations for which the rotational velocity and feed rate are set to 2400 rpm and 0.18 mm/rev, respectively. The root mean square (RMS) values of the steady-state torques, M c , thrust forces, F D   Tables 4 and 5. An increasing scaling factor results in monotonically increasing torques and thrust forces for both the central and peripheral cutting simulations. Moreover, it can be noted that the torque is more affected by mass scaling than the thrust force and, as expected, the peripheral cutting simulation is more sensitive due to the higher velocities involved.
Although the level of the steady-state torques and thrust forces, here given as RMS values, are monotonically increasing with the scaling factor, the evolution of all curves are very similar. For the tested scaling factors, no significant difference in amplitudes of the periodic oscillations could be noted. Moreover, mass scaling with the tested scaling factors has virtually no effect on the predicted temperature fields and chip morphologies.
To improve the computational efficiency while avoiding too large inertial effects from increasing the mass of the Eulerian elements, k m =50 and k m =10 are used for the simulations of central insert cutting and peripheral insert cutting, respectively.

Mesh sensitivity
In [36], guidelines are provided for selecting mesh parameters in CEL analyses of orthogonal cutting. About eight elements across the uncut chip thickness is recommended to obtain a converged chip morphology and cutting forces. It is, however, noted that chip morphologies and RMS values of the cutting forces are very similar also when using four elements across the nominal chip thickness even though larger oscillations on the force curves are obtained. Moreover, it is recommended to use an element size smaller than half the cutting edge radius. Although these guidelines are most likely not completely transferable to the present work, they serve as a basis for selecting the element discretizations described in Sect. 3.1. To avoid excessive computational times, an element height of 30 m was selected for both the central cutting and peripheral cutting simulations. That is, the number of elements across the uncut chip thickness ranges from 4 to 6 for the simulated feed rates of 0.13 mm/ rev and 0.18 mm/rev, respectively. It should, however, be mentioned that the number of elements across the uncut chip thickness is smaller within the mixed central/peripheral zone. Additionally, the element height of 30 m corresponds to about 75% of the cutting edge radius.
A simple mesh sensitivity analysis was conducted to ensure that the discretization of the Euler space presented in Sect. 3.1 was able to produce converged results. These simulations are computational demanding which precludes a more complete mesh convergence study where all element dimensions are allowed to change size in a systematic manner. Therefore, this study was limited to investigate the effect of varying the element height while keeping the remaining element dimensions constant. The mesh presented in Sect. 3.1 is here denoted the base mesh. For the fine and coarse mesh, the element height is set to 20 m and 40 m which corresponds to about 50% and 100% of the cutting edge radius.
In Tables 6, 7, 8, and 9, the obtained RMS values of the torques and thrust forces during steady-state cutting are presented as well as the corresponding percentage difference compared to the results obtained when using the base mesh.
As shown in Tables 6, 7, 8, and 9, the selected element heights give very similar RMS torques and thrust forces. Also, the cyclic variation of the curves are very similar when using the different element meshes. That is, the amplitude of the oscillations of the curves about its mean value during steady state is virtually the same for the different meshes.
For the central cutting simulation with feed rate 0.13 mm/ rev, see Table 6, a decreasing RMS torque is obtained as the   element height is decreasing although this effect is barely visible. For the thrust force curves, an unexpected effect is identified where the largest RMS value is obtained when using the base mesh. The same behaviours are recorded for the central cutting simulation when the feed rate is set to 0.18 mm/rev, see Table 7. The results from the peripheral cutting simulations, see Tables 8 and 9, show virtually no effect on the torque or thrust force curves from changing the element height. From a visual comparison of the field plots produced using the different element meshes, it can be concluded that changing the element height within the specified range does not influence the predicted temperature fields or chip morphologies. This simple mesh sensitivity study indicates that converged results are produced when using the base mesh presented in Sect. 3.1. All simulations presented in the following sections are therefore conducted using the base mesh.

Comparison of experimental and simulated thrust forces and torques
The torque and thrust force are computed in an uncoupled fashion by superposing the results from the central and peripheral simulations. In Fig. 7, the predicted drilling torques and thrust forces are compared with the experimentally measured counterparts. For each feed rate, the solid blue bar indicates the average of the measured RMS values and the orange bar the simulated RMS value. Moreover, the nonefilled bars show the contribution from the respective inserts.
As indicated with the black error bars, i.e. the measured range of RMS values, the repeatability of the experiments is satisfying. As shown in Fig. 7a, an excellent agreement between the measured and predicted torque is obtained for all feed rates with a maximum error of about 3%. The simulations predict that the peripheral insert is responsible for about 70% of the total torque for all feed rates.
The predicted thrust force agrees relatively well with the experiments, see in Fig. 7b, and the simulations show that the central insert is responsible for a somewhat larger part, ∼55%, of the thrust force. However, in contrast to the experiments, the simulated thrust forces are virtually insensitive to the feed rate within the tested range of feed rates. Therefore, the underestimation increases from about 12% for the lowest feed rate to about 20% for the highest feed rate.
Even though the predicted thrust forces are in acceptable agreement with the experimental data, the discrepancy deserves an investigation. As pointed out by the authors in [37], underestimated feed forces is a systematic error in finite element simulations of metal cutting and it is usually considered to stem from suboptimal friction modelling. In their investigation, conducted on 2D turning simulations, the friction model had a significant influence of 20-50% on the predicted feed force when comparing results obtained with frictionless contact, Coulomb friction(without a shear stress limitation), shear friction and various hybrids of thereof. In [10], Arrazola and Özel investigated how the cutting force and feed force were influenced by the choice of friction model in 2D plane strain ALE simulations of orthogonal cutting. With the friction coefficient =0.5 held constant, simulations were performed with both the Coulomb friction  model and the sticking-sliding (with a shear stress limit) friction model with various shear stress limits. It was concluded that both the cutting force and the feed force were increasing with an increasing shear stress limit and that this effect was most pronounced on the feed force. However, when the shear stress limit was increased beyond about half the static yield stress, both the predicted temperature fields and forces were similar to those obtained with the Coulomb friction model. For the present simulations, the influence of the friction model on the predicted torque and thrust force is evaluated for simulations with feed rate 0.18 mm/rev. It is clear that changing the friction parameters affect the central and peripheral simulations differently. As described later in Sect. 4.3, the inserts are subjected to different contact conditions and hence also different friction mechanisms are dominant. In Table 10, the results obtained with the sticking-sliding friction model presented in , see Fig. , are compared to those obtained with the Coulomb friction model with various friction coefficients. Since the shear stress limit was set to ∼57% of the static yield stress, the sticking-sliding model and Coulomb model both with =0.4 would yield similar results if they were adopted in the 2D ALE simulations described in [10]. However, for the present simulations, the effects of removing the shear stress limit are substantial and the predicted torque for both the central and peripheral insert increases by ∼18%. For the thrust force, however, the contribution from the peripheral insert increases by ∼12% while an ∼ 8% reduction is noted for the central insert. Therefore, the total effect of removing the shear stress limit is that the torque increases significantly by ∼18%, while the thrust force only increases by just over 1%. Furthermore, as shown in Table 10, simulations performed with the Coulomb model and other values of prove unsuccessful to generate improved prediction capability of the torque and thrust force.
Also for the sticking-sliding model, the influence of on the torque and the thrust force was investigated. Simulations were conducted with the feed rate set to 0.18 mm/rev and in Table 11, the results obtained with =0.4, see Fig. , are compared to those obtained using other friction coefficients. As expected, both the torque and thrust force are monotonically increasing with increasing friction coefficients. However, the influence of the coefficient of friction is very weak because of the shear stress limit. Furthermore, it can be noted that increasing the friction coefficient has a larger effect on the thrust force than on the torque.  These results indicate that it is not possible to remedy the missing thrust force, while maintaining the predicted torque, only by modifying parameters governing the friction models investigated here. It should, however, be noted that these types of constant friction models may not sufficiently represent the complex contact conditions between the micro-geometry tool edge and the workpiece, see [34]. In their investigation, the pressure-dependent shear friction model was implemented in a 3D model of longitudinal hard turning and proved to be more successful in predicting the process forces.
Indeed, it is difficult to identify a specific modelling choice or model parameter as the reason for the thrust force being underestimated. Rather, several factors influence the thrust force and many of them can be addressed to enhance the force prediction. For example, the authors in [37] investigate the effects of different modelling parameters on the feed force in 2D simulations of turning processes and the tool edge geometry followed by tool ploughing effects were found to have the strongest effects on the feed force prediction. Moreover, modelling of the initial flank wear that takes place within the first seconds of an experiment also contributes by about 20-50% increase of the feed force. Therefore, the authors argue that the inserts should be modelled based on the scanned tool geometry after the experiments rather than before the experiments have been conducted. It should therefore be noted that the present work does not take these key aspects into account. The insert was modelled based on scanning process presented in Sect. 3.2 and the scanning system had a resolution of 0.01 mm. This may be too low to reproduce the tool edge radius with a high precision. In addition, a higher mesh density of the tool cutting edge and workpiece may be required to model the complex contact conditions. Moreover, the initial wear was disregarded since scanning was performed on fresh inserts. The relevance of these aspects on the present simulations is indicated by the fact that the predicted thrust forces are virtually insensitive to the feed rate. That is, the ploughing effect seems to govern the predicted thrust forces and this is also reasonable when considering the relatively small uncut chip thickness as compared to the tool edge radius.
Another obvious aspect that is crucial for accurate prediction of the torque and the thrust force is the set of JC constants. In [38], the authors investigate the effect of the selection of JC constants on the cutting force and feed force in orthogonal turning simulations by testing five different sets of JC-parameters from the literature commonly used to simulate cutting of AISI 316L. It was shown that changing the set of JC constant can affect the feed force dramatically, with an increase of over 50%, while only modest effects were recorded on the cutting force with an increase of ∼7%. In this work, the JC constants were taken from the literature and that may constitute a partial explanation for the missing thrust force. In [25], the same JC constants were used in CEL analyses of longitudinal turning. Also in that work, the simulated feed force were substantially smaller than the experimental ones while the simulated and experimental cutting forces were in good comparison.

Comparison of experimental and simulated chip geometries
Examples of the chips produced during the experiments are displayed in Fig 8. The central chips are spiral shaped and the peripheral insert produce chips with a wavy pattern with just a few waves that promotes easy chip evacuation. These two types of chip morphologies are typical for indexable drilling with two inserts. Chatter vibrations were present during all experiments which is also reflected in the chip morphologies. All chips produced are somewhat serrated and this is more pronounced in the thicker part of the chips. The serration of chips cannot be modelled using the present modelling approach of a CEL analysis without damage evolution.
To make a comparison of the simulated and experimental geometry of the central chips, the outer diameters of the central chips are measured. These measurements give radii of curvature of about 4.0 mm to about 4.5 mm for the highest and lowest feed rate, respectively. As for the experimental chips, the radius of curvature of the simulated chips decreases with an increasing feed rate. In Fig 9, the simulated central chip with a feed rate of 0.18 mm/rev is shown. The outer radius of curvature is measured by transferring the 3D geometry to a CAD program, as shown in Fig. 10. A plane is created near the top of the chip and the three points are chosen on the outer rim of the chip to produce a circle. For these measurements, it should be emphasized that the result depends to some extent on the selection of the three points. The radius of curvature of the simulated chip is measured to approximately 1.4 mm which corresponds to only ∼35% of the experimental value.
The simulated peripheral chip for feed rate 0.18 mm/ rev is shown in Fig. 11. The simulated peripheral chips are found to be slightly cone-shaped, but it is not possible to make an accurate comparison with the experimental chips. This is due to the CEL formalism and the sticking selfcontact that precludes modelling of further chip formation once the chip makes contact with the cut surface of the workpiece.

Simulated temperature distribution
In all simulations, the workpiece reaches steady-state temperatures near the cutting zone and in Fig. 12, the predicted temperature distributions in the workpiece is shown at the final stage of the simulations with feed rate 0.18 mm/rev. As shown in Fig. 12a for the central insert cutting simulation, the temperature along the cutting zone close to the hole centre is about 250 • C and, as expected, it increases along the radial direction, i.e. as the local cutting speed increases. The temperature increases rapidly along the radial direction and about 5 mm from the hole centre, the cutting zone temperature reaches about 750 • C. It can  be noted that this is the radial distance from the hole centre where the mixed central/peripheral zone starts as shown in Fig. 5. Within the mixed zone, in which the nominal chip thickness gradually decreases towards zero, the temperature is fairly constant.
As shown in Fig. 12b for the peripheral cutting simulations, the workpiece temperature increases monotonically along the cutting zone from about 650 • C to about 850 • C. At the peripheral corner edge, the temperature reaches about 950 • C. The higher temperature at the corner edge can be partly explained by the somewhat higher cutting velocity but it is assumed that this is mainly explained by the more complex cutting conditions in this region. Along the corner edge, the nominal chip thickness decreases to zero and the plastic deformation and friction processes are different which results in a larger heat generation and different temperature distribution in the insert, chip and workpiece [39].
In Fig. 13a, the temperature distribution in the cut surface and the chip at the final stage of the central insert cutting simulations are compared for various feed rates. As shown, the temperature distributions are virtually the same for all feed rates except for at the outer rim of the chip where the temperature is substantially higher for the lowest feed rate 0.13 mm/rev. A possible explanation for this observation is that the finite element mesh cannot resolve the very small chip thickness in this region. As shown in Fig. 13b, also in the peripheral cutting simulations, the temperature distributions in the cut surface and the chip are virtually the same for all feed rates. It can therefore be concluded that the heat generation and the heat fluxes to the chip and workpiece are independent of the feed rate. Figures 14, 15, and 16 visualize the simulated temperature history and distribution in the inserts. In Fig. 14, the temperature history at nodes along the coordinate s as shown in Fig. 5 is visualized for the feed rate 0.18 mm/rev. During the remainder of this section, s indicates the coordinate at a path taken along the cutting edge through the hot spot temperature areas. That is, selected nodes along the path may be located on the flank face, rake face or on the edge radius depending on where the hot spot temperature is located for the specific s. Even though steady-state temperatures are not completely attained, the rate of temperature rise has slowed down considerably. Therefore, the insert temperatures at the final stage of the simulations are considered to approximate the steady-state temperature distributions. In Fig. 15a, the simulated central insert temperature at the final stage of the simulation with feed rate 0.18 mm/rev is visualized. Along the entire cutting edge, the hot spot temperature is distributed between the edge radius and a small part of the flank face. From the innermost corner edge and halfway to the outer corner edge, the hot spot temperature increases from ∼200 • C to ∼700 • C and, along the remainder of the cutting edge, the hot spot temperature is fairly constant. As shown in Fig. 15b, the hot spot temperature distribution along the cutting edge is similar for all feed rates. A comparison with Fig. 5 reveals that the highest hot spot temperature of ∼700 • C is reached at s ≈ 5 mm which coincides with the start of the mixed central/peripheral zone.
Within the mixed zone, the temperature is fairly constant. Only the temperature close to the outer corner edge increases with the feed rate.
For the peripheral insert, see Fig. 16a, the temperature increases monotonically along s from ∼450 • C to ∼950 • C at the start of the outer nose radius. In the mixed central/ peripheral zone, i.e. for s ≲ 3 mm, the hot spot temperature is located similarly as for the central insert, i.e. along the cutting edge radius. However, for s ≳ 3 mm, the hot spot temperature moves gradually more to the rake face of the insert. This can be explained by the increasing cutting speed along s as well as the larger chip thickness for s ≳ 3 mm. As shown in Fig. 16b, the temperature on the edge radius and flank face increases rapidly along s in the mixed central/ peripheral zone, whereas the hotspot temperature zone for s ≳ 3 mm is more homogeneously distributed along the rake face. Lower feed rates show similar hot spot temperature distributions. For a feed rate of 0.13 mm/rev, no temperature rise is recorded close to the inner corner radius. In this region, the chip thickness is very small and, for the lowest feed rate, the mesh discretization is too coarse to detect contact between the insert and workpiece.
The hot spot temperature distributions along the cutting edge of the central and peripheral insert clearly show that contact conditions and friction mechanisms act differently on the two inserts. Along the central cutting edge, the friction mechanism is more pronounced on the edge radius and the flank face. However, along the peripheral cutting edge, contact friction on the rake face is more significant. This is in-line with the discussion in [40] regarding wear mechanisms in conventional machining which outlines that the cutting speed has been shown to determine which friction mechanism is dominant, i.e. flank face friction is more pronounced at low-speed machining whereas rake face friction is more significant in high speed machining.

Conclusion
In this paper, experimental and FE-modelling investigations on an indexable drilling process with various feed rates are presented. Two CEL models are used to model each insert's cutting process individually and the torque and thrust forces from the simulations are superposed. Mesh convergence was demonstrated and mass scaling factors were selected for each simulation such that the drilling torque, thrust force, temperature fields and chip geometries were virtually the same as predicted in the corresponding simulations performed without mass scaling. The resulting drilling torque and thrust force at various feed rates were compared to their experimental counterparts. An excellent agreement of the drilling torque is obtained for all simulated feed rates where the maximum error is about 3%. A reasonable agreement is obtained for the simulated thrust forces versus the experimental ones where the maximum error is about 20%. Opposed to the experimentally measured thrust forces, the simulated thrust force is not significantly influenced by the feed rate. The effects of the friction coefficient and the frictional shear stress limit on the simulated torque and thrust force were investigated and it is concluded that the missing thrust force cannot be addressed by only changing the friction model parameters. Furthermore the simulated chip geometries were compared to the experimental ones. The simulated chips produced by the central insert showed severe discrepancy with the experimental ones that had a three times larger radius of curvature. In addition, the simulated temperature fields were analysed and it was shown that the temperature distributions in the workpiece, chip and insert are virtually insensitive to the feed rate within the tested range.
Author contributions All authors contributed to the study conception and design. The numerical work was performed by Daniel Svensson. Preparation of the experiments were performed by Tobias Andersson and Andreas Andersson Lassila. The experiments and data collection were mainly performed by Andreas Andersson Lassila. The first draft of the manuscript was mainly written by Daniel Svensson but parts was also written by Tobias Andersson. All authors read and commented on previous versions of the manuscript. All authors read and approved the final manuscript.
Funding Open access funding provided by University of Skövde. This work was supported financially by the Swedish Knowledge Foundation through the project CoSim (dnr: 20160298).

Conflicts of interest
The authors have no relevant financial or nonfinancial interests to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.