Chatter prediction for uncertain parameters

The occurrence of chatter in milling processes was investigated in this study. The prediction of the stability lobes of metal cutting processes requires a model of the cutting force and a model of the dynamic machine tool behavior. Parameter uncertainties in the models may lead to significant differences between the predicted and measured stability behavior. One approach towards robust stability consists of running a large number of simulations with a random sample of uncertain parameters and determining the confidence levels for the chatter vibrations, which is a time-consuming task. In this paper, an efficient implementation of the multi frequency solution and the construction of an approximate solution is presented. The approximate solution requires the explicit calculation of the multi frequency solution only at a few parameter points, and the approximation error can be kept small. This study found that the calculation of the robust stability lobe diagram, which is based on the approximate solution, is significantly more efficient than an explicit calculation at all random parameter points. The numerically determined robust stability diagrams were in good agreement with the experimentally determined stability lobes.


Introduction
Chatter vibration is one of the strongest limitations in maximizing the material removal rate of the metal cutting processes in industrial applications. Large amplitude chatter vibrations are undesired because they lead to bad surface finish, increased tool wear, and noise [1]. Particularly in milling processes, the stability lobe diagram (SLD) is a very powerful tool for the optimal selection of process parameters and cutting conditions. For example, the stability lobes specify the boundary between the chatter-free (stable) cutting and (unstable) chattering motion in the parameter space spanned by the spindle speed and depth of cut [1]. The SLD computation requires the modeling and parameterization of the process force and the dynamic machine tool behavior [2]. Moreover, the model parameters can be determined by measurements. However, any measurement procedure is subject to uncertainties [3,4]. Additionally, the dynamic behavior of the real machine tool during the process is influenced by the current operating state and environmental conditions. Consequently, a repeated measurement to determine the SLD can provide different results and a repetition of the same machining operation may result in different stability behavior. Thus, for a robust process design, it is important to consider the parameter uncertainties during the calculation of the SLD.
Methods of robust stability prediction with respect to parameter uncertainties have already been reported by several authors. Kong et al. [3], Kong and Yu [4], Sims et al. [5] and Hamann et al. [6] have used fuzzy set theory to determine the confidence ranges of the computed stability diagrams. Sims et al. [5] reported that explicit fuzzy interval-based arithmetic is impractical in the analysis of machine tool vibrations. Instead, numerical methods with a specific sample of uncertain parameters have been used [5,6]. Park and Qin [7] have used the edge theorem combined with the zero exclusion method in robust stability analysis. This guarantees that stability transition does not exist in a polytope within a parameter space spanned by the boundary values of uncertain parameters. Totis [8] used a probabilistic approach to describe parameter uncertainties. Hajdu et al. [9] investigated the robust stability behavior of turning processes by considering the uncertainties directly in the frequency response function (FRF) of the machine tool structure, instead of considering the model parameter uncertainties. These studies revealed that a very significant numerical effort is required in robust stability analysis, and that the calculation time increases as the number of the uncertain parameters increases. A simplified model has been frequently used in the prediction of the SLD. Moreover, it is necessary to use very efficient numerical methods for the stability analysis (e.g., as in Refs. [5,8]) and to select only a few representative sample points from the set of uncertain parameters [5,7].
There are different methods of calculating the SLD, and they can be categorized into time domain and frequency domain methods [10]. Efficient time domain methods include the spectral element method [11] and the Chebychev collocation method [12]. On the other hand, the zeroth order approximation (ZOA) in the frequency domain is very fast [13]. However, in this approach, only the average time periodic coefficient matrix is considered, and this may result in significant deviations in the SLD calculated from the exact SLD [14]. An improvement is possible by considering multiple frequency components [15]. This socalled multi frequency solution (MFS) leads to results that are as accurate as those obtained by comparable time domain methods [14]. The main advantage of the frequency domain methods is that a large number of structural eigenmodes can be considered without a significant increase in the calculation time, because the FRFs are used as a model of structural dynamics. Consequently, for industrial applications that are often characterized by multiple dominant eigenmodes, the MFS is the most efficient and accurate approach towards computing the SLD.
In this study, we used an efficient implementation of the MFS to determine the confidence levels in the SLD. The effect of the uncertain parameters was investigated using the approach proposed for the ZOA [16]. In particular, in this approach, the dynamic behavior must be determined only for a small number of sample points in parameter space, and an approximate solution is used to characterize the dynamic behavior in the remaining parameter space. The approximate solution is then used for the efficient calculation of the robust stability behavior. This paper is organized as follows. The basic concept of the chatter prediction for uncertain parameters is introduced in Sect. 2. The numerical method for the SLD calculation of a single point in parameter space is presented in Sect. 3. The details of the approximate solution and its benefit in comparison to the exact solution are presented in Sect. 4. The application of the proposed method to the cutting tests is described in Sect. 5 and the conclusions are presented in Sect. 6.

