Generic Mathematical Model for Efficient Milling Process Simulation

The current challenge in metal cutting models is to estimate cutting forces in order to achieve a more accurate and efficient machining process simulation and optimization system.This paper presents an efficientmathematical model for process simulation to evaluate the cutting action with variable part geometries of helical cutters and predict the cutting forces involved in the process. The objective of this paper has been twofold: to improve both the accuracy and computational efficiency of the algorithm for cutting force estimation in peripheral milling. Runout effect and the real tool tooth trajectory are taken into account to determine the instantaneous position of the cutting flute. An expression of average chip thickness for the engaged flute in the cut is derived for cutting force calculations resulting in a more efficient process simulation method in comparison with previous models. It provides an alternative to other studies in scientific literature commonly based on numerical integration. Experiments were carried out to verify the validity of the proposed method.


Introduction
An ever more demanding machining industry requires better accuracy and higher productivity.The major aim is to produce components in an efficient way from the very beginning and to avoid delays in production by performing costly machining tests; however, the reduction of setup times in production sometimes becomes a difficult endeavour.
Machining processes can be improved by strengthening the capacity to predict cutting process performance, which is often needed for automation or optimization of machining processes.In this sense, a cutting process estimation model should be the means of real time identification of the actual machining process in order to analyse the effect of any change on the variables involved in the process.
Milling is often the main cutting operation in the machining of dies and moulds where complex shapes with varying slopes, corners, thin walls, and so forth have to be machined.The cutting tool follows a trochoidal movement where the cutting edges engage in and disengage from the material periodically, resulting in an interrupted cutting process.Unlike other machining processes, milling has specific characteristics which make process modelling more difficult.One of these is the geometric characterization of the cutting edge.The variety of tool shapes is too wide to find a generic geometric model which properly represents milling tools.A classical approach in literature has been to develop milling models for each cutter shape.Moreover, cutting edge geometry and the preparation technology used to achieve this geometry have a strong influence on cutting forces and chip formation [1].A generalized geometric model for solid and inserted milling cutters was first presented by Altintas and Engin [2].This model has subsequently been used by other researchers like Merdol and Altintas [3], Wan et al. [4], and Soori et al. [5], among others.
As a consequence of the helix angle, there is a nonuniform distribution of chip thickness along the cutting edge.Every The relation between chip thicknesses and cutting forces is defined through a function named the specific cutting force coefficient or the specific cutting pressure.There are two main approaches available in order to evaluate these coefficients.Firstly, there is the orthogonal to oblique transformation method, as proposed by Armarego [6], where the coefficients are evaluated based on shear stress, shear angle, and friction angle.As this method required abundant orthogonal cutting tests, Wan et al. [7] designed an algorithm to determine shear stress, shear angle, and friction angle directly from a few milling tests.Secondly, the cutting force coefficients can be directly obtained from milling tests for a specific tool geometry and workpiece material.Both methods require experimental identification.
In some mechanistic models, the cutting action can be described separately by two force coefficients related to the shearing mechanism at the rake face and the ploughing mechanism at the nose and flank face of the cutting edge, as proposed by Junz Wang and Zheng [8] and Larue and Altintas [9].Wan et al. [10] proposed the inclusion of a third parameter related to the bottom edge effect of the cutting edge.
It is well known that the accuracy of the prediction for a given machining process simulation model depends not only on the mathematical model used to predict the physics of the machining process, but also on the accuracy of the experimentally determined cutting force coefficients.In this sense, the evaluation of chip thickness and cutting force coefficients are crucial for mechanistic modelling of cutting forces.
The function which best estimates the specific cutting force coefficients is the potential function which takes into account the size effect of the cutting edge.However, to obtain the total force acting on a flute, the solution of the integral is compromised by this nonlinear term of the integrand.In order to avoid this inconvenience, they are frequently estimated as constants or as linear functions of chip thickness.Initial research used average cutting coefficients obtained from measured average cutting forces, to predict cutting forces.Kline et al. [11] developed a cutting force model, where the specific cutting forces are correlated with the machining conditions through a second order polynomial model.Based on measured average cutting forces as well, Gradišek et al. [12] calibrated the cutting force coefficients for a general end mill, considering them constant.
In an attempt to determine force coefficients regardless of cutting conditions, instantaneous cutting force coefficients were introduced.Ko and Cho [13] expressed the instantaneous cutting force coefficients as a Weibull function of the uncut chip thickness.Wan et al. [14] proposed the cutting force coefficients as power functions of instantaneous chip thickness.These coefficients and the runout parameters were simultaneously obtained from a single milling test.Wang et al. [15] proposed an improved method to identify the instantaneous cutting force coefficients by using surface errors.A disadvantage of these models is the extensive testing that is required for the determination of the input parameters.
In the majority of cutting force models, it is common practice to digitize the engaged flute into small disc elements and evaluate the cutter-part engagement condition for each disc in order to determine the cutting forces acting on the flute.However, a densely digitized cutting flute may unnecessarily increase the computation time.
However, this does not occur when analytical models are used.Bhattacharyya et al. [16] presented closed form expressions for cutting forces, where the entry/exit angles of the cutting flute were expressed as Heaviside unit step function and Fourier series formulation.Although closed form expressions are used, it is very time consuming because it requires computing two Fourier series for both the leading and trailing points of the cutting flute.In their study about stability in milling, Jin et al. [17] pointed out the loss of accuracy when a discretized cutting flute is used to estimate the cutting forces; therefore they used an analytical method instead.
Another problem which arises from modelling is that although the mechanics of milling processes and chip generation are well established when they are oriented to steady-state cutting processes, cutter-workpiece intersection domains are not so simple to determine when transient states are issued.These are, however, the most common cutting operations in milling.
It can be concluded from this that there are still challenges in achieving computationally efficient but accurate modelling of cutting forces in machining processes.Accurate predictions usually lead to long simulation times which make them impractical for industry whereas when rough analysis is used, the approximation quality is of limited value [18].This paper presents an efficient mathematical model to evaluate the cutting action with variable part geometries of helical cutters and predict the cutting forces involved in the process.The objective of this paper has been twofold: to improve both the accuracy and computational efficiency of the algorithm for cutting force estimation in peripheral milling.
Compared with existing models and methods, the main contribution of this algorithm is that it enables us to simulate milling processes when cutting conditions change continuously, without the requirement of a long computing time.No such work has been found in literature.This is achieved by means of a cutting force estimation model based on closed form expressions, in an attempt to avoid numerical integration.This procedure also allows for the possibility of extending research topics into other areas of interest within the general framework of machining processes.

