Predicting regenerative chatter in milling with hardware-in-the-loop simulation using a dexel-based cutting model

This article proposes a new model for simulating the interaction between cutting process and machine tool in real-time. The purpose of the model is to be coupled with a real CNC (by using hardware-in-the-loop simulation) in order to consider process forces and to predict regenerative chatter vibrations during virtual commissioning. Therefore a dexel-based workpiece model with adaptive resolution is used for the computation of the chip thickness respectively the cutting forces based on the actual machine tool position and the machining progress on the workpiece. Several simulation experiments are performed to validate the model and to analyze its numerical limits, such as computational accuracy and efficiency. The capability of the model to predict chatter is proven by comparing the simulated critical depth of cut with an analytical solution of the stability lobes. Therefore the dynamics of the machine tool were approximated as a single degree of freedom (SDOF) oscillator. A concluding analysis of the real-time factor confirms the model’s ability to be integrated under hard real-time requirements and with cycle times of just a few milliseconds which are typical of CNCs.


Introduction
The prediction of the machine tool behavior via simulation in consideration of the cutting process has been one of the most important research topics in production engineering within the last decades. The challenges of simulating the entire machine tool behavior are characterized by a diversity of different physical phenomena (e.g. flexible multibody systems, non-linear process forces, heat transfer) which have to be partially modeled in detail depending on the individual problem as described by Altintas et al. [1]. Especially process simulation places high demands on the utilized simulation tools. For the simulation-based prediction of process forces and process stability, e.g. chatter stability in milling, expert knowledge of the process itself and its simulation is required as shown by Brecher [2]. The prediction of the dynamical stability of cutting processes such as turning, milling and drilling processes has been investigated since the 1950s by Tlusty [3], Tobias [4] and Merrit [5]. Analytical approaches were introduced by Altintas in numerous publications. A good description of an analytical method for predicting the dynamical stability of milling processes is found in [6,7]. Numerical methods were presented and improved by Henninger [8] and Chanda [9].
The application of simulation models in connection with real control systems has become widely popular in many fields of research and development in the recent years. Starting with the hardware-in-the-loop-simulation for virtual commissioning of control systems by Pritschow [10], up to simulation-based control functionalities such as condition monitoring or model-predictive path planning algorithms by Verl [11]. Further results were presented by Röck using reduced structural dynamic models [12] and process machine interaction models [13]. In all cases, a time-deterministic and highly efficient computation of the machine behavior is required. The interaction between cutting process and machine tool has a decisive impact on the process results and cannot be neglected in a simulation used for a process-oriented set-up. The actual cutting forces depend highly on the velocity profiles generated by the computer numerical control (CNC) and vary significantly subject to the CAD/ 1 3 CAM process, the control parametrization and the chosen CNC system. This leads to different mechanical excitations of the controlled machine tool and, consequently, to different process results which can only be predicted via a control-coupled simulation in real-time. Figure 1 shows a hardware-in-the-loop simulation scheme for virtual commissioning of real CNC systems.
In this paper a new approach for modelling cutting processes for real-time simulation in order to consider process forces and predict process instability during virtual commissioning is presented. The main benefit of the proposed model over analytical approaches is that it can be used for virtual commissioning of a real CNC program with changing cutting conditions (e.g. transition from zero to full immersion when entering the workpiece, variable feed rates, variable tool penetration depths). In contrast, analytical approaches are only capable of calculating the stability for the static case (settled state) without considering changing cutting conditions (only full immersion with constant feed rate and constant penetration depth).
To make this possible, process forces under real dynamic cutting conditions have to be taken into account which in consequence requires a model that keeps track of the workpiece's current state. Hence a dexel-based workpiece model with adaptive resolution is developed for the computation of the chip thickness based on the actual machine tool position and the machining progress on the workpiece. The proposed approach in this paper is limited to planar milling processes but will be extended to geometrically more complex cutting processes in future work.

Dexel based cutting model for computation of the chip thickness
The geometry of the workpiece is approximated using a defined arrangement of so-called dexels. A single dexel represents a (bounded) ray passing through the workpiece and is geometrically described as a static segment of the straight line with ⇀ r b,i ∈ ℝ 2 being its start point and ⇀ r e,i ∈ ℝ 2 being its end point. Depending on the placement and the length of the individual dexels, differently shaped workpieces can be modeled. The workpiece model in this paper is limited to two dimensions for simulating planar milling processes. In this case, the workpiece is represented as a parallel arrangement of dexels on a plane and each of the tool's cutting edges is a point projected onto this plane.
The path of the j-th cutting edge ⇀ e j created by the tool rotation is approximated as a linear interpolation between its current position ⇀ e j,n at time t n and its previous position (1)  (Fig. 2). Thus, the chip thickness h j,n of the j-th cutting edge is the arithmetic average of all h ij,n over the cutting dexels with indices For material removal the start point of the dexel has to be updated to If the cutting edge path ⇀ p j,n is located between two dexels, h ij,n cannot be calculated directly by the intersection point ⇀ c ij,n . This problem occurs if the spindle speed or the time h ij,n step size of the simulation is small. For varying tool paths, feed rates and spindle speeds it is not possible to predict all cutting edge positions in advance to be able to initialize a sufficient number of dexels located at the right positions before the simulation starts. It is also not feasible to fill the whole workpiece with just enough uniformly spaced dexels without losing computational efficiency. The only way to solve this problem is to adapt the density of dexels during runtime by inserting a new dexel exactly on the current position ⇀ e j,n of the cutting edge if no dexel is available. This approach converges to an optimized distribution of dexels within the workpiece and simultaneously improves the approximation of the chip thickness especially in areas where high dexel density is needed. Figures 3 and 4 show the evolution of the number of dexels during the penetration of a tool with N = 10 teeth and a radius of r T = 0.01[m] in a rectangular workpiece with initially N D = 20 dexels. The number of dexels increases until the tool has completely penetrated into the workpiece (Fig. 4).

