Application of interpolation methods for the determination of position-dependent frequency response functions for the simulation of 5-axis milling processes

The occurrence of chatter vibrations in 5-axis milling processes is a common problem and can result in part failure, surface defects and increased wear of the cutting tool and the machine tool. In order to prevent process vibrations, machining processes can be optimized by utilizing geometric physically-based simulation systems. Since the modal parameters of the machine tool are dependent on the position of the linear and rotary axes, the dynamic behavior of milling processes can change along the NC path despite constant engagement conditions. In order to model the pose-dependent modal properties at the tool tip, the frequency response functions (FRFs) were measured at different locations of the workspace of the machine tool for various poses of the rotary axis of the spindle. To take the varying compliance within the workspace of a machine tool into account in a geometric physically-based milling process simulation, different interpolation methods for interpolating FRFs or parameter values of oscillator-based compliance models (OPV) were applied. For validation, the resulting models were analyzed and compared to measured data. In OPV interpolation, the individual oscillation modes were interpolated in their respective characteristics based on the oscillator parameters (eigenfrequencies, modal masses and damping values). In FRF interpolation, however, there was no differentiation between the modes, resulting in a wrong interpolation. It can therefore provide good results when only a small shift of the eigenfrequencies is expected, as in case of the analyzed machine tool, with only small movements of the translatory axes.


Introduction
Chatter vibrations in milling processes can affect, e.g., the surface quality, increase the wear of the machine tool components and decrease the tool life [1,2]. In order to achieve stable milling processes, several approaches for analyzing and optimizing the machining process can be utilized [3]. For instance, analytical or geometric phyically-based simulation systems were used for the optimization of the process parameter values concerning the dynamic behavior of the tool-spindle-machine system [4][5][6].
The dynamic properties of machining centers change when the position of the rotary and linear axes is varied [7]. Hung and Lin, e.g., showed the influence of the position of a bi-rotating milling head on the process stability [8]. Furthermore, the dynamic behavior during machining processes is influenced by the condition of the spindle [9]. In fiveaxis milling, the finishing process is often conducted using spherical or toroidal cutters. During the finishing process, the tilt and lead angle of the rotary axes change constantly, which directly affects the inner forces of the bearings of the spindle [10], which are often ball bearings, e.g., in highspeed machining centers [11]. As a consequence, the axial preloading of the ball bearings influences the radial stiffness of the ball bearings [12]. Lee and Altintas, as well as Kono et al. determined an alteration of the compliance magnitude of up to 40 % resulting from a variation of the tilt angle [10,13]. Moreover, the position of the linear axes can also influence the dynamic properties of the machining center [7]. 1 3 In order to identify the modal properties of a toolspindle-machine system at the tool tip, impact hammer tests can be conducted [14]. On the basis of the measured impact force and the resulting response of the system, FRFs are determined which describe the dynamic behavior of the tool center point (TCP) at the measurement point (MP) in the workspace of the machine tool. Based on the FRFs, stability lobe diagrams can be determined in order to identify axial depths of cut and spindle speeds leading to stable and unstable processes [15,16]. For the determination of these diagrams, compliance models based on uncoupled, damped harmonic oscillators can be used in geometric physically-based simulations (GPS) [17,18] to model the dynamic compliance in machining processes.
To consider the position-dependent dynamic behavior at the tool tip, the FRFs have to be determined for different locations in the workspace of the machine tool and different positions of the rotary axes, respectively tilt and lead angles. Furthermore, a multidimensional compliance model is needed in order to consider not only the influence of the linear position [19] of the tool tip, but also the influence of a spatial movement of the tool in the workspace. For practical reasons, the measurement of FRFs at the tool tip is only reasonable for a limited number of different poses. The modal parameter values for the remaining poses on an NC path could be approximated by an interpolation method.
In previous studies, for instance, the barycentric interpolation method was used to determine the modal parameter values for the simulation of workpiece [20] and tool vibrations [19]. In [19], the dynamic behavior at the TCP was modeled describing a movement of the linear (X,Y,Z) but not the rotary axes. In [21] also the movement of the rotary axes was considered by utilizing a Machine Learning (ML) approach and the linear interpolation for the prediction of FRFs and the oscillator parameter values of a compliance model, respectively. In order to calculate the FRFs at the tool tip and determine compliance models for several poses in the workspace, the interpolation and prediction of FRFs as well as of modal parameter values of the oscillator based compliance models (OPV) were analyzed. Compared to the ML approach, the linear interpolation approach was easier to implement and required a shorter computing time. Resulting FRFs of both, the ML approach and the linear interpolation, were compared to measured FRFs by the least-squares error. The results of the linear interpolation were less accurate by 25 % and 31 % compared to the ML approach. However, the deviation occured primarily for the modes with low compliance amplitudes. The main eigenfrequencies with a high compliance amplitude could be modeled with both, the ML and the linear approach. Hence, both approaches were suitable for simulating chatter vibrations and calculating stability lobe diagrams.
In this paper, the linear interpolation is further investigated using different approaches in terms of interpolating both, the FRFs and OPV, as in [21]. In the first approach, FRFs are interpolated and oscillator-based compliance models are fitted according to the interpolated FRFs. In the second approach, the OPV eigenfrequency, damping value and modal mass are interpolated. Therefore, OPV previously needed to be calibrated based on the FRFs measured at selected MPs. In this approach, a system consisting of coupled components, i.e., cutting tool, spindle, guides and machine structure, is represented by independent oscillators whose modal properties are interpolated. Since the parameterization of the OPV is not distinct, several different ways exist to model a similar or even the same FRF. Thus, even if the FRFs were similar between the MPs, there could be a high divergence between the OPV. Nevertheless, this procedure could save computing time in contrast to the FRF interpolation and subsequent parameterization of OPV for each simulation step. Therefore, it was considered and evaluated by means of the interpolation of the axis position-dependent dynamic properties of an exemplarily selected machining center. For the analyzed machine tool, three dimensions were taken into account for the interpolation methods, since two linear axes and the tilt angle of a fork head spindle needed to be considered for calculating the stability limit. A GPS for modeling milling processes taking into account position-dependent modal properties is described in Sect. 2. In Sect. 3, the experimental procedure and experimental results are presented. Moreover, the principles of three different interpolation methods are introduced. The results of the interpolation of the FRFs or, respectively, OPV, are presented in Sect. 4. Furthermore, the GPS is used in order to calculate stability lobe diagrams on the basis of the interpolated and the measured FRFs as well as the interpolated OPV of the compliance models. The influence of the interpolation methods and approaches on the resulting stability lobe diagrams are discussed in Sect. 5. In Sect. 6, a conclusion and an outlook are given.