Basic concept of robust stability analysis
The proposed approach is based on the analysis of the dynamic behavior of the cutting process by an efficient implementation of the MFS, which is presented in detail in Sect. 3. The exact solution through the multi frequency approach is only calculated for an appropriate set of data points in the given parameter space, and the results are used to generate an approximate solution of the MFS. An optimal choice of data points and details regarding the interpolation scheme are described in Sect. 4.
Specifically, to predict the stability of metal cutting processes, the interactions within the closed loop of the coupled subsystems machine tool and process are relevant (see Fig. 1). The shape of the stability lobes, i.e., the limiting axial depth of cut a p versus the spindle speed, depends on the parameters of the dynamic model of the machine tool and on the process parameters. The typical parameters of the machine tool model are the mass, damping, and stiffness values of the dominant modes, while the process parameters include the cutting force coefficients, entry and exit angles of the cutting teeth, and geometric properties of the cutting tool. The first step is to select the relevant uncertain parameters and to quantitatively characterize their distribution. It is useful to consider only the most important uncertainties because each additional uncertain parameter increases the dimension of the parameter space by one, and the computing time for the robust stability analysis increases with the parameter space dimension.
Various analysis methods have been developed for the SLD computation at a single data point in the parameter space. The analysis in the frequency domain can be carried out by considering the transfer function G 0 of the open loop, which was reported in Refs. [2,13] where x is the angular frequency; s is the time delay equivalent to the tooth passing period; matrix A contains the directional factors, and G is the transfer function of the machine tool model describing the dynamic displacements at the tool center point (TCP) in response to a force at the TCP. The term ðe Àixs À 1Þ represents the so-called regenerative effect, which appears because the dynamic chip thickness is equal to the difference between the dynamic displacements at the previous cut at time t À s, and the displacements at the present cut at time t. The transfer function of the process is given by the product a p ðe Àixs À 1ÞA, where the cutting force coefficients and the transformations between the machine tool coordinates and process coordinates are summarized in matrix A. For the stability analysis of the process the Nyquist criterion can be applied to the open loop transfer function G 0 . For this purpose, the system coordinates must be decoupled by the eigenvalue decomposition of matrix G 0 , and the stability lobes can be obtained from the individual eigenvalues of matrix G 0 . The calculation and eigenvalue decomposition of the open loop transfer function G 0 is responsible for the main part of the numerical effort. To save computing time, the basic idea is to calculate G 0 and its eigenvalues only at selected data points and use them for an approximate solution (see Fig. 2). The open loop transfer function G 0 is calculated at selected points of the parameter set and stored (left). A random set of points in parameter space is chosen and the stored data are used to approximate the stability lobes for each parameter point (top right). The confidence levels can be calculated from the probability of stable and unstable cuts at each spindle speed and depth of cut (bottom right). In a direct approach without approximation, the number of data points for a representative sample of the parameter space increases exponentially with its dimension [8]. In Sect. 4, we discuss how an accurate approximate solution can be constructed with significantly less sampling points by considering that the effect of many parameters on the dynamic behavior of the process is independent of each parameter. In other words, the change of the dynamic behavior for a variation of multiple parameters can be derived from the behavior wherein each parameter is varied separately. The approximate solution can then be used to determine the SLD in a very efficient way for a much larger number of parameter points without losing much accuracy. Therefore, the method is suitable to the calculation of a robust SLD with the confidence levels of the stability behavior.