Cutting Condition Modelling of Milling Tools
The first step in machining process simulation is the identification of the intersection between the cutter and the workpiece in order to determine the material that the tool actually removes.The cutting edge-part engagement will lead to chip thickness variation and the resultant cutting force calculation.Since the designed algorithm intends to be of a general nature applicable to variable geometries, the axial and radial depth of cut need to be evaluated at discrete intervals.For a given helical flat-end mill, the position of any point on the cutting edge  in the cut is determined by an angle   , which is measured clockwise from the positive -axis in a normal plane to the tool axis.The corresponding instantaneous chip thickness ℎ  is the radial difference between two consecutive trajectories of subsequent cutting flutes and varies periodically with the position angle   : is the feed per tooth.
Figure 1 shows the action of one flute on a plane that represents the unfolded surface of cut.The cutting edge is a straight line inclined by the helix angle   and it moves from left to right.The angle   indicates the movement of the leading point from the beginning of its cutting action.In the figure, three different phases depict the interaction of the cutting flute with the workpiece: two transient phases at the entry/ exit of the cutter and the uniform cutting zone.
Following the scheme shown in Figure 1, the flute starts the cut from the side AB, the point A being the starting point and AB the programmed depth of cut   .As the flute advances, the length of the flute engaged in the cut increases to its maximum value.This is the transient phase at the entry of the cutting flute into the material.Its value depends on the flute engagement angle  pr which is the projection on -plane of the cutting flute in contact with the workpiece material: In the uniform cutting zone, the length of the cutting flute remains constant while it moves along it.Once the leading point reaches , which is defined by the exit angle   , the cutting edge length gradually decreases until the cutting flute is completely out of the cut.This is the transient phase at the exit of the cutting flute.
The surface of the cutter in contact with the workpiece is defined by the difference between the exit and entry angles, named the engagement angle of the cutter  en .Depending on whether  en <  pr or  en ≥  pr means that one cutting flute or more than one may be cutting at the same time.
In order to identify the intersection domain between the cutter and the workpiece, it is necessary to determine the starting point where the flute begins the cut and the point where the same flute exits from the cut.
For the proposed algorithm, the workpiece and cutter geometry variables such as the tool diameter , the helix angle   , and the number of flutes   , along with the machining parameters, are input data.It calculates the positions of the cutter axis considering the real tooth trajectory and in case its effect is noticeable on cutting forces, the cutter runout is considered as well.In this case, the variation of cutting point radius for parallel axis offset is expressed as follows: It is assumed that tool runout is mainly due to the excentricity of the tool-toolholder-spindle assembly, as described by Kline and DeVor [19].In this case, tool excentricity is characterized by means of two parameters that represent the value of the tool offset  and its angular position  with regard to a reference flute ( = 1).The path of motion of a point on the cutting edge which intersects the workpiece will be a trochoidal trajectory.Subsequently, the intersection between the cutting edge and the workpiece is calculated, determining the part of the cutting flute that actually engages with the material;  1 is the leading point and  2 is the trailing point of the engaged cutting flute.This leads to the chip load calculation and the corresponding cutting forces.
The relationship between the position of the entry/exit angles throughout the advance of the cutting edge ( en ) and the flute engagement angle will result in six different machining cases which cover any possibility of cutting operation.A coordinate system on the spindle of the machine tool will be the reference in order to position the cutting flutes and the workpiece boundary.In this classification, two main subgroups can be distinguished, depending on whether there is a uniform cutting zone or not.The latter consists only of transient states, without a uniform cutting zone.In these cases, once the flute has entered into the cut, it begins the exit before completing its entry.Figure 3 shows the pseudocode for the classification algorithm of the six possible cases.