Simulation of milling processes taking position-dependent dynamic properties into account
In this section an overview of the position-dependence of the dynamic properties of machining centers and the simulation of milling processes using a GPS is given. A focus is set on the compliance models. In order to analyze the dynamic behavior of milling processes and to identify the stability limit, a GPS system, developed at TU Dortmund University, was used [22]. GPSs are time-discrete simulations for determining, e.g., process forces, process dynamics and resulting surface location errors in milling processes by integrating and combining different physically-based models with a geometric process model [5,23]. Geometric models based on, e.g., the Constructive Solid Geometry (CSG) technique [24] can be used to describe the tool and workpiece and, thus, to calculate the uncut chip shape for each tooth feed [22]. The uncut chip shape is scanned with a ray-based model representing the position of the cutting edges and calculating the uncut chip thickness in refined discrete time steps [25,26]. The process forces resulting from the tooth engagements are calculated using the empirical cutting force model by Kienzle [27]. This force model needs to be individually parameterized for the applied combination of cutting tool, material and tilt angle [18]. The specific cutting force parameter values used in the process simulations are shown in Table 1. Process forces can lead to deflections of the compliant system consisting of the machine tool, spindle, cutting tool and workpiece. These vibrations during the milling process cause a modulation of the chip thickness. This leads to a recurring excitation of the dynamic system and, thus, can cause regenerative chatter [16].
For the determination of the dynamic behavior of a toolspindle-machine system, the OPV eigenfrequency f 0 , modal mass m and damping coefficient were fitted using an evolutionary optimization algorithm to reach the least sum of squared differences between the magnitude functions of the measured and modeled FRFs [18]. As the OPV of the compliance model were fitted by an evolutionary algorithm [18], multiple, equivalent solutions can be found instead of a globally optimal set of OPV. Hence, several different combinations of the OPV can be used to model the FRFs. Therefore, diverging modal masses or damping coefficients can occur for similar FRFs of different MPs.

