A technique of experiment aided virtual prototyping to obtain the best spindle speed during face milling of large-size structures

The paper presents an original method concerning vibration suppression problem during milling of large-size and geometrically complicated workpieces with the use of novel way of selecting the spindle speed. This consists in repetitive simulations of the cutting process for subsequent values of the spindle speed, until the best vibration state of the workpiece is reached. An appropriate method of obtaining a computational model, called a modal approach, consists in identifying the parameters of the workpiece model created using the Finite Element Method (FEM). Thanks to the results of the identification of the modal subsystem obtained by the Experimental Modal Analysis (EMA) method, it can be stated that the parameters obtained from the experiment and delivered from the computational model have been correctly determined and constitute reliable process data for the simulation tests. The Root Mean Square (RMS) values of time domain displacements are evaluated. The efficiency of the proposed approach is evidenced by chosen technique of mechatronic design, called Experiment Aided Virtual Prototyping (EAVP). The proposed method is verified on the basis of the results of the experimental research of the relevant milling process.


Introduction
Tool-workpiece relative vibrations are the main reason of various problems occurring during a process of machining of large-size structures [1]. Some examples of these problems are reduction of the quality of the machined surface, reduction of overall machine tool productivity, increased tool wear and, in extreme cases, even the destruction of a tool or a workpiece [2]. Regeneration phenomenon, the kinematic excitation caused by repetitious immersion of cutting edges into the workpiece and the internal feedback interactions in the machine-holder-workpiece-tool system are recognized as the most important causes of vibration during milling, especially chatter vibration [3][4][5]. There are many vibration suppression and avoidance techniques developed over the years [3,6] however, many of them have been applied for research purposes only. In industrial practice, in order to avoid vibration, ''safe'' parameters (e.g. small depths of cutting, low feed and tool rotation speeds) for milling process have to be selected. Moreover, because of a lack of ability of making a lot of ''material'' experiments during production process, methods of searching for the required technology ought to be effective. Although many of the proposed approaches usually avoid chatter vibration, vibrations caused by other causes, such as an unbalanced tool or workpiece flexibility, are still present. The machining process efficiency is often limited as well. In order to increase efficiency while retaining surface quality and minimizing vibration level, including chatter avoidance, the method of spindle speed selection with the use of Experiment Aided Virtual Prototyping (EAVP) is proposed. This method is developed based on experience from the previous authors' works [7][8][9]. It belongs to the wider group of vibration reduction methods that utilize spindle speed variation [10,11], or matching the spindle speed to the selected properties of the machining process and the workpiece [12][13][14]. The advantage of these methods is that they generally need neither modifications of the machine tool structure nor utilization of sophisticated equipment, like for example: active damping methods [15][16][17][18][19], semiactive methods [20,21] or dynamic feedforward control applied towards heavy machine tools [22].
The methodology of vibration suppression through spindle speed variation met a lot of successful applications. For example, one way to reduce vibration level during both roughing and finishing industrial operations is to determine the optimal spindle speed [4,7,23,24]. Moreover, the approach reflects its convenience in practical use, because of necessary interference only with technological parameters; all conditions of the workpiece fixing remain unchangeable. A promising previous method of vibration reduction could be based on adjusting constant spindle speed to the optimal phase shift angle between subsequent tool passes [12]. The same results were obtained by the authors on a basis of the assumption, that the whole work of instantaneous cutting forces along accompanying cutting layer thickness has to be minimized [7,23]. Taking into account the influence of the dynamic properties of the workpiece on the amplitude and frequency of vibration implies the need to create a ''map'' determining the optimal spindle speeds at various points of the workpiece surface clamped on the machine tool [8,9], in accordance with the condition of minimum work of the cutting forces. However, the above approach concerns flexible workpieces having small dimensions. In case of large workpieces, the so called basic optimal spindle speed map does not meet the assumed requirements and is only the first approximation; the best spindle speed should be sought under dynamic conditions of the milling process. The other important drawback is that the approach considers only natural frequencies accompanying the identified poles of the system but omits the more important and intensive influence of harmonic frequencies of the excited vibrations. The latter means that due to more complicated nature of the large-size milling process vibration it cannot be tied only to the chatter phenomenon.
A dedicated solution described in the paper is the proposed original technique of the Experiment Aided Virtual Prototyping (EAVP). Determination of favorable states of vibration on a basis of reliable calculation models during large-size milling processes is a principal indication for solving the subjective task [25,26].
The paper was the result of a thorough reorganization and significant extension of the material presented in [27]. It is organized as follows. After the introduction containing a description of various methods of vibration reduction in machine tools, a model for computer simulation is proposed, with particular emphasis on cutting dynamics and description in hybrid coordinates. Then the original method of selecting the best spindle speed using the full FEM model of the workpiece and modal analysis is presented. The next part shows the results of simulation and experimental research on milling selected surfaces of a large-size object. The elaboration ends with conclusions and a list of contemporary bibliography.

