Application of mechanistic force models to features of arbitrary geometry at low material removal rate

This paper presents a workpiece discretization method to apply existing cutting force models to predict the forces generated during low material removal rate robotic machining operations of features with arbitrary geometry. Two machining operations along a straight edge which are modelled using this feature discretization method are shown, a chamfer pass on a sharp corner and the removal of a trapezoidal cross section. The workpiece features are measured using a high-resolution laser profile scanner to obtain the volume of the features to be removed. The identified features are discretized into rectangular sections such that the cutting force models can be applied to predict the cutting forces. A linear and an exponential mechanistic model which relate tool immersion and feed rate to the cutting force are applied to the scanned workpiece features. The linear and nonlinear models show good agreement with the measured data, with the exception that the linear model occasionally over predicts the forces depending on the radial depth of cut.


m
Generic axis N Number of teeth on discretized cutter r Tool radius t, r, a Tangential, radial, and axial axes for orthogonal cutting element V f eed Linear feed velocity x, y, z Cartesian coordinates or axes δ Tool helix angle ψ Angular offset of segments in discretized cutter φ Angular position of cutting tool φ st , φ ex Start and exit immersion angles P i Single profile of point cloud aligned normal to toolpath V i Set of vertices defining approximated chamfer area v i,a , v i,b , v i,c Vertices of the approximated chamfer area i Index adoption in industrial machining applications. Compared to traditional CNC machining, robotic manipulators offer larger workspaces, smaller physical footprint and increased dexterity at lower cost [1]. These advantages come at the cost of greatly increased compliance compared to traditional machine tools. While the cutting forces generated during bulk material removal operations may exceed the capacity of typical robotic manipulators, low material removal rate (MRR), and therefore low contact force applications such as deburring and polishing present a potentially useful application for robotic manipulators [1,2]. The compliance of the robotic manipulators makes them much less robust to the disturbances introduced by the stochastic geometry of the burrs. A number robotic de-burring strategies have been proposed in the literature that seek to improve the robustness of the robotic manipulator with some form of force or hybrid force-position control strategy [1,[3][4][5][6]. More recent efforts have included adaptive controllers that change the desired tool position and contact force in real time [7]. While some form of force control strategy is necessary to maintain stable contact with the workpiece [4], these methods of control are highly dependant on how accurately the forces can be predicted as a function of process parameters. Therefore, a process model with more accurate cutting force predictions can also be employed to improve the robustness of the system's overall control. This paper builds upon the authors' previous work [8] which showed that mechanistic force models could be applied to estimate the cutting forces generated through the machining of low volumes of aluminum. However, in the previous work, these models were only shown to be valid for machining rectangular cross-sections. For many low MRR operations such as chamfering and deburring, the cross-section of material being removed will be irregular and as a result the models as presented in Miller et al. [8] will be unable to estimate the forces. Thus, the goal of the current study is to present a novel method to apply the mechanistic cutting force models to features of arbitrary measured geometry. A cutting force model is required to convert a desired depth of cut or surface finish into a desired contact force for the force-feedback controller to maintain. The accuracy of this force model is critical to achieving the desired depth of cut, as errors is the commanded cutting force can lead to unacceptably deep or shallow cuts or, in the worst case, damage to the tooling or robot. This force modelling can be accomplished either implicitly or explicitly. For models with an implicit cutting force model, a range of acceptable cutting forces is collected experimentally and used to tune the control parameters to achieve a desired depth of cut and surface finish [3,9,10]. The implicit method is the simplest to implement, however, it requires the model to be retrained for every new workpiece tool combination and each new desired finish. The examples in the literature that do model the cutting force explicitly typically use a simple analytical model, for example in Her and Kazerooni [4] the tangential cutting force is directly proportional to the material removal rate. In Hsu and Fu [5] the contact force is modelled as an impedance using the relative position of the tool and edge of the part. For finishing operations with low material removal rate, such as deburring, where the goal is to remove a highly variable feature, the uncertainty in the feature geometry can contribute significantly to the cutting force, compared to traditional machining operations where the relative cutting force contributions of the uncertainty in the workpiece dimensions are typically small. Therefore, in order to accurately predict the machining forces generated, the features need to be measured directly. Obtaining an accurate measure of the surface can be achieved in a variety of ways, modern laser scanning technology in particular offers a fast and accurate method to measure burrs directly [11].
While few explicit deburring force models are present in the literature, there are a number of strategies for predicting the cutting forces generated by conventional milling operations. These existing models can be roughly divided into four categories: analytical, finite element methods (FEM), artificial intelligence and mechanistic models. Analytical models use the theories of solid mechanics to build a physically motivated model to predict the forces generated by cutting operations. Finite element models represent the workpiece as a mesh of discrete elements that are related together through a simple analytical relationship [12][13][14]. Both analytical and FEM models are a poor fit for this application as they require detailed knowledge of workpiece material properties that are difficult to estimate accurately, additionally, they are typically very computationally expensive to evaluate. Artificial intelligence and machine learning methods are increasingly common in the literature, however, these types of models are difficult to implement and require specialized knowledge to train and maintain which would limit its applications in industry as it would require additional training or hiring of staff. In addition, artificial intelligence models have a larger number of parameters and therefore require much larger training sets. This leaves mechanistic models as the best candidate for this application.
Mechanistic models use empirical parameters to tune a physically motivated basis function [15]. The empirical parameters capture the effects of both tool geometry and workpiece material properties and must be re-computed for each new combination of tool and workpiece material. Among the most commonly found in the literature are the mechanistic models proposed by Altintas and Keinzel. The Altintas model proposes a linear relationship between the instantaneous chip thickness h and cutting forces F along the axis of interest m, such that, where F t,r,a are the instantaneous forces in the tangential, radial and axial direction respectively, h is the instantaneous uncut chip thickness and a is the axial depth of cut, K (t,r,a)c are the empirical parameters that represent the forces generated by the mechanical shearing of the metal and K (t,r,a)e are empirical parameters that represent the forces generated by all sources of friction independent of the cutting action such as ploughing and edge rubbing [16]. Linear mechanistic models have shown good results when applied to discretized workpieces, as shown by the work of Schnös et al. [17] and Meirs et al. [18]. In both papers a linear mechanistic model is applied to a voxelized approximation of a known CAD-CAM tool path to accurately estimate the machining forces generated. Another popular mechanistic model in the literature is the exponential model proposed by Keinzel [19]. This exponential model proposes a relationship between the cutting forces along a given axis m, and the uncut chip thickness h such that where K n and β are the lumped empirical coefficients. The Keinzel model is less commonly used due in part to its nonlinearity. However, with current computing technology, the model can be evaluated in real time and the few papers that have used the model showed results that were comparable to or surpassed the accuracy of the Altintas model [19][20][21]. This paper is divided into six sections. The first section is a brief introduction to robotic deburring and cutting force modelling. The following section demonstrates the method through which the mechanistic force models shown in Eqs. (1) and (2) can be applied to a helical end mill. Section 3 details the method through which the mechanistic models can be applied to scanned features of arbitrary geometry. Section 4 presents the equipment and procedures used to validate the force model. The final two sections present the results of the experiments and offer concluding remarks.

Model form for helical cutting
The models as presented in Eqs. (1) and (2) are valid for individual cutting elements whose edge is orthogonal to the direction of cut, in order to apply them to milling operations they need to be integrated along the helical cutting edge of the end mill. The derivation presented in this section is adapted from Altintas et al. [16] and Adem et al. [19]. Since the models are developed for orthogonal cuts the fist step is to discretize the helical tool. The discretization of a 2-flute helical end mill in Fig. 1 is accomplished by first dividing the tool into j disk shaped cutting elements of height dz with i teeth, where i is the number of flutes on the tool. To approximate the helix, each disk is offset from the previous by an angle where where δ is the helix angle of the tool and D is the tool diameter. The instantaneous chip thickness h i,j for i th tooth of the j th cutting element is expressed as, where f is the feed per tooth in mm rev × tooth and φ i,j is the immersion angle of the i th tooth of the j th cutting element is given by where φ is the immersion angle of the first tooth of the first cutting element, N is the number of teeth per cutting element. It should be noted that Eq. (4) is a simplified model for the chip thickness that ignores cutter run-out and tool deflection. Substituting h = h i,j into Eqs. (1) and (2) for the linear model yields The relation of each element to the next [8] and for the exponential model where dF m j,i is the force in the m-axis for the i th tooth of the j th cutting element. In order to model the forces generated by the whole tool, the forces from each cutting element are transformed from their local (t, r, a) frame into the (x, y, z) tool frame with the transformation matrix ⎡ where the X-axis is along the direction of cut, the Y-axis is normal to the cut and the Z-axis is along the tool axis. The transformation matrix also contains a reflection about the X-axis and Y-axis to align with existing conventions for milling forces. Since the cutting elements now share a common coordinate frame, the forces can be integrated along the Z-axis where F m is the instantaneous force in the m direction as a function of the immersion angle φ and the feed per tooth where m ∈ x, y, z. The axial integration can be done analytically or numerically depending on the geometry of the cutting tool. In order to obtain the average forceF m , the instantaneous force is integrated over one full revolution of the cutter [16]. Figure 2 shows the definition of the limits of integration φ st and φ ex which are the immersion angles where the cutter enters and leaves the workpiece respectively. With the limits Fig. 2 Diagram showing the calculation of immersion angles for up milling operation [8] of integration defined, the mean force per revolutionF for the linear model which has an analytical solution is given as and the exponential model which must be integrated numericallȳ The force acting on the robot at each time-step can then be estimated by inputting the current depth of cut a and feed per tooth f into Eqs. (10)-(12) for the linear model and Eqs. (13)-(15) for the exponential model. The empirical cutting, frictional and lumped parameters of the models (K rc , K re , K tc , K te , K ac , K ae , K t , K r , K a , β) were identified using the simplex search method outlined in Adem et al. [19]. The training method consists of several full immersion tests conducted at a variety of feed rates, the parameters are then fit to the collected data using a simplex search method to minimize the error between the model output and the measured data. A more detailed account of the training process is available in [8], the results of which are summarized in Table 1 for convenience. While these models provide an accurate force estimate they function under the assumption that the cross section of the material being removed is rectangular. Therefore, in order to apply these models to the features of arbitrary geometry, the features themselves need to be discretized into rectangular segments.

Feature discretization
The following is an explanation of the process by which the identified features can be discretized for use with the cutting force models. Figure 3a shows an example of an arbitrary cross section in the world frame and Fig. 3(b) shows the first step of the process which is to transform the cross section into the tool frame. This transformation is not strictly necessary as a first step; however, the forces do need to be transformed into the tool frame to use as controller inputs and it makes the process easier to visualize.
Once the feature has been aligned to the tool-frame, the feature can then be divided up into n discrete segments of height dZ. For each dZ element the minimum and maximum immersion angles, φ st and φ ex , can then be found by determining the intersection of the segment with the tool envelope. Figure 4 shows the example profile discretized and indicates the three cases used to calculate the immersion angles: Case 1, full immersion; Case 2, partial immersion; and Case 3, overhangs.

Case 1: Full immersion
The tool is fully immersed in the feature, therefore the immersion angles for up-milling are simply given as and for down-milling Case 2: Partial immersion The feature projects only part way into the tool envelope. In this case the immersion angles for up-milling are given by and for down milling, where L seg is the length of the segment within the tool envelope and r is the tool radius.
Case 3: Overhang For the case of an overhanging feature the immersion angles for up-milling are given as and for down milling operations, Once the immersion angles are computed for each segment of the discretized feature, the mean cutting force per revolution, acting on the robot can be calculated bȳ where n is the number of discrete segments in the discretized feature andF m,k (φ st,k , φ ex,k ) is the mean force contribution of the k th segment of the discretized feature as given by Eqs. (10)- (15). With the discretization process in place, machining experiments can be conducted to validate the process. Figure 5 shows the configuration of the robotic manipulator used to collect the machining data used in this paper. The spindle assembly consisted of an NSK NR-3060S spindle powered by a 350-W brushless motor, (NSK EM 3060 J).

Robotic research platform
The spindle assembly was mounted to the DENSO VS-6556W 6-DOF robotic manipulator. An ATI IP60 Gamma force torque sensor was installed between the spindle and the last joint of the robot in order to measure the contact

Experimental procedure
The cutting force data was obtained by machining two features in an Aluminum 6061-T6 coupon. Figure 6 shows a cross section of the material removed during the two operations: Feature 1 is a standard 45 • chamfer with a base width of 1 mm; Feature 2 is created by repeating the tool path on the same coupon with a 0.2-mm Z-axis offset, removing a trapezoidal cross-section of material. The geometry of Feature 1 was selected to represent a typical operation to remove sharp edges. The geometry of Feature 2 was selected to represent a lower MRR finishing operation like deburring. Figure 5 shows the robot carrying out a chamfering operation.These two machining operations were repeated with increasing feed rate to gather a range of machining forces. Table 2 collects the feed rates used for all 16 tests, for each feed rate both a Feature 1 and Feature 2 machining operation were conducted. The feed rates used in the experiment were selected by determining the upper feed rate limit for the robotic platform by conducting tests at progressively higher feed rates until the tool would stall. Once the upper limit was determined a lower limit of 2 mm/s was selected to avoid overly slow tests. The resulting parameter space between the upper and lower limit were then sampled. Once the machining operations were  carried out, additional processing of the laser data was required before it could be discretized.

Feature scanning and identification
Before and after each machining operation, the part was scanned using a Micro-Epsilon ScanCONTROL 2950 profile scanner. Figure 7 shows a sample profile scan taken by the laser before and after each machining pass, with Feature 1 being the difference between the uncut edge and the chamfered edge and Feature 2 being the difference between the chamfered edge and the re-chamfered edge. To reduce computation time, the material removed from each profile P 1 , P 2 , ...P n is approximated as a triangle, with the trapezoidal geometry of Feature 2 defined as the difference between two triangles as shown in Fig. 8. The triangular approximation of the chamfer area is created through the following procedure, whose key features and order of operations are also presented in Fig. 9. As a summary, after the top face is identified from the shadow cast by the laser, the front face can then be located with respect to the top face. The bisector of both the top and front  face is then used to locate the cut face. The triangle defined by the three faces is then used to approximate the chamfer area.
First the leftmost edge of the top surface is identified by the sharp edge created by the laser's shadow, labelled 1 on Fig. 9. The top face is identified by including all points within 3 mm that fall within the standard deviation of the laser, given by the manufacturer as 0.05 mm. In this way, outlier points and points belonging to the cut face are  Fig. 9.
The first point of the front face, shown as 3 on Fig. 9, is located by choosing a point sufficiently far from the last point in the top face so as not to contain the cut face. For the data presented here the front face was defined as all points between −7.5 and −6.5 mm as seen in Fig. 7. A linear best fit line is then created for the front face from the identified points, line 4 on Fig. 9. The approximate centre of the cut face, labelled 5 on Fig. 9, is defined as the point closest to the line bisecting the lines fit to the front and top face. From the approximate centre, points on either side are added until the length of the group of points is equal to 25% of the commanded chamfer width. A linear equation is then fit to the identified points, represented by line 6 on Fig. 9.
The triangle described by the three identified trend-lines define the approximated chamfer area. The coordinates of the vertices of the identified triangle v i,a , v i,b , v i,c are then stored as a single set V i = {v i,a , v i,b , v i,c } to be used to calculate the cutting forces. A flowchart of the chamfer area approximation process is presented in Fig. 10. For Feature 2 the trapezoidal area removed is approximated by first performing the procedure outlined above to obtain a triangular approximation. Then the triangular approximation is truncated along the surface of the chamfered edge as seen in Fig. 8.
This process is repeated for every profile P 1 , P 2 , ...P n generating a set of vertices V 1 , V 2 , ....V n at every 400 µm along the path. As an example, Fig. 11 plots the identified feature area against time for the chamfer, Feature 1, test conducted at 10 mm/s. The ≈ 0.1 mm deviation noticed between 0 and 3 s is caused by a slight warping in the test coupon, bringing one corner of the coupon closer to the robot path, resulting in a slightly increased depth of cut.
To begin approximating the cutting forces, the method to discretize the scanned features described in Section 3 are applied to the approximated features. The approximated feature of each profile is discretized into rectangular  (15) to obtain the mean force per cutter revolution generated by each scanned profile along the part. Figure 12 presents a flowchart of the process used to generate force estimates from the approximated chamfer areas. It should be noted that this procedure is only valid for the 45 • chamfers discussed in this work, for more complex features a more robust feature identification method will need to be  developed. The results of this process are collected in the following section.

Results
In this section, the results of applying the existing mechanistic models to the features discretized following the process described in the previous two sections. In Section 5.1 the results of a single Feature 1 machining test are presented. Section 5.2 discusses the performance of both models at feed rates below 10 μm rev tooth . Following that, in Section 5.3 the behaviour of the Y-axis forces across both models and both features is discussed, in particular FPT for the machining of Feature 1 (chamfer). Error bars at 1 σ . Xaxis values shifted slightly for readability, exact values can be found in Table 2 the linear model's highly variable performance with Y-axis force estimation. The mean error valuesē presented in this section are calculated bȳ where F pred is the force predicted by the model, F mes is the measured force and N is the total number of data points. Table 3 collects the mean error for the two models investigated. For all Feature 1 (chamfer) tests Fig. 13 plots the mean force predicted by the linear model and the mean measured forces against FPT. Figure 14 plots the mean force predicted by the exponential model and the measured mean force against FPT.  Table 2 For the Feature 2 (trapezoid) tests Fig. 15 plots the mean force predicted by the linear model and the mean measured forces against the measured FPT. Figure 16 plots the mean force predicted by the exponential model and the measured mean force against the measured FPT.
For the Feature 1 (chamfer) trials, the linear model shows very poor agreement with the Y-axis forces with a mean error two orders of magnitude larger than that of the exponential model. However, for the trapezoidal tests the linear model outperforms the exponential model by a factor of 10 which. This discrepancy in the Y-axis is discussed in more detail in Section 5.3. To begin examining the results in detail it is useful to look at the results of a single test. Figure 17 plots the forces predicted by the models as dashed lines and the measured force as a solid line plotted against time for the Feature 1 (chamfer) pass conducted at 10 mm/s. The force estimates in Fig. 17 are generated by the identified features, whose approximate areas are plotted in Fig. 11. Looking at the X-axis forces in Fig. 17, the FPT for the machining of Feature 2 (trapezoid). Error bars at 1 σ . Xaxis values shifted slightly for readability, exact values can be found in Table 2 Fig. 16 Measured mean cutting force and exponential model prediction vs. FPT for the machining of Feature 2 (trapezoid). Error bars at 1 σ . X-axis values shifted slightly for readability, exact values can be found in Table 2 forces predicted by both models appear to follow the same general trend as the measured forces. From 0 to 3 s the Xaxis forces decrease slightly from 1 before levelling off at approximately 0.75 N. The identified feature area plotted in Fig. 11 follows a similar trend, the area decreases gradually over the first 3 s before levelling off. This qualitative result indicates that by including the specific geometry of the workpiece features, the models are better able to track the measured forces.

individual machining operation example
The variability in the model predictions in Fig. 17 are due to variance in the size of the identified features which correspond to the large spikes in the identified feature area that appear in Fig. 11. The identified areas of Feature 1 have a standard deviation of 0.03 mm 2 and the identified features of Feature 2 have a standard deviation of 0.032 mm 2 . These deviations in the identified areas are created by noise in the laser data caused by the surface finish of the top and cut faces scattering the laser light. The noise in the laser data can cause the best fit lines created from those faces (lines 2 and 6 on Fig. 9) to be misaligned from their true position and lead to large errors in the estimated feature geometry.  Figures 13,14,15 and 16 show the models' performances degrading at feed-rates below 10 μm rev tooth . Both models see a 200% mean increase in the mean X-axis error below 10 μm rev tooth across both features. The performance loss is due in part to the to the low cutting forces and subsequently reduced signal to noise ratio at the low feed rates.

Model performance loss at low chip thickness
The performance loss at low feed rates could also be explained by micro scale effects that manifest at low chip thickness. At chip thickness in the micron range the elastic recovery of the cut surface and the "ploughing" or plastic deformation of the chip start to have significant effects on the measured force [22,23]. The poor low feed rate performance could indicate that a mechanistic model that takes those micro scale machining effects into consideration might be necessary for very low FPT operations. For potential industrial applications of these models, this is less of an issue as the goal will be to machine any given edge with the largest possible FPT to maximize throughput. Figure 17 demonstrates a consistent issue with the linear model: while it is able to accurately track the X-axis forces, it dramatically over predicts the Y-axis (normal force) generated. The Y-axis forces plotted in Fig. 13 show the linear model over-estimating the Y-axis forces by nearly an order of magnitude.

Normal force estimation errors
The measured Y-axis forces, however, are typically around 0.1 N in magnitude which significantly reduces the signal to noise ratio compared to the X-axis forces, especially at feed rates above 10 μm rev tooth . The R 2 values collected in Table 4 show that the Y-axis forces are much less strongly correlated to the FPT than the X-axis forces in certain cases. A weak correlation to the feed rate and thus material removal rate could make it unsuitable to include as part of an online process model. The goal of the deburring force model is to return the desired contact forces to supply to a force controller to maintain a specific depth of cut, if a given Y-axis contact force can be obtained from a range of possible depths of cut then it would cease to function as a useful control input.
Considering these two factors, it may be the case that only the X-axis force can be reliably correlated to MRR at the feed rates tested using the linear model. Therefore, the linear model's poor Y-axis performance could be irrelevant. The results of Feature 2 in Table 4, however, indicate that in certain conditions the Y-axis forces do show a strong correlation to the FPT. The improved Y-axis force correlation for Feature 2 is also reflected in the Y-axis performance of the linear model, with the linear model's Y-axis error improving by nearly two orders of magnitude between Feature 1 and 2.
Another possible cause for the improved performance of the linear model between Features 1 and 2 is the specific geometry of Feature 2. To observe the effect of feature geometry on the force predictions, Fig. 18 plots the mean force per cutter revolution predicted by the linear and exponential models plotted against the tool immersion percentage at a 1-mm axial depth of cut and a FPT of 12 μm rev tooth . The Y-axis forces on Fig. 18 show that, at the mean depth of cut for Feature 1, the linear model is predicting a much larger cutting force. However, at the mean depth of cut of Feature 2, where the exponential model has

Conclusion
In this paper a method for discretizing scanned work pieces to apply existing mechanistic force models to predict the forces generated through the removal of features of arbitrary geometry was presented. By measuring and discretizing the features being machined at each point along the path, the exponential model was able to predict the resultant cutting forces to within 0.05 N along the X-axis and 0.006 N along the Y-axis for Feature 1 and within 0.12 N along the X-axis and 0.04 N along the Y-axis for Feature 2. The results of the linear model were less conclusive. The linear model's X-axis error was similar to the exponential model's with a mean X-axis force error of 0.24 N for Feature 1 and 0.12 N for Feature 2. For Feature 1, the linear model over predicted the Y-axis forces with a mean error two orders of magnitude larger than the exponential model's; however, for the Yaxis forces of Feature 2, the linear model outperformed the exponential model by an order of magnitude. Additionally, both models saw decreased performance at feed rates below approximately 8 μm rev tooth . This loss of performance at these low feed rates is likely due to the micro scale machining effects that start to become apparent at extremely small chip sizes. Availability of data and material Raw data is protected under IP agreements associated with CRDPJ 514258-17.
Code availability Code base is protected under IP agreements associated with CRDPJ 514258-17.

Declarations
Ethics approval No ethical approval required for the experiments conducted, as no human or animal participants involved.

Consent to participate
No consent required as no participants involved.

Consent for publication
No consent required as no participants involved.

Conflict of interest
The authors declare no competing interests.