Interpolation methods for FRFS and OPV
The pose-dependent FRFs were measured in impact hammer tests for several tilt angles at discrete MPs in the workspace of the analyzed machine tool. Oscillator-based compliance models were fitted afterwards. Different methods for the interpolation of either FRFs or OPV are presented in the following sections.

Experimental procedure for the determination of pose-dependent FRFs
In this study, the impact hammer tests were conducted on a five-axis CNC machining center DMG Mori HSC 75 linear ( Fig. 1) with an impact hammer Type 8202 (Brüel&Kjaer). The impulse was introduced at the tip of a spherical end mill Fraisa X7400 with a cutting edge number of z =2, a diameter of d = 10 mm and a cantilever length of d 1 = 30.4 mm. The resulting acceleration was measured using an accelerometer (PCB 352C23). The machine tool was equipped with a fork-head kinematic spindle, which provides two linear axes, Y and Z, and one rotary axis, B, on the spindle side. Measurements of the FRFs were conducted in 48 poses representing different positions of the B-axis at 17 MPs. The OPV were determined on the basis of the measured FRFs by parameterizing the compliance models, as described in Sect. 2. In order to keep the experimental effort as low as possible, the position of the B-axis was varied with the angles +10°, 0°, −10°, −20°, −45° and −60° only at the selected MPs. The position of the X-and C-axis are related to the machine table and are, therefore, not relevant for the compliance of the tool-spindle-machine system analyzed in this work.

Calculated stability limits
In order to analyze the effect of the pose on the process stability, stability limits were calculated for different poses of the machine tool axes using a geometric physically-based process simulation system as described in Sect. 2. For this purpose, compliance models were parameterised based on the measured FRFs for three exemplary poses. In order to take into account the variation of the effective cutting speed and cutting edge geometry along the cutting edges of the spherical end mill, the tilt angle-dependent model shown in Table 1 was used to calculate the process forces. As described in [29], a shift of the eigenfrequencies ▵ f 0 results in a shift of the asymptotes of the respective stability lobes where z is the number of cutting edges. Thus, for an end mill with z = 2 cutting edges a shift of the eigenfrequency of ▵ f 0 = 10 Hz already results in a shift of the asymptote of the zeroth stability lobe of ▵ n a,0 = 300 1/min . Thus, the pose-dependency of the eigenfrequencies shown for different tilt angles of the B-axis in Fig. 2 led to a shift of the stability lobes for the different poses (cf. Fig. 3).
Comparing the stability limits for the poses G0 and Q0, the lobes were particularly shifted in the low spindle speed range. If not modeled, this position-dependency could lead to process parameter values being chosen which result in an unstable process. For the pose G-30, an increase of the stability limit by a factor of about two was calculated (Fig. 3). This can be explained by the reduction of the specific process forces, especially, the specific cutting normal force, cf. Table 1. This shows that the modeling of pose-dependent dynamic characteristics is necessary for the prediction of stability limits for machine tools with tendentially low and high pose dependence of the dynamic properties.