Cutting forces
The radial and tangential cutting forces acting on each cutting edge depend on the chip thickness h j,n at time t n and can be calculated by: with the radial and tangential cutting coefficients k r and k t and the depth of cut a [6]. After summing up all cutting edge forces and transforming them into xy-coordinates by using the rotation matrix Φ j,n of the j-th tooth we get the total forces acting on the tool (7) F rj,n = k r ah j,n F tj,n = k t ah j,n , with the tool radius r T . Figure 5 shows the acting forces on the cutting edges respectively on the tool center point.

Static force simulation
In the following various simulation results in time domain will be presented. It will be shown that based on the dexel model valid forces for milling processes can be computed. Moreover, the numerical influence of the simulation time step size Δt and the spindle speed Ω will be discussed. To validate the simulated milling process forces it is sufficient and more comprehensible to consider only one tooth ( Fig. 6) with the parameters in Table 1. It should be noted at this point that the main focus of this simulation was not on the analysis of a realworld process, but rather on validating the method. Therefore the values in Table 1 were chosen to represent a convenient parameterization for machining light metals (e.g. aluminium).
To validate the dexel-based cutting forces a comparison against the "exact" analytically calculated forces is performed. Thus the exact chip thickness for one tooth is given by with the feed rate v x , the time period for one spindle revolution rev and the spindle speed Ω. The unit step function g(Ωt) depends on the angular position of the cutting edge Ωt and determines whether the tooth is in or out of cut. The cutting force can then be calculated to (11) g(Ωt) = 1 ← −0.5 < Ωtmod(2 ) < +0.5 g(Ωt) = 0 ← otherwise. Comparing the results of one full revolution − < Ωtmod(2 ) < , the simulated dexel-based cutting forces F x,sim , F y,sim show a good correlation with the analytically calculated exact forces F x,ex , F y,ex . Figure 7 shows the results with different simulation time step sizes Δt. Figure 8 shows the influence of the spindle speed. The left picture demonstrates how a low spindle speed of Ω = 200 rev∕min (v x = 0.001 m∕s ) suppor ts a very high accuracy of the simulation results in consequence of the dexel adaptation. However, the right picture shows the result in case of high spindle speeds, e.g. Ω = 4000 rev∕min (v x = 0.02 m∕s ). The error comes from the gross path discretization which is directly coupled with the time step size and the spindle speed. Reducing the time step size gives more accurate results but at the same time worsens the computational efficiency. Figure 8 highlights the importance of knowing the accuracy of the model in dependence of the spindle speed. To quantify the error we use the absolute differences between the simulated and the exact forces and define the percentual maximum error and the percentual mean error Figure 9 shows the percentual maximum and mean error versus the spindle speed for different simulation step sizes with v x = 0.02 m∕s :   To complete the validation for the milling process it is necessary to compare simulation results with more than one cutting tooth. Therefore the scenario from Fig. 6 and Table 1  For milling processes with N = 4 (N = 6, N = 8, N = 10, … ) constant cutting forces occur if the tool is completely penetrated in the workpiece. In that case the error caused by a gross path discretization due to large spindle speeds or large simulation time step sizes becomes crucial. This can be seen in the middle picture of Fig. 11. These errors are not acceptable and the only way to avoid them is to reduce the time step size, as demonstrated in the right picture of Fig. 11.

Dynamic force simulation
In the previous sections, it was assumed that the machine tool is absolutely stiff and no vibrations occur. In the following the model of the tool has been extended to a single mass oscillator to consider the first natural frequency of structural dynamics.

