Cutting force and vibration prediction of milling processes regarding the nonlinear behavior of cascade controlled feed drives

Milling processes are characterized by their high degree of flexibility regarding shape creation and thus occur in many value adding processes in the electrical and mechanical engineering industry. To facilitate the systematic prediction and the efficient optimization of the process performance, a holistic simulation model of the milling process with cascade controlled feed drives is developed. The integrated model consists of two subsystems, namely the feed drive model and the process model. In the feed drive model, the dynamic models of the electric motors, the serially connected torsional oscillators of the mechanical elements and the cascade control system are established. In this model, the coordinate point of the milling tool is computed at each discrete time step and assigned to the process model. The cutting forces and the corresponding load torques are calculated in the process model and sent back to the feed drive model. In addition, practical validation of the complete system has been conducted both in the time and the frequency domain. The simulation results yield an accurate prediction of the tooth passing frequency with a variance of less than 1%. Moreover, the influence of the cascade control parameters and the nonlinear behavior of the electric motor on the process vibrations have been experimentally verified.


Introduction
The globalization of markets and production is leading to increasing complexity of products and a growing number of product variants. Moreover, high flexibility and reliability in production are mandatory to ensure shortening product life cycles and more demanding quality standards. In this regard, the application of advanced manufacturing technologies as well as their digitalization and automation have become indispensable. The utilization of simulation technology, in particular, continues to gain in importance throughout the entire product life cycle. The modeling of complex systems, e.g. machine tools and machining processes, facilitates efficiency both in time and costs in manufacturing and product development [1].
With its high shaping flexibility, the milling process is of high significance in finishing processes. One of the most challenging aspects in milling are vibrations during the process, which both directly and indirectly increase the wear on the tool and machine components. The resulting degradation in workpiece surface quality and dimensional accuracy is also non-negligible [2].
In order to analyse the dynamic behavior of the milling processes with various parameters, a simulation model predicting cutting forces and vibration of milling processes has been developed and experimentally validated. The milling process model, integrated with the cascade controlled feed drives of the CNC milling machine axes, facilitates the correlative analysis of the nonlinear dynamics of the AC electric motors and the mechanical vibrations of the milling processes.

3 2 State of the art
In Sect. 2, firstly, the present state of research on simulation methods of the motion control system of the CNCmachine, including the modeling of electric motors and feed axes, is briefly described based on the literature. Secondly, drawing on the literature, different approaches of simulating the cutting forces and the vibration of machining processes are characterized and compared.

Simulation models of motion control systems
The physical behavior of the electric motors can be simulated applying differential equations. The solution of these equations can be attained in the time domain or as transfer functions in the Laplace domain. Various models are based on the coordinate system transformation from a threephase system into a rotational system, which simplifies the voltage and current equations with three components to two components [3]. Numerous approaches simulating different electric motor types and their practical implementation have been established in [4]. The modeling methods are carried out, for example, in Matlab Ⓡ and its integrated simulation environment SiMulink Ⓡ [5,6]. As the primary focus of this paper is the development of an integrated simulation system of the milling process with the control behavior, the reader is referred to the literature [7] for a more extensive survey of modeling electric motors.
The feed drive mechanisms can be simplified as a serial connection of several individual masses coupled by springdamper elements. This can be regarded as a multi-mass oscillator [8]. By modeling such oscillators, the dynamic behavior of the feed drive system can be directly derived from the differential equation which describes the motion of the oscillators in the time domain [3]. By establishing the general mechanical transfer function between the loads and motor position, the structural dynamic model of the drive system can be constructed. Dadalau et al. derived in [9] parametric equations to demonstrate the interrelationship between stiffness and geometrical parameters of the screw. In order to verify the influence of lateral dynamics of the screw on the drive's positioning precision, Okwudire et al. developed a hybrid finite element model of ball screw drives [10]. In [11], different compliance models of mechanical drives based on ball-screw and linear motors were presented.
There are further approaches on the software market, which facilitate the simulation of a machine tool and its subsystems. Specifically, the Virtual nC kernel (VNCK) from SieMenS can be used to simulate the operating sequence of a defined NC-program. During the course of the ReffiZ project, VnCk was applied as a platform and supplemented with program extensions which include the simulation of the mechanics and the process as well as an interface to a virtual drive application [12]. The company Pimpel GmbH developed the program CHECKitB4, which is applied to assess the machinability of the pre-designed workpieces as well as the selection of the machine and the tool. As a result, time losses and inaccurate process calculations can be prevented in early development stages [13]. A significant extension of the application scope can be realised through CHECKitB4 GIANT LEAP to simulate the automation and control technology of a virtual machine. This facilitates the testing and evaluation of a pre-designed NC-program to detect possible collisions and to verify the process runtime, the workpiece geometry as well as the motion behavior of the operational system in advance [14].