Computation of open loop transfer function
For the stability analysis of the metal cutting process at a single data point, the dynamic cutting force and dynamic displacements at the TCP are relevant. The dynamic components are considered to be small perturbations of the static or quasi-static displacements and forces; therefore, a linearized model was used in the stability analysis [17]. The dynamic cutting force in the machine tool coordinates can be given by the following equation: where x(ts) and x(t) specify the dynamic displacements of the tool at the previous and present cut, respectively. The Chatter prediction for uncertain parameters 321 matrix B(t) contains the cutting force coefficients and the directional factors related to the transformation between the machine tool and process coordinates (for details see Refs. [17,19]). For the milling processes, these coefficients, and therefore the matrix B(t), are time periodic with a tooth passing period s. Note that the force F(t) in Eq. (2) contains forces from all cutting teeth. Equation (2) is a time domain representation of the lower part of Fig. 1 including the transformations from the machine tool to the process coordinates and vice versa. Since the matrix B(t) is generally time periodic, i.e., B(t) = B(t ? s), a transformation to the frequency domain is not straightforward. In fact, a multi frequency approach is necessary. From the Floquet theory we know that the dynamic cutting force can be decomposed into a periodic part and a remaining exponential part, and expressed as follows [14,18] FðtÞ ¼ Á Á Á e ÀiXt 1 e iXt Á Á Á À Á Fe kt ; where X = 2p/s is the tooth passing frequency,F k the Fourier coefficients of the periodic part of the cutting force, and k the Floquet exponent specifying the exponential behavior of the perturbations. Similar statements hold for the dynamic displacements x(t). We assume a vanishing real part of the Floquet exponent, i.e., k = ix c , which holds, among others, for the processes operating directly at the stability boundaries. The Fourier transformation of Eq.
(2) contains a convolution, where only a discrete set of the frequency peaks appear. The discrete convolution can be constructed as a matrix multiplication and the frequency domain representation of Eq. (2) can be given by the following equation: The vectorx is defined similar toF, i.e., it contains the Fourier coefficients of the periodic part of the dynamic displacements x(t).B k represents the Fourier coefficients of Consequently, we define the multi frequency transfer function G(x) for the frequency domain representation of the machine tool's dynamic behavior [19].
whereĜ(x) denotes the classical FRF at the TCP of a single frequency component. By incorporating Eq. (5) into Eq. (4), we obtain the following frequency domain representation for the following closed loop of the process machine interactionŝ with the open loop transfer function defined in Eq. (1).

Stability lobes and modal assurance criterion
At the stability boundary, the closed loop expressed by Eq. (6) is fulfilled, which means that at least one eigenvalue of the open loop transfer function G 0 (x) is equal to 1. For a numerical implementation, the infinite dimensional matrices in Eq. (6) must be truncated. It is recommended that the dimension is chosen such that the entries in the multi frequency transfer functionĜ(x) cover all of the dominant eigenmodes of the structure. The numerical implementation details of the MFS have been provided in Refs. [18,19]. Here, we define the eigenvalues K(x) = K R (x) ? iK I (x) of the truncated matrix function e Àixs À 1 À Á AGðxÞ. First, the critical frequencies x crit with K I (x crit ) = 0 were determined. These were frequencies with a vanishing imaginary part of at least one of the eigenvalues. Because the metal cutting dynamics are characterized by a time delay system with a transcendental characteristic equation, there are infinitely many points in parameter space, where the real part of one Floquet exponent of the delay system vanishes, or equivalently, where K I (x crit ) = 0 [18]. At the critical frequencies, Eq. (6) is simplified as follows: which means that the stability boundary, i.e., the lowest critical depth of the cut where Eq. (6) holds, is determined by the eigenvalue with the largest real part, namely, K R (x crit ) [17]. This corresponds to the critical frequencies close to the eigenfrequency of one of the structure's dominant modes. Thus, the resolution of the frequency grid must be specified according to the width of the dominant peaks of the FRFs at the TCP. If only the average directional factorsB 0 are considered, Eq. (6) is equivalent to the ZOA. In this case an analytical expression for the eigenvalues K(x) can be obtained [17]. However, for the MFS, the number of the higher harmonics and the dimension of the matrices increase. This results in longer computation time for the MFS, in comparison to the ZOA. On the other hand, the results of the MFS are more accurate than the ZOA, which is especially significant for the cutting processes with strong parametric excitations, such as the low radial immersion milling. Another problem is that, in the MFS, the analytical expressions of the K(x) eigenvalues are generally not available. Instead, the eigenvalues must be determined numerically, which leads to a problem with the identification of the critical frequencies (x crit ) by finding the crossings of the K(x) eigenvalues with the real axis. In particular, from a numerical eigenvalue decomposition, it is typically unclear which K(x k ) eigenvalue at the x k frequency is associated with one K(x k?1 ) eigenvalue at the next x k?1 frequency. This makes it impossible to identify the eigenvalue crossings with the real axis. We solved this problem using the modal assurance criterion (MAC) to relate the eigenvalues at one frequency to the correct ones at the next frequency. In particular, the eigenvector matrix V(x) was calculated along with the K(x) eigenvalues, where each column of the V(x) matrix contained one normalized eigenvector. The MAC matrix of the relationship between the V(x k ) and V(x k?1 ) eigenvectors at the two frequencies of x k and x k?1 is given by where * denotes the conjugate transpose and denotes the elementwise product or the Hadamard product. By rounding the MAC values in the M matrix in the vicinity of one or zero, we obtain a permutation matrix, which can be used directly to accurately connect the numerically determined K(x k ) and K(x k?1 ) eigenvalues from the two frequency steps. From the ordered eigenvalues, the identification of the critical frequencies is possible in the same way as for ZOA, and the stability lobes can be calculated by Eq. (7). The MAC matrix M with the MAC values was also used to identify the correct relationship between the K eigenvalues at different parameter points.