Mathematical Model for Cutting Force Calculation
The cutting force on any point of the cutting edge is distributed nonuniformly along the cutting edge in contact with the workpiece since every point has different angular positions   and according to (1), they will cut different chip thicknesses.Figure 3 shows the nonuniform distribution of chip thickness for the engaged length of the cutting flute  at time .
According to the mechanistic model, the contribution of a differential element on the cutting edge to the total cutting force is where   (),   (), and   () are the cutting force coefficients for tangential, radial, and axial directions, respectively.The cutting force coefficients along with the chip thickness are the primary means to predict the cutting forces.The accuracy and quality of these calculations will depend on the analytical method used to calculate them.
Differential cutting forces are integrated along the engaged length of the previously determined cutting edge.The presence of any nonlinearity may condition the solution of the integrals, whether analytically or numerically.
Common practice is to solve it numerically by dividing the engaged length of the cutting flute into thin axially divided disc elements, as shown in Figure 4, and to evaluate the chip load for each disc.For disc  the tangential force is expressed as follows: Δ ,, () =   () ⋅ ℎ , () ⋅ Δ  () . ( The radial and axial forces have similar expressions.
The cutting forces along , , and  directions can be mathematically expressed as Subsequently, the total cutting force components at any tool rotation angle  can be evaluated by summing the forces acting on all discs and on all flutes of the cutter: ,, () ,  = , , .
The milling force estimation model proposed is based on the average chip thickness for the engaged flute.The average chip thickness is calculated for the flute engagement angle, according to the following expression: With this method, only one iteration is needed for every angular position of the cutting tool.In this sense, the proposed objective of efficiency is achieved, this method being much faster than those using numerical integration techniques without loss of accuracy.The proposed algorithm offers a computationally efficient means to simulate the cutting forces.In addition, this method can potentially be extended to other aspects of milling such as tool deflection, vibrations, or tool wear.Following the scheme presented in Figure 1, during the transition state at the entry of the cutting edge into the workpiece, the average chip thickness for the engaged flute is calculated as follows: As the cutting flute advances, the axial depth of cut in this zone increases steadily until the programmed value is achieved.It is calculated with the following equation: In the uniform cutting zone, where axial and radial depth of cut are kept constant, the average chip thickness is calculated as Once the cutting edge starts exiting from the workpiece, the engaged length of the flute decreases until it is completely out of the cut.In this phase, the average chip thickness is calculated as In a similar way as at entry, the axial depth of cut is calculated as follows: Having determined the flute engaged in the cut, cutting forces from ( 4) can be obtained by taking the integral along the active portion of the cutting edge.The solution is conditioned by the existence of any nonlinearity in the terms of the integrand.One of these terms is the cutting force coefficient and more specifically the function used to estimate it.In order to avoid this inconvenience, they are frequently estimated as constants or as linear functions of chip thickness.In these cases, the solution of the integral is easily found.In the algorithm proposed here and taking into account the size effect, a better fit is given by a potential function of the average chip thickness as defined in (8): where   ,   ,   ,   ,   , and   are the parameters empirically obtained by running a set of milling tests for a specific cutter-workpiece combination.Since the proposed methodology is oriented to variable cutting conditions and in order to avoid numerous fitting tests, the cutting force coefficients are modelled as varying parameters of the position angle, instead of being simply defined as constants or functions based on average cutting forces.
Figure 5 shows the evolution of the cutting force coefficients   and   with chip thickness.It is worth noting the size effect trend of the fitted values, where cutting force coefficients increase exponentially for small values of chip thickness.
In short, the proposed algorithm executes the following steps.
Step 1. Determine the nominal tool path based on the cutter radius and the geometry of the workpiece.
Step 2. Calculate the actual position of the cutter with (3).
Step 3. Calculate the entry/exit angles for cutting flute .
Step 4. Determine the case of cut, as shown in Figure 2.
Step 5. Calculate the average chip thickness for the engaged flute with (8).
Step 6. Calculate the cutting forces for cutting edge .
Step 7. Repeat Steps 3-6 until the cutting process is completed.
As was shown above, the cutting force calculation using the conventional method requires the calculation of the cutting forces for each individual disc, whereas the proposed method requires only one iteration.The greater the number of discs, the greater the difference between both methods.
A comparison between the conventional approach, based on the discretized cutting edge, and the proposed approach, based on the average chip thickness, was made and the results are shown in Figure 6. Figure 6(a) shows the average chip thickness for the engaged flute and the corresponding chip thickness for each disc element.The estimated cutting forces for both models are shown in Figure 6(b).As can be seen, their precision is very similar; however, the ( − 1) reduction of computing time when comparing the proposed model with conventional method is the main advantage of our model.
The accuracy of the proposed model for the cutting force calculation on one flute was shown and now the model can be applied to a more general milling tool and under variable conditions.Figure 7 shows a transient machining operation in milling.The cutter engages the workpiece and during its advancement, a change in radial depth of cut is encountered.This cutting operation can be divided into four different states; the first three are transient states until the uniform cutting is achieved in the fourth state.
Zone 1 represents the start of the machining cut where the flute engages progressively with the material.The engagement arc increases with the advance of the cutting tool.The entry angle is kept constant while the exit angle varies in each cutter revolution and the cutting force profile varies accordingly.Zone 2 begins when the flute encounters a change in radial depth of cut.In this state a special situation is observed.During some rotations of the cutter, the load on the cutting flute is not continuous.As shown in Figure 6, not only does the maximum in   and   vary but also the frequency of the cutting forces varies.This nonuniform load along the cutting flute causes an irregular distribution of the derived stresses on the different sections of the cutter, which may lead to cutter breakage.
In zone 3, the cutting flute withstands a continuous load again.The arc of engagement increases while the chip thickness at the exit of the cutter decreases progressively.The wave length in the abscissa axis slightly increases as well.Finally, zone 4 corresponds to uniform cutting.All parameters are kept constant at this stage.
As is shown in Figure 7, it is possible to estimate the cutting forces precisely in transient states by considering the action of the whole cutting edge.As the position of each   cutting edge involved in the cut is computed at every time interval, it is possible to analyse any small change which may have occurred during these transient states.

Experimental Verification
The proposed algorithm has been evaluated by carrying out a set of transient machining tests which are naturally encountered in milling, some of them being the entry/exit of the cutter into the workpiece, pocket machining, or changes in axial/radial depth of cut.
For this purpose a vertical CNC milling centre was used.Cutting forces were measured with a Kistler 9257B dynamometer, mounted between the workpiece and the machining table.A low-pass filter at a cut-off frequency of 100 Hz was used to facilitate an easier comparison between the measured cutting forces and the estimated ones.
A two-fluted cylindrical end mill with a diameter of 8 mm and helix angle of 30 ∘ was used for a pocket machining test.This experiment brings together two characteristics that point the difficulty of modelling the milling process: variation of the cutting conditions at entry and change in the feed direction of the cutter.Before the cutting tests were performed, cutter runout was measured in an offline tool presetting in the machine tool.The position of the cutting edge with respect to the rotational axis was measured.From this, every position of the cutting edge was estimated.
It is worth mentioning that although measurements capture the cutting force variations from the cutting test, there is a difficulty in matching the measured cutting forces with the estimated ones in this type of experiment.For each cutter revolution, the tool displacement equals the feed per tooth until the intersection between the cutting flute and the part occurs.
In order to compare the measured cutting forces with the predicted ones, it is essential to determine the angular position of the cutting flutes when the intersection occurs.If this is not the case, the simulations would be limited only to obtaining the maximum and minimum values of the cutting forces for varying immersion and tool positions.The cutting forces would not show any detail of the cuts which in fact does occur.
In a CNC machine tool, once the cutter starts linear displacement, it is already rotating.The linear displacement of the table is independent from the rotation of the cutting tool.It is not possible to know in advance the angular position of the cutting flutes at the time when the cutter initiates the displacement.
However, for cutting force simulation, it is necessary to determine the real tooth trajectory when the flute engages the part.This requires the synchronization of the feed rate and the rotation of the cutting tool.
To overcome this obstacle, a laser sensor is used to generate a reference pulse in each tool revolution, which is necessary in order to identify the reference flute and control the entry of each cutting edge.A sketch of the workpiece geometry and the tool paths are shown in Figure 8 along with the simulated and measured cutting forces.The waveform of the predicted cutting forces is consistent with the measured one both in shape and size.Note that although both cutting force waves have the same frequency, the origin is different.In the case of the measured cutting forces, the origin has been preserved from data acquisition.In this experiment, it is necessary to identify not only the actual position of the cutting flute but also the actual portion of the flute that first engages the corner.
The results of the experiment from Figure 8 show variation in the cutting conditions in two different areas which are highlighted in the graph.These areas correspond to the entry of the cutting flute in the corner and when the cutter changes its feed direction.A closer view of these two areas is shown in Figures 9 and 10, respectively.Both entry and exit angles vary in each cutter revolution and the cutting force profile varies accordingly.From Figure 9, it should be noted that during some rotations of the cutter, the load on the cutting flute is discontinuous.Not only does the maximum in   and   vary but the period of the cutting forces and in this case the frequency also vary.
As a result of the synchronization process, there are slight discrepancies between the predicted and the measured cutting forces in Figure 10.Nevertheless, the results show that the proposed method is an efficient way of simulating transient cases.

Conclusions
The main conclusions of this paper are the following.
(i) A new mathematical model for cutting force estimation in peripheral milling is proposed in this paper.Attention has been paid to two opposing aspects, the prediction accuracy and the computational efficiency of the proposed algorithm.Usually accurate predictions lead to long simulation times which make them impractical for industry while rough predictions lead to inaccuracies in cutting force simulation.It has been possible to find a trade-off solution with this proposed algorithm.
(ii) Since the real tooth trajectory for the cutter is taken into account, the entry/exit angles for the cutting flute are determined accurately.The identification of the intersection between the cutting flute and the workpiece allows us to determine the material that the tool actually removes.The simulation runs in a series of small time steps along a cutter revolution.At each time interval, the position of each cutting edge involved in the cut is computed so that any small change in part geometry can be detected along the tool path.This allows us to detect any variation in the period and frequency of the cutting forces.
(iii) A novel expression for the average chip thickness of the engaged flute is derived in order to avoid numerical integration.Closed form formulations are obtained for each force component resulting in a computationally efficient mathematical model when compared with numerical models which are mainly used in milling process research.
(iv) Compared with existing models and methods, the main contribution of this algorithm is that it enables us to simulate milling process when cutting conditions change continuously, without the requirement of a long computing time.
(v) According to the results obtained from simulations and cutting tests, the validity of the model is proven.The simulations and experimental results showed agreement when changes in radial depth of cut and in the feed direction occur.

2 Mathematical
Problems in Engineering cutting point on the flute has different chip loads and, in consequence, different cutting forces.

Figure 1 :
Figure 1: Different phases of the engaged cutting flute.

Figure 2 :Figure 3 :
Figure 2: (a) Pseudocode for the different cases for the engaged flute.(b) Example case 3.

Figure 4 :
Figure 4: Cutting edge discretization in  disc elements.

Figure 7 :
Figure 7: Schematic transient machining cut and simulated cutting forces.

Figure 9 :
Figure 9: A closer view of the cutting forces in Figure 8 at the entry of the cutter in the corner.

Figure 10 :
Figure 10: A closer view of the cutting forces in Figure 8 when the cutter changes feed direction.