Cutting process dynamics
Dynamic analysis of a face milling process of a flexible workpiece (Fig. 1) has been performed, based on the following assumptions [23,28].
• The spindle, together with the tool fixed in the holder, and the table with the workpiece, are separated from the machine tool structure. The remaining parts of the milling machine are recognized as ones whose influence can be neglected [4,7,29].
• Flexibility of the tool and flexibility of the workpiece are considered. The latter especially concerns a large-size flexible workpiece [30]. • Coupling Elements (CEs) are applied for modeling the cutting process dynamic interaction [23,28]. It is a phenomenological model taking into account the properties of the processed material, as well as the relationship between the instantaneous values of the components of cutting forces and the geometry of the cutting layer. Modeling using CE well reflects the physical nature of the material being processed. • An effect of first pass of the edge along cutting layer causes proportional feedback, and the effect of multiple passes causes delayed feedback additionally. The latter makes it possible to include the multiple regenerative effect in the calculation model, which is one of the potential causes of loss of tool-workpiece vibrations stability.
For instantaneous contact point between the chosen tool edge and the workpiece (idealized by CE no. l), proportional model of the cutting dynamics is included [23,28]. Based on this model, the cutting forces depend proportionally on the instant thickness of the cutting layer h l (t), as well as on the instant width of the cutting layer b l (t); both differ in time. According to the direction of the action, we separate cutting force component F yl1 acting along nominal cutting velocity v c , cutting force component F yl2 acting along cutting layer thickness, and additionally in contrast to previous approaches [23,31]-cutting force component F yl3 acting along cutting layer width (Fig. 1). The individual components of the cutting force are described by the following equations: b D (t)-desired cutting layer width, b D (t) = a p (t)/cos j r , Db l (t)-dynamic change in cutting layer width for CE no. l, h Dl (t)-desired cutting layer thickness for CE no. l, h Dl (t) % f z cosu l (t), Dh l (.)-dynamic change in cutting layer thickness for CE no. l, k dl -average dynamic specific cutting pressure for CE no. l, l l2 , l l3 -cutting force ratios for CE no. l, as quotients of forces F yl2 and F yl1 , and forces F yl3 and F yl1 , s l -timedelay between the same position of CE no. l and of CE no. l-1, a p (t)-desired depth of cutting, j r -edge angle, f z -feed per edge, u l (t)-angular position of CE no. l [10]. Zero cutting forces in Eqs. (1)-(3) correspond to the loss of contact between the cutting tool edge and the workpiece, which causes the first effect of nonlinearity of the model, due to the limitation of the cutting force characteristics. As a result, these equations describing cutting forces for CE no. l in the case of such a three-dimensional model can be represented in matrix notation: or using the abbreviated notation: where F l t ð Þ-vector of cutting forces of CE no. l,  |fflfflfflffl ffl{zfflfflfflffl ffl} dynamic change in the width of the cutting layer Db l (t). Vector (5) can also be described in six-dimensional space, i.e.,: where: The above considerations exemplify a significant progress with respect to the previous ones [7,31], because now the non-linear feedback interactions are additionally included.
In order to simplify further notation, relationship (6) takes the form: where: As a result of a milling process modeling, a hybrid system, which consists of the separated subsystems, is obtained (Fig. 1). These subsystems are: • modal subsystem, i.e., a stationary model of the Finite Element Method (FEM) of a flexible workpiece supported by Elastic-Damping Elements (EDE), which moves at the desired feed speed v f . At first, the subsystem is idealized as a set of tetragonal 10-node Finite Elements (FE) [7,23] and has a large number of degrees of freedom. However, after modal transformation [7,23], the behavior of this subsystem is described by a vector of its modal coordinates a, whose number is in practice much smaller than the corresponding number of degrees of freedom. Therefore, when we consider a finite number of subsystem normal modes mod, we define its dynamic properties using: This is also called the stiffness modal matrix; • structural subsystem, i.e., a non-stationary discrete model of a rotating milling cutter with a given spindle speed n (i.e., a flexible finite element as Euler-Bernoulli Bar (E-BB) no. e [7, 23] having a local coordinate system x e1 , x e2 , x e3 ) and the cutting process (i.e., Coupling Element (CE) no l [23,28] placed in the momentary position of the ''active'' cutting edge [23]). The cutting edges are ''active'' when they come into contact with the workpiece and the others are called ''inactive''. The subsystem's behavior is described by a vector of its generalized coordinates q. The dynamic properties of the decoupled structural subsystem (i.e., the E-BB modeling the tool itself) are determined by the matrices of inertia M, damping L and stiffness K; • abstractive connecting subsystem as a conventional contact point S between tool and workpiece. Its generalized coordinates are related to other coordinates using time-dependent constraints equations [23,24]. The latter allows us to eliminate these generalized coordinates from the description of the behavior of the hybrid system.