Utilized interpolation methods
In this subsection, different interpolation methods for the interpolation of FRFs or OPV are introduced. The interpolation of OPV is only reasonable if each compliance model represents the same eigenmodes, even if their compliance amplitude becomes negligibly low in some poses. By unambiguously indexing the oscillators according to the corresponding eigenmodes, a misallocation of the oscillators and, thus, an interpolation of oscillators representing different eigenmodes can be prevented. This is especially important if the eigenfrequencies of multiple eigenmodes cross when changing the pose. A simple interpolation method utilized for the interpolation of the OPV or FRFs was the Nearest Neighbor Interpolation (NNI). It assigns the value of the closest known point, i.e. MP, to the interpolation point [30]. Due to different distances between the MPs on the different axes, the closest MP is determined as shown in Fig. 4. After the determination of the nearest MP, the dynamic behavior of the nearest neighbor is selected, i.e., the OPV or the FRF was adopted for the interpolation point R, respectively.
To keep the experimental effort as low as possible, every MP was investigated for B = 0°, but the variation of the angle of the B-axis (B ≠ 0°) was only conducted for selected MPs. In order to take the influence of the B-axis position on  the modal parameters into account, the NNI was conducted in one step if B = 0° and in two steps in case B ≠ 0°. If B = 0°, the closest MP was identified calculating the shortest direct distance using the Y-an Z-coordinates of the MP (see Fig. 4). If the required tilt angle was B ≠ 0°, the algorithm monitors whether measured data for different tilt angles were available at the nearest MP or not. If no data for B ≠ 0° could be identified, the model searched for a closest MP that contained compliance models and FRFs for B ≠ 0°. Afterwards, the dynamic behavior of the determined pose was selected for the required axis-positions (B, Y, Z). Hence, when planning the measuring points and poses, it should be taken into account that an uneven distribution of the examined poses can lead to the choice of a distant measurement in this interpolation method.
The barycentric interpolation was already used in previous studies in order to calculate the compliance of a workpiece [20] and tool [19] at different positions. This method uses a tetrahedral approach for the interpolation of the OPV or FRFs. The compliance model for the current position R (Fig. 5) is interpolated from the nearest MPs. Depending on the MP distribution, two or three MPs are sufficient, if the interpolation point is positioned on a straight line or a triangle connecting the MPs in the 3-dimensional search space, respectively. In case the distribution of the MPs does not provide a straight line or triangle as a connection, a tetrahedron with the vertices N, O, P, Q which encloses the interpolation point is used. The coordinates of the measurement and interpolation points are given in generalized barycentric coordinates. Thus, non-equidistantly distributed measuring points on the various axes of the machine do not cause any weighting of the individual axis coordinates. Figure 5 shows the principle of the barycentric interpolation for the case a tetrahedron is used, while equations 1-3 present the calculation of the eigenfrequency f 0 , modal mass m and damping coefficient of each individual oscillator for the interpolation point R.

The Weighted Nearest Neighbor Interpolation (WNNI) is a combined method of the linear interpolation and the NNI.
In order to calculate the OPV or the FRF of the required pose (RP), the four nearest MPs enclosing the RP are identified. The OPV and FRFs of the selected MPs are weighted according to their distance from the interpolation point. Therefore, the individual distances from the interpolation point to the MPs are summed up and the distance for each MP i is divided by the total distance to determine the weighting factor W i , shown in Eqs. (4) and (5): For the interpolation based on the tilt angle, the linear interpolation is conducted. If the required tilt angle B r ≠ 0 and nearest MPs do not contain data for B ≠ 0, the closest MPs providing FRFs for B ≠ 0 are chosen. First, for each of the four nearest MPs, the FRFs or OPV are interpolated for the required tilt angle B r , as shown in Fig. 6. Afterwards the interpolation is conducted for the position of the interpolation point R.The FRF is defined as the compliance amplitude |H| for each frequency step f j in the considered frequency range, see Eq. (6). Here, |H(f j )| is the compliance amplitude at the point f j and f min and f max are the lower and upper limits of the considered frequency interval. The interpolation of the FRFs at the involved MPs is shown in Eq. (7). For the required FRF, the interpolated data from the MPs are weighted as shown in Eq. (8). The interpolation of the OPV is conducted analogously to Eqs. (7) and (8).
B u and B o represent the two B-axis positions where measurement data exist and the required B-axis position B r is located in between. The compliance amplitude |H| at the MP i is interpolated over the measured frequency range in order to determine the FRF at B r . In Fig. 6, B u is represented by B 0 and B o is represented by B 1 . It can be suitable to not interpolate all OPV, eigenfrequency f 0 , damping coefficient and modal mass m, especially, when a parameter varies significantly due to the evolutionary optimization algorithm. In this study, the eigenfrequency and damping coefficient were interpolated by the described approach. The modal mass was selected from the nearest neighbor, since it is not changed with a pose variation, but only displaced. This procedure can produce faulty models, if the OPV have strongly deviating parameters. In the present OPV, the eigenfrequencies and damping coefficients were similiar, while the modal masses deviated by a factor of about 10 due to the used OPV fitting method. Hence, the deviation caused by the interpolation of all three modal parameters would have caused a larger error than a possibly unfavorably selected modal mass from the nearest neighbor.