Robust stability analysis with approximate solution
In the robust stability analysis, it was necessary to determine the stability of the cutting process for a representative set of points in parameter space. To save computing time, we propose an approximate solution based on the exact stability result calculated by Eq. (7) at several selected data points only. The approximation required a continuous function for the varying parameters in the range of the considered values. However, this was not fulfilled for the axial depth of cut a p . For example, if the stability at the neighboring data points was determined by the different K(x) eigenvalues, an approximation based purely on the limiting axial depth of cut a p at the neighboring points would result in inaccurate results. Therefore, we did not only use the values of a p for the approximation, but rather used the full information regarding the K(x) eigenvalues as the basis of the approximation.

Construction of approximate solution
In this study, we focused on the parameter uncertainties of the machine tool's dynamic behavior and did not consider the uncertainties of the force coefficients. The FRFs of the dynamic displacements at the TCP in theĜ(x) matrix can be expressed as the superposition of the individual harmonic oscillator FRFs. Three parameters are required to describe a single oscillator. In particular, theĜ ij element of the FRF matrix can be given bŷ with the modal compliance N m , damping ratio D m , and tuning ratio g(x) = x/x 0m between the frequency x and the eigenfrequency x 0m of the mth oscillator. Thus, the variation of the modal parameters affects the elements of the machine tool structure's transfer function and results in the change of the K eigenvalues. As mentioned above, the K(x) eigenvalues were calculated numerically, which meant that the effect of the parameter changes on the stability behavior could be calculated analytically. We defined the vector p and specified a point in parameter space. More precisely, the elements of p were a specific realization for all of the considered uncertain parameters p i . In other words, p contained the modal compliances, damping ratios, and eigenfrequencies of all dynamic machine tool model oscillators. The approximate solution was based on the interpolation and extrapolation of the K(x) eigenvalues of the open loop transfer function between and beyond a few selected sampling points. Here, one sampling point was specified by the parameter vector p 0 , for which all individual parameters p i adopted their expectation values p i;0 . The eigenvalues were calculated for the parameter point p 0 . A variation Dp in the parameter space resulted in the change of the eigenvalues from K(x; p 0 ) to where R(Dp) specifies the relative changes of the eigenvalues depending on the parameter variation Dp. The correct relationship between the eigenvalues K(x; p 0 + Dp) and K(x; p 0 ) was obtained from the MAC matrix, as described in Sect. 3.2. The relative changes of the eigenvalues depending on the parameter variation Dp were approximated by a second order polynomial, as follows To calculate the polynomial coefficients, the K(x; p) eigenvalues were calculated at the two additional grid points p j À and p j þ for each dimension j of the parameter space. The additional grid points were selected according to the standard deviations r i of the ith parameter and could be specified by p j À;i ¼ p j þ;i ¼ p 0;i for i 6 ¼ j and p j AE;i ¼ p 0;i ð1 AE r i Þ for i = j. A piecewise linear approach was also implemented as an alternative to the proposed second order approach. This provided slight advantages in terms of calculation time increased the approximation error. In principle, an extension to higher order or spline polynomials was possible.
The approach expressed by Eq. (11) did not consider any interaction between the individual parameters. For the examples considered here, the interaction between the parameters of the different modes could be largely neglected. This was based on the fact that the FRFs at the TCP in Eq. (9) were a superposition of the individual oscillator FRFs, even though the resulting K(x; p) eigenvalues may have depended nonlinearly on theĜ ij elements of the FRF matrix. Nevertheless, the interactions between the parameters of a single mode were relevant. This concerned the mixed variation of the natural frequency and the damping ratio of one oscillator in particular, while the effect of the corresponding modal compliance was independent of the changes in the other parameters. This was attributed to the fact that the compliance appeared only linearly as a proportional factor in Eq. (9). This was advantageous in defining the oscillators by means of mass, stiffness, and damping values, where more mixing terms appeared. To consider the effects of the mixing terms on the approximate solution, Eq. (11) was extended where the index m in the second product adopts only values corresponding to the damping ratios and the index m ? 1 corresponds to the eigenfrequency of the same oscillator. Thus, R m;mþ1 (Dp) specified the relative changes of the eigenvalues, owing to the parameter interactions of the mth oscillator. These are given as To identify the four coefficients in Eq. (13), the eigenvalues had to be calculated explicitly at four additional parameter points with the method described in Sect. 4. For these four additional grid points, we used the four combinations that could be obtained by p 0 þ p m AE þ p mþ1 AE . To determine all of the approximate solution coefficients, explicit calculations were necessary at the operating point p 0 , the two grid points for each of the three individual parameters of one mode, and the interaction terms of four grid points per mode. Thus, the number Z of the explicit calculations that were required to determine the coefficients of the approximation solution is given as where M is the total number of modes. Note that the dimension of the parameter space is D = 3M. Because that many mixing terms are negligible, the number of the required calculations to determine an approximate solution in the 3M dimensional parameter space scaled only linearly with the dimension. In particular, for a large number of uncertain parameters, this was an extreme reduction in comparison to the 3 D number of the calculations that were required by the method proposed in Ref. [8], where three grid points were used for each parameter and all parameter interactions were considered. From Eq. (14) and the fact that the eigenvalues are typically complex, it follows that 20M parameters must be stored for a single eigenvalue at one spindle speed. This value must be multiplied by the number of eigenvalues and the dimension of the frequency grid because it was necessary to store the behavior of each K(x; p 0 ) eigenvalue of the open loop transfer function for each frequency point x k . A reduction in the required memory space would be possible if the approximate solution was not constructed for all eigenvalues, but rather for only a small part of the eigenvalues. The selection of the relevant eigenvalues could be carried out, for example, according to the real part of the eigenvalue at the operating point p 0 because only the eigenvalues with the smallest real part could determine the stability behavior, as discussed in Sect. 3.2.