Dynamics of flexible details in hybrid system coordinates
Vector of deflections of CE no. l is expressed as a function of vector of generalized coordinates q and vector of modal coordinates a. Hence the following relationship is obtained [7,23]: where -vector of hybrid coordinates of the hybrid system, T l (t)-transformation matrix of displacements' vector q from coordinate system x e1 , x e2 , x e3 of E-BB no. e, to coordinate system y l1 , y l2 , y l3 of CE no. l [23,28,31], W l (t)-matrix of constraints between displacements' vector in modal coordinates a, and displacements in coordinate system y l1 , y l2 , y l3 of CE no. l [23,31]. After transformation of the vector of force interaction of CE no. l (Eq. 14) to hybrid coordinates, we obtain: As the result of the hybrid system's consideration, the matrix equation of dynamics of non-stationary model of the milling process in hybrid coordinates will have the following form [7,23,24]: where i l -number of ,,active'' coupling elements [23,28].
Simulations in the time domain, based on Eq. (19), take into account the most important non-linear effects observed in real milling operations, i.e., the relative vibration of the tool and the workpiece along the thickness and width of the cutting layer, loss of contact of the cutting edge with the workpiece, multiple regenerative effect, through successive change of the geometrical position of the momentary connection of the tool with the workpiece. In this way they demonstrate the advantages of performing time simulations instead of the traditional linear stability analysis.
In order to identify the modal model of the flexible workpiece [which is part of Eq. (19)], the matrix of normal modes W and the matrix of corresponding natural angular frequencies X of the modal subsystem must be determined. The separation of the modal subsystem from the entire non-stationary structure allows us to reduce the finite element model to several modes, the number of which depends on the importance and the need to select modes for further analysis. As a result, the size of the model is significantly reduced.
Normal modes are unchanged over time because it is assumed that the modal subsystem remains stationary during the cutting process. For this reason, normal modes W and angular frequencies X can be identified by: • computer software for calculating eigenfrequencies and corresponding normal modes of discretely For example, impact modal tests are performed on the workpiece actually installed on the machine tool's table.
Both approaches are recommended due to the need for mutual verification of the results obtained.
The number of modal coordinates is generally much lower than the number of generalized FEM coordinates, therefore calculations should be performed much more efficiently. The dominant vibration modes can be different at different places of the manufactured surface, especially for workpieces of large size and complicated shapes. Therefore, in the case of a modal model computed analytically using the FEM model, it is usually most convenient to use the set of the first few elastic modes. It is assumed that the lower order modes are more important in terms of quality and quantitative properties of the cutting process simulation. However, in the case of the modal model obtained directly from the experiment, it is assumed that the modes identified on the basis of the Frequency Response Function (FRF) measured only at driving points lying on the manufactured surface dominate.