Single degree of freedom (SDOF) vibration model
To validate the dynamic forces and the ability to predict the chatter phenomenon in milling a SDOF model which covers the tool vibrations in feed direction is used. It is described as a single mass oscillator with a spring-damper element between the tool center point and the tool holder excited by the cutting force (Fig. 12). The equation of motion can be given by where x denotes the tool position and F x the cutting force in feed direction computed by the dexel model. u represents the target position of the tool (commanded by the CNC), 0 the natural eigenfrequency, D the damping ratio and k 0 the stiffness of the spindle-tool system.
Because of the hard real-time requirements in a CNC-coupled hardware-in-the-loop simulation, a very efficient time integration method-the double shot semi-implicit Euler method-is used to solve Eq. (18). With x n = x 0 and v n =̇x 0 as the initial condition, the following integration scheme is executed for each time step t n+1 = t n + Δt: Ω, a, … ), The explicit scheme in Eq. (19) requires only one evaluation of F x per time step, which is the computationally most expensive part of the calculation. This leads to a very fast solution and enables the use of small simulation time step sizes Δt.

Transient results
The transient simulation shows the process forces while the tool penetrates into the workpiece. Based on the earlier introduced parameters (see Table 1), the extended parameter list in Table 2 is used in the following.
The simulation results are shown in Fig. 13 with a = 0.01[m] (process stable) and in Fig. 14 with a = 0.03[m] (process unstable). Because of the tool's symmetry with N = 10 the force F x becomes constant if the tool has completely penetrated into the workpiece and the process is stable (Fig. 13). In case of instability the force F x , and therefore also the tool position x , starts oscillating with the chatter frequency c ≈ 530 rad∕s which is close to the natural frequency 0 = 500 rad∕s of the SDOF model. In that case the tool jumps out of cut and the force periodically becomes zero (zero cut)-regenerative chatter occurs (Fig. 14).

Process stability lobes
The reliable prediction of chatter vibrations as shown in Fig. 14 is very important for the virtual commissioning of milling processes. The following results verify the predictability of stability lobes with the SDOF vibration model combined with the dexel-based cutting force model from Fig. 12. Again, the simulation scenarios shown in Figs. 13 and 14 are used, but this time with the parametrization from Table 2. The results are compared with the analytically calculated stability lobes based on the method presented by Altintas and Budak in [6]. The graph in Fig. 15 shows the depth of cut a over the spindle speed Ω with the analytically calculated critical depth of cut a lim (stability lobes) and stable/unstable simulation results based on the dexel model (Fig. 16). Figure 15 shows a very good correlation between the dexel-based approach and the analytical solution for spindle speeds between Ω = 200 rev∕min and Ω = 600 rev∕min . A detailed representation of the forces F x with a spindle speed of Ω = 600 rev∕min is demonstrated in Fig. 17.

Real-time capability
The real-time capability is represented by the real-time factor (RTF) which is defined as

Fig. 12 SDOF vibration model
where t sim denotes the simulation time interval and t CPU denotes the CPU time consumed to execute the simulation time interval. For real-time simulation, the RTF has to be less than 1. Figure 17 shows the RTF over the time for the simulation of the milling scenario shown in Fig. 12 with Ω = 2000[rev∕min] and Δt = 0.001[s] . The RTF is nearly constant over time. Only at the start of the simulation, when the tool touches and enters the workpiece, the computational effort is varying due to the adjustment of the number of dexels. Table 3 shows the RTF and number of dexels N D with different parameters measured after ten seconds of simulation time.
The parameters which have an influence on the RTF are the spindle speed Ω , the simulation time step size Δt To demonstrate the method's boundaries in terms of computational capacity on a certain hardware setup, Table 3 also covers parameter settings resulting in RTFs greater than 1.

Conclusion and outlook
In this article a new approach for simulating planar cutting processes and predicting chatter vibrations in real-time was introduced. The cutting forces were simulated by using a  dexel-based workpiece model for the calculation of the chip thickness. The number of dexels is adapted at runtime and converges to an sufficiently accurate spatial discretization of the workpiece. Thereby it is ensured that the chip thickness can also be calculated for varying cutting conditions and for arbitrary simulation time step sizes. Several scenarios for static and dynamic force simulation were validated by comparing the results with analytical solutions. Beyond that, a quantitative analysis of the method's computational efficiency via the real-time factor (RTF) confirmed its ability to be executed under hard real-time requirements. Fig. 16 Dexel-based cutting force validation near the stability limit with the critical depth of cut a lim = 16. 6[mm] at the exemplary spindle speed Ω = 600 rev∕min (see Fig. 15)

Fig. 17
Real-time factor ( RTF ) of an exemplary milling simulation One of the main advantages of using geometry-based models simulated in the time domain is that more complex applications with changing cutting conditions and arbitrary workpiece geometries can be considered. Currently the proposed model is being adapted on an out-of-round turning process whose first results are shown in Fig. 18. Out-ofround turning is characterized by permanently changing cutting conditions and unsteady process forces that can only be modelled using numerical approaches such as the here proposed method.
Another improvement will be the extension of the model for 3D processes coupled on multi degree of freedom (MDOF) machine dynamic models. Fig. 18 Simulation of an outof-round turning application as an example of a process with permanently changing cutting conditions and unsteady forces