Interpolation of position-dependent modal properties
In order to analyze the suitability of the investigated interpolation methods, the results of interpolating FRFs and OPV were validated for exemplary poses of the machine tool (Table 2). In Figs. 7 and 8, the results for both, interpolated FRFs and compliance models with interpolated OPV, are compared for the interpolation methods presented in Sect. 3. The linear interpolation of the FRFs (Fig. 7) led to a lower deviation from the measured FRFs than the NNI in terms of both, the compliance amplitude and the phase shift. Especially the important eigenmodes with a high compliance amplitude, which typically limit the stability boundary, were interpolated with a higher accuracy. Since the measured FRFs had a noisy phase shift, the FRFs calculated with both interpolation methods also exhibited this behavior. Subsequent  In terms of the interpolation of compliance models (Fig. 8), the linear interpolation and barycentric interpolation of the OPV led to a better correlation with the measured FRFs in the frequency areas with high compliance compared to the NNI and WNNI. All interpolation methods led to differences in the FRFs at the eigenfrequency f 0 = 2070 Hz in X-direction as well as between f = 2420 Hz and f = 2950 Hz in Y-direction. This can be explained by the fact that fixed numbers of oscillators were used when parameterizing the oscillator models for the MPs used. A number of six oscillators in X-and four oscillators in Y-direction were modeled which represented the eigenmodes considered decisive for process stability owing to their high compliance amplitudes. In contrast, the eigenmodes at f 0 = 2420 Hz (X-direction) and f 0 = 2950 Hz (Y-direction) were not modeled since they were not considered critical for the calculation of process stability due to their comparatively low compliance amplitude.
Using the NNI, the selected OPV in X-direction resulted in a higher amplitude at the eigenfrequency f 0 = 1500 Hz due to the high pose influence of this eigenmode (c.f. Fig. 2). In Y-direction, a good correlation with the measured data could be achieved, as a low pose influence was identified between the interpolation point and the nearest measurement point.
The FRFs calculated by WNNI of the OPV showed a representation of the eigenmode at f 0 = 1500 Hz in X-direction due to an interpolation of the damping value. The significantly higher compliance amplitude at f 0 = 1320 Hz in Y-direction can be explained by the interpolation of the compliance values and selection of the modal mass of the nearest neighbor. A different parameterization of the FRFs at the surrounding MPs due to the different characteristics of their eigenmodes resulted in this deviation.
For the barycentric interpolation of the compliance models for the poses G0 and G-30, two measuring points G+10 and G-10, respectively, G-20 and G-45 were identified, on which connecting straights the interpolation points were located. Thus, the OPV for these poses could be interpolated using two known compliance models each. For the point G0, four MPs which span a tetrahedron around the interpolation point were selected.
The interpolated OPV of the compliance models resulting for these poses are given in Table 3. Between the varied poses in both directions there were differences particularly in the eigenfrequencies f 0 between f = 1400 Hz and f = 1600 Hz. Furthermore, the deviation of the damping factors and modal masses m of the oscillators between different MPs was higher than between different B-axis poses at one MP (compare Q0 and G0 vs. G0 and G-30). As shown in Fig. 8 for pose G0, these interpolated compliance models showed a good correlation with the measured FRFs in terms  of the compliance amplitude. The phase shift showed a high deviation between the measured FRFs and the interpolated FRFs which arised from the low number of oscillators used in the initial parameterization of the compliance models at the MPs [18].
. error measure, the error of the eigenfrequency with the highest compliance amplitude is presented. For the barycentric interpolation as well as the NNI and WNNI of the OPV, the mean error values were higher than for the linear interpolation and NNI of the FRFs, which resulted from the fitting of the OPV using a