SLD for uncertain parameters
To carry out a robust stability analysis it was necessary to analyze the stability behavior with respect to the parameter changes. For the calculation of the confidence levels in the stability chart, it was necessary to calculate the probability of a stable process for a fixed spindle speed and depth of cut. If a certain probability distribution was assigned to each uncertain parameter, this may have been possible through an analytical approach because the approximate solution was specified by a multivariate polynomial. The problem with the analytical integration was that the integration boundaries over the parameter space were not explicitly known. In particular, they were given implicitly, depending on the real part of the dominant eigenvalues. Therefore, we used a Monte Carlo simulation for the calculation of an SLD with confidence levels. An adequate number of sampling points in parameter space was generated from the probability distribution of uncertain parameters. The real and imaginary parts of the eigenvalues were calculated from the approximate solution of each parameter point at each spindle speed. The limiting axial depth of cut a p was determined for all parameter realizations with the method described in Sect. 3.2. A robust SLD could be obtained after sorting the resulting value of a p for each spindle speed.

Computational time and approximation error
The linear approximate solution for the ZOA method has already been presented in Ref. [16], where it was applied to a process with a large number of uncertain parameters. The number of uncertain parameters was 162. In this study, the approximate solution was extended using the MFS for the explicit calculation of the stability boundary and a second order polynomial approach for the approximation. New numerical results were available for a reference process with a medium number of uncertain parameters. The number of uncertain parameters for this process was 24. The numerical effort of both methods and the two reference processes are discussed below. The reference processes with the medium and large number of uncertain parameters were labeled as M and L, respectively. Tables 1 and 2 summarize the results obtained by the ZOA method. The presented algorithms were implemented in Matlab. The explicit calculation of the SLD at a single parameter point using the ZOA method for the reference process L lasted 0.1 s (Intel Core i7 2 600 K, 3.4 GHz; 8 GB RAM). The advantage with regard to the computation time for the approximate solution was relatively low and decreased for a smaller number of uncertain parameters. Regarding the occurrence of an approximation error and the increased storage requirement, the approximate solution was not necessarily favorable in comparison to the explicit calculation.
For the MFS, the choice of the Fourier order was important with regard to the calculation time. The ZOA algorithm computed the eigenvalues of a two-dimensional matrix at each frequency. Moreover, the dimension increased for the MFS. The investigations regarding the influence of the Fourier order revealed that, for the considered cutting conditions, a Fourier order of h r = 6 represents a good compromise between the accuracy and the calculation time. More details on the selection of the Fourier order for the MFS can be found in Sect. 3.3 of Ref. [18] and Fig. 3 in Ref. [19]. With a Fourier order of h r = 6, the matrix of the transfer functions of the open loop had 26 dimensions. Thus, the approximate solutions for the 26 loci had to be determined. Tables 3 and 4 summarize the results of the numerical effort for the reference processes L and M, respectively. In comparison to the explicit calculation, the approximate solution reduced the calculation time to approximately 9% of the linear approach or to 14% of the polynomial approach. Generally, only some of the 26 loci exerted a potential effect on the stability limit. By considering only the six most important loci, the calculation time was further reduced with only slightly decreased accuracy. For the polynomial approach with only the six most important loci, the calculation time decreased to approximately 8% of the required computation time in the explicit solution. In contrast to the ZOA method, the advantage of the approximate solution was also significant in the reference process M with a medium number of uncertain parameters. This occurred because the calculation time for the MFS was essentially determined by the Fourier order and not by the number of uncertain parameters, owing to the fact that, for the MFS, the numerical eigenvalue decomposition was responsible for a large part of the computation. The computation time for the explicit calculation of an individual and randomly chosen realization of the parameter vector p was only slightly different for the two models; namely, 4.8 s and 4.3 s for the reference processes L and M, respectively. According to Eq. (14), the numerical effort of the approximate solution decreased approximately linearly with the number of the model parameters. The best result could be obtained for the reference process M, where the calculation time with the approximate solution was reduced by up to 0.7% of the calculation time for the explicit solution.
In comparison to the explicit calculation, the disadvantage of the approximate solution was the error of the stability result. Figure 3 shows the relative error of the approximate solution for the 1 000 parameter points of the randomly determined model parameters. For the reference process L, the calculations were performed with the MFS. As expected, the approximate solution that was implemented with a second order approach led to a lower error than that of the approximate solution with the linear approach. The maximum error of the second order polynomial approach was less than 4%.