Selection of the best spindle speed based on simulation of cutting process and full FEM model of the workpiece
In the proposed technique of Experiment Aided Virtual Prototyping, based on the modal approach, the following procedure of spindle speed selection is suggested (Fig. 2).
1. Creation of the FEM model of flexible workpiece including its support modelled as a set of Elastic-Damping Elements (EDEs), and computation of normal modes. 2. Modal parameters' estimation with the use of the EMA techniques.

Verification of the FEM model using the Modal
Assurance Criterion (MAC) [32].
where W a a -vector of normal mode no. a calculated from the FEM model, W e b -vector of normal mode no. b identified by modal tests. 4. Once the FEM model of the whole workpiece is assumed to be valid, dominant normal modes are selected in order to be used in modal subsystem of the hybrid model of the milling process. 5. Determination of cutting process simulation parameters (k dl , l l2 , l l3 ) by conforming milling simulations performed for nominal (standard) spindle speed of the process to the results of real cutting performed for the same speed. Obtaining the above compliance condition is based on a comparison of the Root Mean Square (RMS) values of tool-workpiece relative displacements, obtained from simulation and the implementation of the cutting process. 6. Simulations of the milling process for various spindle speeds from the selected range. Cutting process parameters determined in p. 5, are applied. For a desired set of the simulated milling processes, vibration levels are observed (for example, the RMS value of tool-workpiece relative displacements), and as the result, the best spindle speed is selected. 8. Real milling with the best spindle speed selected.
The advantage of this method is that the exact model of the entire workpiece and supports is used for simulation. This can be especially important in the case of complicated part geometry. On the other hand, however, using the full model allows us to specify all normal modes that are necessary to model the dynamic behavior of the entire workpiece. It may include modes that are not really relevant to modeling machined surface behavior. This can unnecessarily complicate the model and thus extend the simulation time.
An important feature of the developed method based on EAVP technique, from the point of view of prospective industrial applications, is that in order to select the best spindle speed, the required ''material'' experiments of the actual process were carried out only twice, i.e., for standard process parameters and improved.

The workpiece
The experimental research concerned investigating the dynamic behavior of a large workpiece (total size 2061 9 1116 9 540 mm) made of STW22 03M steel (Fig. 3a) and clamped on a table of the MIKROMAT 20 V portal milling center. The workpiece was selected from the common production program of one cooperating industrial company. For a purpose of the research the following measurement equipment setup was applied [33]: 15 DJB A/120 V accelerometers with measurement range ± 75 g; PCB 086C03 modal hammer, range ± 2224 N; National Instruments PXI-8106 controller with NI PXI 4496 24-bit simultaneous DAQ card working under the LabView 2016 RT environment. All signals were sampled with frequency at least 5 kHz (during cutting experiments) or 15 kHz (during modal tests). Measurements were performed using custom developed software. There were two surfaces milled (Fig. 3b). For surface 1, full face milling was performed first by the tool moving from the left (i.e., starting from a vicinity of accelerometer 22) to the right. Down milling was performed next by the tool moving in the opposite direction (i.e., starting from a vicinity of accelerometer 25). These two passes formed one complete operation. Length of the surface was 1778 mm, and width-57.5 mm. The milling was performed by the 4-edge Sandvik face milling cutter R390-044C4-11M060, /44 mm, cutting edge angle j = 90°. For surface 2 only one pass of down milling was performed by the tool moving from the left (i.e., starting from a vicinity of accelerometer 18). Length of the surface was 1789 mm, width-55 mm. Milling was performed by the 11-edge Sandvik face milling cutter R390-125Q40-17H, /125 mm, cutting edge angle j = 90°. The use of tools with cutting edge angles j = 90 o does not require taking into account in the calculation model the phenomenon of tooling system bending, typical for relatively large cutter diameters [34].
The measuring points were selected in such a way as to be able to record measuring signals mainly along the milled surfaces. Due to the limited number of accelerometers in experimental research, there were 9 of them, i.e. 5 on the first surface and 4 on the second, placed at equal intervals, and the remaining-in other less important places of the workpiece (Fig. 3b).