Simulation of machining processes
In order to model the machining process holistically, various simulation models and levels have to be coupled via a systematic approach. In this section, a brief overview of different methods of the material removal simulation of cutting processes is presented. Subsequently, diverse process modeling techniques for cutting forces and process vibration prediction are highlighted.
The material removal simulation is applied to determine the tool engagement on the workpiece material, including the chip formation as well as the varying depth of cut a p and the cutting width a e . Physical and mechanical characteristics, such as the cutting forces, the torques, the vibration behaviors as well as the workpiece shape changes, can be derived accordingly [15]. To model three-dimensional objects, solid modeling techniques are commonly applied. Altintas et al. described various solid-model-based systems, which is beyond the scope of this paper and the reader is directed to [16] for further details. In this review paper, moreover, discrete modeling techniques, namely the Voxel model [17], the Dexel model [18] and the Z-Buffer-based system [19] are clarified. Damm verified in [20] that a Dexel model is characterized by a reduced memory demand compared with a Voxel model and an improved accuracy as to the z-buffer model. Recent researches showed that by extending the Dexel model with the finite element simulation, it is possible to assess the thermo-mechanical shape and dimensional change of the workpiece [21,22]. Additionally, analytical methods are also adopted in tool engagement simulations. Tunc and Budak presented in [23] an analytical approach to estimate the workpiece surface information with the assistance of the cutter location file. Despite the short computation time of analytical methods, their application is limited to certain types of machining [16]. Due to the simple implementation of discrete models, they are widely used for tool and workpiece simulation [24]. On the contrary, solid modeling techniques are more suitable for applications with high accuracy requirements [25].
In order to compute the forces and vibrations of the machining processes, analytical methods are applied. In the simulation approaches of [26,27], the depth of cut is determined considering both the static and dynamic parts which are resulted from the tool geometry and the process kinematics as well as the dynamic vibrations. Kaymakci et al. [28] developed a unified cutting force model by transforming the forces into cutting edge coordinates for turning, drilling and milling processes. To simulate the process states at discrete time intervals, the time marching methods are first presented in [29]. Lazoglu et al. demonstrated in [30] the time-discrete calculation of cutting forces for five-axis milling operations. Surmann and Biermann developed in [31] a time-domain simulation system with reduced computation time, in which a model of surface formation of each tooth feed is derived using the CSG (constructive solid geometry) modeling.
Despite extensive approaches and applications on simulating the control system of machine tools and machining processes, knowledge about the interaction of the motion control parameters and the machining process dynamics remains limited. An individually developed simulation model for the representation of the motion control system of the machine tool and the cutting process yields several advantages over external software solutions. Firstly, it provides more flexibility in implementing and modifying functionalities when necessary. Secondly, the comprehensive access to the assigned variables and the transparency in evaluating applied functions are indispensable for the basic scientific research. Last but not least, the developed model achieves the analysis of interdependencies among control and process parameters in an isolated system without being disturbed by irrelevant influences, such as environmental surroundings and the operating temperature of the mechanical parts.