Dynamic machine tool behavior
The measurements were carried out on a DMU 80 eVo linear (please kindly check it) (DMG / Mori Seiki). To compute the stability lobe diagrams, the dynamic behavior at the tool center point (TCP) was represented by an oscillator model. The modal parameters were identified by the measured FRFs. The maximum speed of the spindle was 18 000 r/min. It was found that the dynamic behavior of the spindle was speed dependent; therefore, the measurements were performed on the rotating spindle. An impact hammer carried out the excitation,  during which, the frictional forces could have occurred in the contact zone owing to the relative movement between the rotating target and the hammer tip. Thus, transverse forces, which can influence the measurement result, may have been introduced. The measurements made here used a steel hammer tip. In previous investigations, it was shown that the frictional transverse forces under these conditions could be neglected [20]. However, the measurement of the displacement had to be carried out without contact between the sensor and the tool. Thus, capacitive sensors were used. A further limitation was that the measurement of the rotating spindle with the milling tool could not take place at the TCP. Therefore, a cylindrical dummy tool was used. For the non-rotating spindle, it was possible to measure the FRFs for the actual tool. Figure 4 shows the comparison of the measured FRFs for the dummy and original milling tool at zero spindle speed (n = 0). The figure also shows the cross FRFs for the dummy tool, for which two dominant eigenfrequencies occurred at approximately 2 900 Hz and 3 700 Hz. There also existed two dominant modes for the actual tool. However, the frequencies changed to 3 200 Hz and 4 100 Hz, respectively, and the resonance amplitudes also changed significantly. Therefore, using the speed of the dummy's dependent modal parameters resulted in errors in the calculated stability limit. The measuring of the rotating spindle with a dummy tool provided the speed dependent modal parameters, that is, the natural frequency f Dummy,m (n), damping ratio D Dummy,m (n), and resonance compliance N Dummy,m (n) (see Fig. 5). The measurement of the actual tool provided the corresponding parameters, but only for the spindle speed of n = 0. Therefore, the speed dependent natural frequency of the tool, namely f Tool,m (n), became f Tool;m ðnÞ ¼ f Tool;m ðn ¼ 0Þ f Dummy;m ðnÞ f Dummy;m ðn ¼ 0Þ : The damping ratios D Tool,m (n) and resonance amplitudes N Tool,m (n) of the tool for the spindle speeds of n [ 0 were calculated in the same way. This simple approach towards determining the changing FRFs at the tool center point of the varying spindle speed was based on the following considerations. The dominant eigenfrequencies of the system with the original tool and dummy tool were close to each other (see Fig. 4). Therefore, it was assumed that the dominant mode shapes of the vibrating spindle-holder-tool system and the system with the dummy tool were comparable (see Fig. 6). Since the changing eigenfrequency of the rotating tools was mainly caused by the changing stiffness of the spindle bearings [21,22], the authors concluded that the ratio of the eigenfrequency of one of the rotating spindle's eigenmode to that of the static spindle did not significantly depend on the detailed geometry of the tool (actual tool or dummy tool).
In the stability analysis, only the two dominant modes were relevant. Therefore, the parameters of the two harmonic oscillators were chosen such that the FRFs of the oscillators approximated the measured FRFs around the resonant peaks. This procedure was carried out separately for the direct compliances G xx and G yy and for the cross compliances G xy and G yx because a separate fitting of the FRFs would result in more accurate results [23]. Consequently, the entire dynamic machine tool model consisted of eight oscillators with three parameters for each oscillator, which resulted in 24 uncertain parameters. This model corresponded to the reference process M that was mentioned in Sect. 4.3.
Because we are interested in the parameter uncertainties of the dynamic machine tool behavior, the FRFs were measured under different test conditions, and at least 10 measurements were carried out in each case. The investigated test conditions were as follows: (i) Acceleration sensor applied with wax and impulse hammer excitation by hand; (ii) Non-contact measuring, capacitive transducer, and impulse hammer excitation by hand; (iii) Non-contact measuring, capacitive transducer, impulse hammer excitation by means of a device for the better reproducibility of the force pulse; (iv) Non-contact measuring, capacitive transducer, automatic tool change between each individual measurement.
To determine the standard deviations of the modal parameters, only the direct FRFs in the x-direction were recorded. A total of 90 measured FRFs with various test conditions were measured and evaluated. The parameters of the oscillators were fitted for each FRF measurement. A normal distribution was assumed for the modal parameter uncertainties. The mean value and standard deviation were estimated for each individual parameter. The best reproducibility was obtained for the natural frequencies of the system, while the parameter distribution was much broader for the damping ratios and resonance amplitudes (see Fig. 7). The uncertainties represent measurement errors, errors of the parameter identification algorithm, and physical changes in the parameters of the actual structure.