Standard parameters
The standard parameters used for milling process in the common production scheme of cooperating industry, in case of the presented workpiece were: • for surface 1n = 1300 rev/min, v f = 600 mm/ min, a p = 1 mm; • for surface 2n = 560 rev/min, v f = 1233 mm/ min, a p = 1 mm.
The vibration level during milling with the above mentioned parameters is treated as a reference to compare the results of the method used to search for the best spindle speed. It is worth noting that the vibrations during milling with standard parameters were generally very low, so further reduction of vibrations was a significant challenge of the research.

Modal identification and spindle speed selection
The experimental modal tests of the workpiece were performed. 15 global modes were identified with the use of Eigenvalue Realisation Algorithm (ERA) [23,35]. The FEM model of the workpiece and its support was tuned to correlate it with the identified modes. With the use of appropriate FEM software i.e., PERMAS and MEDINA [23], calculation model of the workpiece, composed of 12 511 156 nodes, 7 764 895 10-node tetragonal finite elements, 2158 bar elements (RBAR) and 423 rigid bodies (RBE), was obtained. 15 natural frequencies and accompanying normal modes of the workpiece mounted in the supports were calculated. 9 of the modes identified by the ERA conformed with modes calculated by the FEM (MAC values are over 0.8 and the corresponding computed frequencies are in bold, see Table 1). For these modes modal damping coefficients were selected in accordance with the identification results. For the rest of the modes, damping coefficients were estimated (italic letters in Table 1). Approximate time to obtain the expected compatibility of the modal model of the workpiece using a laptop computer equipped with an Intel Ò Core TM i7-6700HQ CPU 2.60 GHz processor and having RAM 32 GB, was about 240 min. The vectors of mass-normalized normal modes for points lying on the machined surface were then exported to the original AMIKRO4 cutting process simulation software running in the 64-bit MSYS2 MinGW environment on the same type of portable computer. Simulations were performed for a set of the defined spindle speeds using FEM model based on 15 modes having modal parameters as are shown in columns 4 and 5 ( Table 1). The above means that both the FEM model and the modal models were constructed and then subjected to experimental validation in the real range of interesting frequencies of 500 Hz. Subsequent step of the modal procedure of searching for the best spindle speed using EAVP when milling large-size workpieces (Fig. 2) is to create nonstationary calculation model of the face milling process of the indicated surfaces. The reason is that a compliance ought to be obtained between the RMS values of simulated vibrations' displacements and the measured ones during standard production conditions. Of course, one can also consider using other quantities to validate milling process simulation models. However, it should be noted that the simulated models are non-stationary and strongly non-linear. Therefore, other conformity assessment measures, based e.g. on the frequency analysis of steady vibration states, do not apply here. RMS is the best measure for assessing displacements in the examined issue.
The corresponding RMS values of simulated plots, adjusted to milling at standard parameters (see p. 4.2) were: for full milling of surface 1-0.000287 mm, for down milling of surface 2-0.002971 mm. As a result of the selection, the following parameters of the hybrid model are adapted to milling simulations to meet the RMS values compliance condition: for surface 1, k dl = 500 daN/mm 2 , l l2 = 0.4, l l3 = 0.42; and for surface 2, k dl = 500 daN/mm 2 , l l2 = 0.4, l l3 = 0.24. The approximate time of determining values of these parameters was: for surface 1-5 min., for surface 2-4 min. The speed range was 1300…1500 rev/min for surface 1, and 550…800 rev/min for surface 2. The lower limit is the standard spindle speed used for milling the selected workpiece, and the upper limit is the maximum allowable speed for the selected tool and the maximum cutting speed for the insert. Spindle speeds below 1300 and 550 rev/min were not simulated as the intention of the performed research was not only to reduce vibration level but also to increase process efficiency. For every spindle speed, vibration level was observed and three indices were calculated, i.e.,: • RMS 95% -RMS of tool-workpiece relative displacements calculated for 95% of the whole cutting time. The remaining 5% consist of transient effects of the tool entrance and exit out of the workpiece; • A max -maximum amplitude of tool-workpiece relative displacements calculated for the same period as for RMS 95% ; • RMS 95%MR -RMS of tool-workpiece relative displacements calculated for 95% of the whole cutting time (similarly to RMS 95% ), but relatively to the mean value of the considered vibrations (MRmean related). Thus, this indicator is equivalent to standard deviation of vibrations. This can be  Table 3) interpreted as an indicator of the vibration level after correcting the surface displacement caused by tool deflection in the workpiece. This corresponds better to the way how vibrations are measured during the real milling process, because during acceleration measurements ''static'' components of the machined surface deflection are omitted.
The values of the abovementioned indices for simulated spindle speeds are presented in appropriate figures. According to the simulation results, spindle speed n = 1380 rev/min (v f = 637 mm/min) was selected for surface 1 as the best (Fig. 4) and was applied for milling process. Two exemplary plots of simulation results are presented in Fig. 5. For surface 2 as the best spindle speed n = 700 rev/min (v f =-1541 mm/min) was selected (Fig. 6) and applied for milling. Additionally, n = 800 rev/min was also applied for milling in order to present an alternative spindle speed, as it was shown by simulations.   Figs. 4, 6). It allows to determine appropriate conditions at which the workpiece vibration level approaches minimum.