Simulation of milling processes considering the position-dependent dynamic properties
In order to investigate the influence of the changed modal properties of the system and of the interpolation method, stability lobe diagrams were calculated using compliance models based on the measured and differently interpolated FRFs as well as directly interpolated compliance models. Therefore, a geometric physically-based milling process simulation as described in Sect. 2 was used. A spindle speed range of n = 3000 1/min -6000 1/min and a depth of cut of a p = 0.01 mm -0.80 mm, which is reasonable for machining the material ASP 2012 (1.3397), was investigated. Representative calculated stability limits for the poses Q0, G0 and G-30 are shown in Fig. 10. The compliance models based on the measured FRFs and the interpolated compliance models were utilized to calculate stability limits. Afterwards, the resulting stability diagrams were compared and the results were evaluated for the different interpolation approaches and methods.
Using the barycentric interpolation of the compliance model, the calculated stability limits showed a good correlation with the references, which were calculated with compliance models based on the measured FRFs. This confirms a reasonable interpolation of the eigenmodes crucial for process stability, i.e., the eigenmodes with the highest compliance amplitudes (cf. Fig. 9).
The stability limit calculated for the FRFs interpolated by WNNI from its four nearest neighbors also closely matched the references. The used interpolated FRFs also showed a good correlation of main eigenfrequency to the measured FRF and a low mean error of the compliance amplitude (cf. Fig. 9). The deviation of the phase shift was comparable to the barycentric interpolation of the OPV.
The interpolation of the compliance models from the four nearest neighbors using WNNI led to a shift of the stability lobes in view of the spindle speed. This shift resulted from a significantly higher error of the eigenfrequency with the highest compliance amplitude compared to the other interpolation methods (cf. Fig. 9). For the pose G0, the calculated stability limit was slightly lower (Fig. 10b)), since the compliance amplitude of the interpolated compliance model was higher (cf. Fig. 8).

Conclusion and outlook
In this paper, two general approaches for the linear interpolation and NNI of the dynamic behavior of 5-axis machine tools are presented. The linear interpolation of FRFs and the interpolation of the OPV of oscillator-based compliance models were compared. The calculated FRFs were validated by comparing them with measured FRFs for three exemplary poses in the workspace of a machining center. Moreover, a geometric physically-based process simulation was used to determine stability lobe diagrams for each interpolation method and each approach for exemplary poses. The barycentric interpolation method as well as the WNNI of the OPV showed a good correlation with the measured data regarding the calculated FRFs as well as the stability limit. The results were comparable to an ML-based prediction method presented in [21] in terms of the error of the mean compliance amplitude and phase shift as well as the calculation of the eigenfrequency with the highest compliance amplitude.
For the NNI, a high deviation between the measured data and the calculated data can occur in case that the MPs show a strongly different dynamic behavior. The direct interpolation of the FRFs using WNNI also resulted in a good agreement between the interpolated FRFs and calculated stability boundaries with the measured reference FRFs and stability lobe diagrams, respectively.
A disadvantage of the interpolation of the FRFs is the need of a subsequent parameterization of the compliance model for each time step of the GPS while the pose is continuously varying. Thus, this parameterization must be automated and requires a significantly higher computational effort than the OPV interpolation. Hence, it can cause high calculation times and can be an additional failure source. In contrast, in the interpolation of the OPV, the oscillator fitting is conducted in a preprocessing step for a limited number of MPs, which allows a manual validation of the OPV, especially of the eigenfrequencies with the highest compliance amplitudes.
In this study, the different interpolation methods used to determine position-dependent FRFs were validated based on measurements conducted on a compact 5-axis machining center. To validate the interpolation methods for FRFs with a higher position-dependency, future studies should also incorporate experiments conducted on machine tools with a larger workspace and different axis configurations and, thus, a higher variation of the dynamic behavior among the workspace.