Experimental set-up for cutting tests
The simulation and experimental determination of the SLDs was conducted for the up-milling of an aluminum AA7075 workpiece with a cylindrical end mill. The radial and tangential cutting force coefficients k r and k t were identified in preliminary cutting tests. The uncertainties of the cutting force coefficients were not investigated. The process parameters were: carbide endmill, tool diameter d = 12 mm, helix angle g = 30°, number of teeth Z t = 4, workpiece aluminum AA7075, radial immersion a e = 3 mm (25%), tangential force coefficient k t = 804 N/ mm 2 , radial force coefficient k r = 171 N/mm 2 .
To experimentally determine the SLD, the workpieces were clamped onto a three-channel Kistler Type 9257B dynamometer. As can be seen in Fig. 8, the pure peripheral milling was performed without the engagement of the secondary cutting edges. In other words, the workpiece width specified the axial depth of the cut. Moreover, the axial depths of the cut with a discretization of 0.5 mm and rotational speeds ranging from 11 500 r/min to 16 000 r/ min with a spacing of 250 r/min, were investigated by cutting tests. At each spindle speed n, the axial depth of the cut increased iteratively until the process behavior changed from stable to unstable. After each individual cutting test, an automatic tool change was performed. In particular, the used milling tool was placed in the tool magazine of the machine tool and was clamped again. All of the experiments were performed with the same tool. However, the wear of the tool was not investigated.
The process forces that were recorded during each machining test were used to identify the process stability. In the case of a stable process, the rotational frequency f n , tooth passing frequency f TP , and their harmonics, were found in the power spectra of the process forces. When an unstable process occurred, significant vibration components were observed at the chatter frequencies (f chatter ) and were identified by not being a multiple of the rotational frequency [24]. A process was declared as unstable when the ratio of the magnitude at the chatter frequency to the magnitude at the tooth passing frequency was Mag (f chatter )/ Mag(f TP ) [ 0.005, which turned out to be a meaningful choice in the cutting tests (see Fig. 9). Figure 10 shows the results of an experimental SLD consisting of the stability result for each individual cut. Strictly speaking, it would not be very meaningful to present a stability boundary between the stable and unstable cuts because each individual experiment was subject to uncertainty. In fact, the objective of the cutting tests was to determine the occurrence probability of an unstable process. This could be obtained if various cutting tests were performed at each rotational speed and an axial depth of cut. In this study, four test series were carried out, but they were not sufficient for the calculation of a representative probability. Instead, a presentation through four stability boundaries was chosen, as shown in Fig. 11.

Stability lobe diagrams
Each of the four cutting test series was carried out within one day. However, between the test series, several days or even several weeks passed and the experimental   setup was removed from the machine. In all experiments, the same tool and holder were used. To reproduce the overhang length as much as possible, the tool was always inserted into the tool holder with the smallest overhang length. In all cutting tests, dry machining was applied.
For the considered cutting process, the calculation by the ZOA method already provided sufficiently accurate results. However, if the radial immersion or the helix angle of the mill was reduced, for example, the accuracy of the ZOA decreased and the MFS became more advantageous. Figure 12 shows the confidence levels calculated with the ZOA method, and the four experimentally determined stability limits. The experimentally determined SLDs concerned most spindle speeds within the 95% confidence intervals of the numerical computation. For the spindle speeds of n & 15 500 r/min, the experimental stability boundary was outside the 95% confidence intervals. A similar result was obtained only if the modal parameters of the non-rotating spindle were used in the calculation of the confidence levels (see Fig. 13). A robust process design from the robust SLD of the simplified model of the nonrotating spindle model shown in Fig. 13 resulted in the lower axial depths of the cut and, therefore, to lower productivity. Moreover, the agreement of the experimental results with the confidence levels of the speed-dependent machine tool model shown in Fig. 12 was better than the agreement with the confidence levels from the non-rotating spindle model shown in Fig. 13. This means that the approach towards determining the speed-dependent dynamic machine tool model parameters, which is expressed by Eq. (15), was reasonable. Figure 14 shows the calculated confidence levels for the MFS. As has already been noted above, the example process was not an actual low immersion milling. Therefore, in comparison to the ZOA calculations, the differences in the calculated confidence intervals were expected to be small. However, it was shown that the approximate solution presented here could also be successfully applied to the MFS.

Conclusions
The proposed method enables the efficient calculation of confidence levels for the stability behavior of metal cutting processes. The method was based on the stability analysis in the frequency domain. The main idea was the Fig. 9 Cutting force spectra of stable and unstable cutting test; test series 2; n Spindle = 12 250 r/min Fig. 10 Experimental cutting test results for one test run with boundary between stable and unstable cuts Fig. 11 Experimental stability boundaries from four cutting test series (run 2 corresponds to Fig. 10) development of an approximate solution, which required explicit stability calculations at only a small number of parameter points. The number of calculations that was required to construct the approximation increased only linearly with the space dimension of uncertain parameters. This is a significant advantage in comparison to other robust stability analysis approaches. For the presented reference process, the maximum error between the approximation and the explicit calculation was kept below 4% using a second order polynomial approach for the   approximate solution. The method is particularly advantageous for a stability analysis with the MFS and for a large number of uncertain parameters. The confidence levels could be calculated very efficiently from the approximate solution using the Monte Carlo method. The knowledge of the confidence levels allowed the robust design of milling processes. The standard deviations of the modal parameters, which were determined in a series of tests, resulted in a good representation of the actual cutting process uncertainties. However, it turns out that meaningful simulation results were only possible if all relevant test conditions in the cutting process were considered during the experimental determination of the probability distributions of the parameter uncertainties. In the case of high speed cutting, the speed dependent dynamic behavior of the spindle is one example of such changing conditions.
Another potential application of the robust stability analysis is the adjustment of the model parameters during machining operations. The dynamic machine tool behavior, and therefore the model parameters, exhibited variation because of tool wear, thermal behavior, and other effects. Owing to the necessary standstill of the machine, it was very costly to determine the changes in the dynamic machine behavior by explicitly conducting frequency response function measurements. Instead, the additional sensors and the evaluation of the control/machine data could continuously provide information without interrupting the machining operations. In general, this database was insufficient for reconstructing the modal parameters of the dynamic behavior. However, with the approximate solution of the robust stability analysis, the reconstruction of the modal parameters may be possible. In particular, a set of initial model parameters can be provided by conducting measurements. The sensor and machine data evaluations provide information about changing signals such as noise or process forces. With regard to process stability, fast algorithms can be applied to identify the parameter changes responsible for the changes in the measured signals. Moreover, the chatter stability can be reassessed based on the new parameter set, and this will enable the adjustment of technological parameters under changing conditions.