Real milling results
Milling operations were performed for both workpiece surfaces and various spindle speeds selected according to standard parameters and the abovementioned method. In the Table 2 selected spindle and feed speeds are presented. In Tables 2, 3 and 4 sign ''Axx'' denotes given accelerometer's number (see, Fig. 3b). In the Table 3, RMS values of displacements for Fig. 7 Vibrations of workpiece during full milling of surface 1 at standard parameters (see Table 2) Fig. 8 Vibrations of workpiece during full milling of surface 1, modal approach (see Table 2) performed milling operations observed in measurement points during passing the tool over the given sensor, are presented. For example, value for sensor A22 was calculated for time period from the beginning of the milling operation (excluding the time period while the tool was entering into the material) to the instant of time when the tool was in the middle between positions of sensor A22 and A23 (see, Figs. 7, 8 for an illustrative example and Fig. 3b for sensors positions). Displacements values were calculated by double integration of the measured accelerations. Table 4 presents the same data but as the relative values. The latter helps to notice that the approach provides better results towards vibration suppression. Vibration reduction is denoted with the ''-''sign. First, milling with the standard parameters was performed and its results are recognized as the reference for subsequent tests. The presented methods of the best spindle speed selection was examined. Chosen results of the measurement are shown (Figs. 7, 8).
The application of the modal approach to milling surface 1 allows EAVP to choose the best spindle speed for which the simulation results are overestimated in relation to the measured average displacements (Fig. 4, Table 3), and for surface 2-the results of down milling are somewhat underestimated (Fig. 6, Table 3). There is currently no explanation for this phenomenon. Most important, however, is that in reality this approach gave milling results better than standard types.

Conclusions
The efficiency of the proposed method of vibration suppression has been proven by selecting the best spindle speed when milling large-size structures using a dedicated technique of Experiment Aided Virtual Prototyping (EAVP). Thanks to the results of the modal approach, and thus-the identification of the modal subsystem obtained by the ERA method, it can be concluded that the parameters obtained from the experiment and provided from the calculation model were correctly determined and constitute reliable process data for simulation tests. This is confirmed by 9 natural frequencies identified for the spectrum band up to 500 Hz.
In general, the proposed method ensures and also allows to predict-satisfactory results of the reduction of tool and workpiece vibration, in some cases up to 50%, especially in the case of down milling operation. The advantage of the modal approach is the accuracy of the results of the computer prediction of the best spindle speed, which makes them attractive for practical applications in an industrial environment. However, creating a full FEM model in this approach means spending a lot of time both building the model and adapting it to the results of experimental modal tests. In the case of larger or more complex components, it can take many hours of labor, and thusgenerate unnecessary economic follow-ups. Tests are limited to treated surfaces only. Nevertheless, reliable dynamic behavior of the machined surface will be identified when all modes affecting a particular surface are observed during this limited modal test.