Cascade controlled milling process simulation
The focus of this work lies on an individually established two dimensional milling process model with the integrated motion control system of the feed drives. Based on the results of the preliminary experiments, the influence of the motion control system on the dynamic behavior of the milling process can be substantially represented by a simplified process model simulated on a two dimensional plane of the physical x-and y-axes of the machine. In addition, as the z-component of the cutting forces in planar milling, namely the passive force, is considerably smaller than the x-and y-components, it is sufficient for the purpose of this work to consider the cutting process as being two-dimensional. In the following section, the components of the cascade control system of the feed drives, the milling process calculation and their interfaces are described in detail.

Cascade control system modeling of milling machine feed drives
To establish the motion control system of the milling machine feed drives, firstly, the coordinate system transformation is illustrated. Secondly, the equations for modeling the electric motors are specified and their implementation is presented. The mechanical elements for the realization of the feed axis functions and the cascade control system are characterized subsequently. The individual components are then combined into a holistic system which is applied to update the position of the tool center point (TCP) for the milling process calculation at each discrete time step in Sect. 3.2.
First of all, to construct the dynamic models of the electric motors, the current and voltage values of the U-V-W three phase system are transformed into a rectangular coordinate system (COS) with two axes. In order to represent the stationary processes of an asynchronous motor, the stator fixed -reference frame can be applied by performing the Clarke transform (Eq. 1) [32]: For the demonstration of a synchronous motor, the phase values are to be transformed into a rotating reference frame using the Park transformation matrix [33]. Here, the coordinate axes d, q are aligned with the rotor flux of the machine. In Eq. 2, the transform of the three phase voltage values u U , u V , u W into the equivalent voltage values in a rotating reference frame is given. represents the rotating angle of the vector to follow the frame d, q attached to the rotor flux. Equation 3 gives the corresponding inverse transform: The Park transform applied in this work is simplified with the help of the Clarke transform (see Eq. 4). The matrix from Eq. 5 returns the inverse transformed components of (1) the stator fixed COS. Accordingly, the Park transform can be regarded as a series connection of the Clarke and rotational transforms [3]: This subsystem is implemented in a generic form so that it can be utilized to transform both voltage and current values without any additional adaption. Secondly, the permanent-magnet synchronous motor (PMSM) as the machine feed drive is described mathematically and configured. This can be divided into four subsystems respectively for the calculation of: • the rotor flux linkage Ψ of the d-and q-axis, • the current value i of the d-and q-axis, • the mechanical variables, namely the torque and the angular velocity and • the electrical rotor parameters for the Park transform.
Moreover, the implementation of the differential equations is presented. The voltage expressions define the d-and q-axis components of the stator voltage. The first summand indicates the voltage drop over the winding resistor R s , which is proportional to the stator current component i d,q . The second summand represents the induced voltage due to the changing magnetic flux Ψ d,q . The last term reflects the rotor angular velocity e dependent stator voltage, which is influenced by the flux linkage of the other axis. The equations are used to calculate the flux linkages between the rotor and the stator of the synchronous machine, where L d,q and Ψ PM represent the stator inductance and the permanent magnet flux linkage respectively [34]. Since the systems of Eqs. 6-9 are mutually dependent, feedback loops are implemented in the simulation model. The calculated Ψ d,q and i d,q are applied to determine the motor torque with the equation (9) and Ψ q = L q ⋅ i q where Z p is the number of poles [35]. The resulting mechanical angular velocity m of the motor shaft can be derived from the equation with J as the moment of inertia and M L as the torque at constant speed which is derived from the process forces from the milling process model. Subsequently, the rotor angular velocity e and the corresponding rotor position for the Park transform (Eqs. Last but not least, the control of the synchronous motors (see Fig. 1) and the coupled mechanical elements is constructed as a cascade control loop. The position controller, as the primary controller of the feed axis, drives the set point n set of the secondary controller, namely the speed controller. The actual position value x act and speed value act is derived from the aforementioned mechanics model. The tertiary controller is used to regulate the motor current and thus the torque, which are required to execute the PMSM model [36].

Milling process model with cutting force and vibration prediction
In the simulation of the milling process, the flow chart of the milling process model is shown in Fig. 2. The file in G-Code format is firstly interpreted into a matrix with discrete time steps, the coordinate points of the set toolpath, the spindle speed, the feed rate and a logical value of the cutting state. The time steps are variable and are adjusted to divide the designed tool path into an integer number of time steps. Both the rotational speed and the feed rate are taken into account. The workpiece contour is defined as a polygon. The defined tool path coordinate points are forwarded to the feed drive model as the set points for the position controllers of the feed axes. Subsequently, the calculated actual position values are fed back to the milling process model for the computation of the actual TCP, the cutting edge coordinates as well as the angles between the cutting edges and the x-axis. During the chip thickness specification, the intersection point between the line drawn from the TCP to the cutting edge and the current workpiece geometry is firstly computed at each time step. The chip thickness is derived from the distance between the current cutting edge and the intersection point of the previous step. The specific cutting force k c and the cutting force F c are calculated using the equations Here k c1.1 is k c with a chip width and thickness of 1 mm. h is the chip thickness calculated in the chip thickness specification. b is the chip width. m c is the gradient of the specific cutting force in a given range of chip thickness. K , C V , K st , K ver , K ss and K kss are the correction factor for the rake angle, the cutting speed, the chip upsetting, the wear occurring during machining, the cutting material and the cooling lubricant respectively [37,38]. Subsequently, considering the cutting edge angles, the cutting force in x-and y-direction F c x,y can be computed. The torque required to drive the load in x-and y-direction are calculated as with as the ball screw efficiency, P as the screw lead in m and the total axial force where is the coefficient of friction of linear guide, m is the mass being moved in kg and g is the gravitational acceleration. Thereupon, the torque at constant speed needed in the feed drive model can be calculated as the sum of the torque to drive the load M D , the torque due to preload M P and the torque due to friction of support bearings and seals M F (see Eq. 17): (14) and   Fig. 3. Here the implementation of a climb milling process applying a tool with four cutting edges is illustrated as an example. Each cutting edge is modeled as a polygon of four points, namely the TCPs and the corresponding cutting edge coordinate points of the previous and the current step. The cutting edge polygons are then subtracted from the workpiece polygon to generate the updated workpiece geometry. Including the TCPs at both time steps in the cutting edge polygons ensures that no minute residual polygon remains after the cutting operation with a high feed rate. Subsequently, the process model with cutting force computation has been calibrated by conducting a series of experiments with different workpiece dimensions and locations on the work table respectively. The complete simulation system can be exploited to analyze the interactions between the motion control system of the machine with the control parameters and the machining process dynamics, such as process forces and vibrations. The experimental validation for the integrated model, which is within the scope of this work, is presented in Sect. 4.

Experimental validation and discussions
The experiments have been conducted in the DNM 500 vertical machining center from Doosan Machine Tools Co., Ltd., which is located at the Institute of Resource and Energy Efficient Production Systems at the Friedrich-Alexander-Universität Erlangen-Nürnberg in Germany. The kinematics model of the machine is shown in Fig. 4a. To validate the simulated process vibrations, the tri-axial acceleration sensor W356B11/NC and the Apollo light measuring system are implemented to record the vibration amplitude during (17) the milling process. The schematic representation of the workpiece clamping in the milling machine is illustrated in Fig. 4b. The machine vice is rotated by 45 • clockwise so that the coordinated motion control of the two feed axes can be analyzed with a specified straight milling tool path. The workpiece of 60 mm in both length and width is made of the stainless steel X2CrNiMo17-12-2.
In order to investigate the influence of the control parameters on the 2-D process dynamics, a process with a relatively high feed per tooth has been designed and implemented. To allow for a sufficient 2-D model representation, tools with vertically overlapping cutting edges are not applicable. Therefore, the corner milling cutter M4132-040-B16-05-09 mounted with five inserts SDMT09T320-F57 WSP45S from the company Walter AG has been chosen for the experimental tests. The dimensions are shown in Fig. 5a. Figure 5b indicates the sequence of the validated end face milling process with the cutting tool with five edges, the tool path and the workpiece geometry during the process. The corresponding process parameters of the three milling path sections A, B and C are listed in the Table 1. The fixed-step size of the integrated simulation model is configured as 0.1 ms. The values of the parameters used in Eq. 13 are listed in Table 2 [39,40].
When the corner radius r of the cutting insert is unequal zero, the area of the uncut chip cross section A (see Fig. 6a) can be calculated as which is identical as the area of the chip cross section in Fig. 6b, where the corner radius of the insert equals zero. As a result, the corner radius of the insert does not affect the computation of cutting forces using the Kienzle formula (Eq. 14).
Firstly the comparison of the vibration amplitude of the simulation and the practical experiments in time domain is shown in Fig. 7. Due to the large amount of data, only 1.0 s  Fig. 7). By demodulating the vibration signals with envelope analysis, the significant process-relevant frequencies can be extracted. For an improved data interpretation, the frequency domain analysis of the vibration envelopes has been carried out implementing the Fast-Fourier-Transform (FFT). The corresponding vibration amplitudes of the path section A with a duration of 5.54 s are compared in the frequency domain in Fig. 8. Due to the fact that at higher frequency ranges vibration peaks do not occur, the diagram is limited to 500 Hz. The blue curve in Fig. 8    can also be observed in the simulation result (see Fig. 8 orange curve). Based on theoretical calculations, the eigenfrequencies of the simulated mechanical elements, namely the ball screw, the serially connected mechanical elements as well as the feed axis X and Y with the corresponding masses, are 425.9 Hz, 1154.2 Hz, 279.9 Hz and 199.1 Hz respectively. Due to the limited data quantity of the simulation and the dominance of the process vibrations, the eigenfrequencies of the mechanic elements are difficult to identify. However, the motor resonance frequency at 19.4 Hz is to be observed both in the experimental and the simulation results.
To further verify the impact of the feed drive cascade control on the process vibration behavior, the position controller gain K v is selected and modified. The simulation results of the frequency analysis of the vibration envelopes with different K v -factors are demonstrated in Fig. 9. It can be observed that, with a higher K v -factor, characterizing higher dynamic of the position controller, the vibration amplitudes at the motor resonance frequency 19.4 Hz decrease. On the other hand, the vibration amplitudes at the tooth passing frequency and its harmonic 66.6 Hz and 133.3 Hz increase with higher K v -values.

Summary and outlook
This paper presents a new holistic simulation model of the milling process with cascade controlled feed drives. In the process model, the calculated cutting forces are applied to compute the load torques, which are necessary for the feed drive modeling. Considering both the dynamics of the mechanical elements, simulated as a serial connection of torsional oscillators, and the nonlinear behavior of the electric motors, the position of the TCP is fed back to the process model at each discrete time step. In addition, the workpiece geometry is updated at each step by subtracting the polygon shape of the tool cutting edges from the workpiece polygon.
To validate the simulation model, practical experiments have been carried out on a milling machine. The results indicate that the simulation model yields an accurate computation of the tooth passing frequency with only less than 1% deviation. To improve the model precision, further component extensions should be developed in simulating the main spindle dynamics and in calculating the cutting forces which contain the Fourier frequency components.
Furthermore, by varying the position controller gain, the distribution of the vibration amplitude in the frequency domain changes. Thus, the influence of the cascade control system and the nonlinear electric motor behavior on the process dynamics can be verified. Further research focus will lie on the systematic investigation of the interrelationship of the control parameters of the electric motor and the machining process dynamics. The feasibility of reducing process vibration through the improved configuration of the control parameters should be more comprehensively examined.
Funding Open Access funding enabled and organized by ProjektDEAL. 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/.