Experimental and Numerical Dynamic Identiﬁcation of an Aerostatic Electro-Spindle

Featured Application: Aerostatic spindles are used in high-speed micromachining applications. The of the work to validate the developed non-linear numerical model through the proposed identiﬁcation procedure and the performed experimental tests. Abstract: This paper proposes a method to experimentally identify the main modal parameters, i.e., natural frequencies and damping ratios, of an aerostatic spindle for printed board circuit drilling. A variety of methods is applied to the impulse-response function of the spindle in the presence of zero rotational speed and different supply pressures. Moreover, the paper describes the nonlinear numerical model of the spindle, which consists of a four-degree-of-freedom (DOF) rigid and unsymmetrical rotor supported by two aerostatic bearings. The main goal of the work is to validate the developed non-linear numerical model through the proposed identiﬁcation procedure and the performed experimental tests. The comparison proves satisfactory, and the possible sources of uncertainty are conjectured. This proposes a cheap and simple experimental identiﬁcation procedure for aerostatic spindles. Impulse and static tests are performed to identify the natural frequencies, damping ratio and static stiffness evaluated at the nose of an electro-spindle for PCB drilling. The natural frequencies and damping factors are measured by considering the impulse-response function of the system, whereas the stiffness at the nose spindle is detected through the application of static forces. These spindle features are evaluated at with different supply Moreover, the describes numerical model of which consists of a 4 DOF rigid and unsymmetrical rotor supported by two aerostatic bearings. The equations of motion of the rigid rotor are solved simultaneously with the isothermal time-dependent Reynolds equation through Euler’s explicit method to obtain the lubricant pressure distribution. The main goal of the work is to validate the developed non-linear numerical model through the proposed identiﬁcation procedure and the performed experimental tests.


Introduction
Gas bearings, owing to their properties, play an essential role in turbochargers, gyroscopes, turbo blowers, gas turbines and micromachining spindles. In particular, micromachining spindles are widely used in surface finishing, printed circuit board (PCB) drilling, micro-processing, micro-milling, micro-grinding and, in general, for micromachining materials with low shear resistance. This is because compared to oil and rolling bearings, gas bearings can reach higher rotation speeds while simultaneously reducing the amount of machine maintenance required.
The aim of research on spindles supported by gas bearings is to develop systems with excellent dynamic stability and stiffness, along with very low runout of the tool, even at very high rotation speeds. Poor damping is one of the main disadvantages of journal gas bearings, which can compromise the accuracy of high-precision machining. Moreover, friction power losses in motors and bearings can be a non-negligible source of heat at high speeds. The contribution of thermally induced deformations of a machine tool can exceed 50% of the total machining error [1], which is very significant in precision milling. To reduce these problems, an efficient and precise refrigeration system is required. Recent developments in manufacturing technologies have made it possible to realize innovative high-precision devices that can even be miniaturized. Research has therefore been enriched by numerous experimental contributions, some of which are discussed below. Machine tools of 4 mm and 1/8" diameters used for micro-milling can now be supported on aerostatic bearings up to 400 krpm and 450 krpm [2,3]. Experiments on This paper proposes a cheap and simple experimental identification procedure for aerostatic spindles. Impulse and static tests are performed to identify the natural frequencies, damping ratio and static stiffness evaluated at the nose of an electro-spindle for PCB drilling. The natural frequencies and damping factors are measured by considering the impulse-response function of the system, whereas the stiffness at the nose spindle is detected through the application of static forces. These spindle features are evaluated at zero rotational speed and with different supply pressures. Moreover, the paper describes the non-linear numerical model of the spindle, which consists of a 4 DOF rigid and unsymmetrical rotor supported by two aerostatic bearings. The equations of motion of the rigid rotor are solved simultaneously with the isothermal time-dependent Reynolds equation through Euler's explicit method to obtain the lubricant pressure distribution. The main goal of the work is to validate the developed non-linear numerical model through the proposed identification procedure and the performed experimental tests. Figure 1a,b show the investigated electro-spindle. During tests, the system was kept in horizontal position by means of a nylon block and a clamp. The spindle is driven by a permanent magnet (PM) synchronous motor that is placed between the rear and front aerostatic journal bearings. A sketch of the rotor-bearings system is shown in Figure 2. The bearings have slightly different diameters: 25 mm (front bearing) and 22 mm (rear bearing). Their axial length is 30 mm. They are supplied by means of two series of supply orifices of diameter d s = 0.15 mm, located at 7 mm from the borders of the bearings. Each series includes 12 holes equally spaced along the circumferential direction. Downstream each supply hole, a small pocket of dia. 1.1 mm and depth 0.2 mm is present. Considering the manufacturing tolerances, the radial air-film thicknesses are h f = 17.5 ± 1.5 µm and h r = 21.5 ± 1.5 µm on the front and rear bearings, respectively. Table 1 lists the mass properties of the rotor. The spindle is equipped with a water-cooling system to ensure temperature stability during operation. textile spindle [27] and a mesoscopic spindle of 10 mm diameter [28]. However, the experimental campaigns to validate the models are not always time-and cost-effective. This paper proposes a cheap and simple experimental identification procedure for aerostatic spindles. Impulse and static tests are performed to identify the natural frequencies, damping ratio and static stiffness evaluated at the nose of an electro-spindle for PCB drilling. The natural frequencies and damping factors are measured by considering the impulse-response function of the system, whereas the stiffness at the nose spindle is detected through the application of static forces. These spindle features are evaluated at zero rotational speed and with different supply pressures. Moreover, the paper describes the non-linear numerical model of the spindle, which consists of a 4 DOF rigid and unsymmetrical rotor supported by two aerostatic bearings. The equations of motion of the rigid rotor are solved simultaneously with the isothermal time-dependent Reynolds equation through Euler's explicit method to obtain the lubricant pressure distribution. The main goal of the work is to validate the developed non-linear numerical model through the proposed identification procedure and the performed experimental tests. Figure 1a,b show the investigated electro-spindle. During tests, the system was kept in horizontal position by means of a nylon block and a clamp. The spindle is driven by a permanent magnet (PM) synchronous motor that is placed between the rear and front aerostatic journal bearings. A sketch of the rotor-bearings system is shown in Figure 2. The bearings have slightly different diameters: 25 mm (front bearing) and 22 mm (rear bearing). Their axial length is 30 mm. They are supplied by means of two series of supply orifices of diameter = 0.15 mm, located at 7 mm from the borders of the bearings. Each series includes 12 holes equally spaced along the circumferential direction. Downstream each supply hole, a small pocket of dia. 1.1 mm and depth 0.2 mm is present. Considering the manufacturing tolerances, the radial air-film thicknesses are ℎ = 17.5 ± 1.5 μm and ℎ = 21.5 ± 1.5 μm on the front and rear bearings, respectively. Table 1 lists the mass properties of the rotor. The spindle is equipped with a water-cooling system to ensure temperature stability during operation.    Figure 2 shows the cartesian reference system, , used to define of mass and its position with respect to the bearings. The origin, , is loc edge of the front journal bearing. Axis is along the axial direction and d front to the rear bearing. The axial coordinates of journal centers are in and for front and rear bearings, respectively.

Shaft-Displacement Detection
The rotor displacement is detected by means of two couples of cap located at 1 = −21 mm (front plane) and 2 = 91.5 mm (rear plane). The c the shaft is located at = 41 mm. The shaft center-of-mass displacemen tions about x and y axes are obtained from the sensor readings on front ( planes ( 2 , 2 ), as below:  Figure 2 shows the cartesian reference system, Oxyz, used to define the rotor center of mass and its position with respect to the bearings. The origin, O, is located at the front edge of the front journal bearing. Axis z is along the axial direction and directed from the front to the rear bearing. The axial coordinates of journal centers are indicated with z f and z r for front and rear bearings, respectively.

Shaft-Displacement Detection
The rotor displacement is detected by means of two couples of capacitance probes located at z 1 = −21 mm (front plane) and z 2 = 91.5 mm (rear plane). The center of mass of the shaft is located at z G = 41 mm. The shaft center-of-mass displacements and the rotations about x and y axes are obtained from the sensor readings on front (x 1 , y 1 ) and rear planes (x 2 , y 2 ), as below: The shaft displacement at the spindle nose (z = z nose = −40 mm) can be extrapolated using Equation (2):

Static Stiffness on the Nose
The static stiffness on the nose was evaluated by extrapolating the displacement produced through the application of a 30 ± 1 N load at z = z nose . A dynamometer was employed to impose the load, and the spindle displacement at z = z nose was computed by using Equation (2). Table 2 lists the resulting stiffness, measured at different supply pressures.

Experimental Characterization at Null Rotational Speed
The damped natural frequencies, ω d , and the damping factors, ζ, were experimentally investigated by means of impulse response tests. A vertical impulse was applied to the nose in correspondence of the tool artifact, and the spindle trajectory was detected by the capacitive sensors (capaNCDT CS05 by Microepsilon) (see Figure 1a). The tests were performed at null rotational speed and absolute supply pressures p s = 0.3, 0.5, 0.7 and 0.9 MPa. Each test was repeated four times, and signals were sampled with sampling frequency f s = 50 kHz.
The spectra of signals y 1 and y 2 , measured by sensors on the front and rear measuring planes, were evaluated through fast Fourier transform (FFT) in order to visualize the number of modes captured and their frequencies. Figure 3 shows one of the computed spectra (at 0.5 MPa supply pressure), while the related time signals are depicted in Figure 4. In all the investigated cases, only one rigid mode shape is visible. Meanwhile, the higher frequency can be attributed to the flexural mode of the spindle since it does not depend on the supply pressure of the air bearings. As can be seen from Figure 4, the first mode shape is conical since y 1 and y 2 present a 180 • phase shift. 0.7 4.0 0.9 4.6

Experimental Characterization at Null Rotational Speed
The damped natural frequencies, , and the damping factors, , were experimentally investigated by means of impulse response tests. A vertical impulse was applied to the nose in correspondence of the tool artifact, and the spindle trajectory was detected by the capacitive sensors (capaNCDT CS05 by Microepsilon) (see Figure 1a). The tests were performed at null rotational speed and absolute supply pressures = 0.3, 0.5, 0.7 and 0.9 MPa. Each test was repeated four times, and signals were sampled with sampling frequency = 50 kHz. The spectra of signals 1 and 2 , measured by sensors on the front and rear measuring planes, were evaluated through fast Fourier transform (FFT) in order to visualize the number of modes captured and their frequencies. Figure 3 shows one of the computed spectra (at 0.5 MPa supply pressure), while the related time signals are depicted in Figure  4. In all the investigated cases, only one rigid mode shape is visible. Meanwhile, the higher frequency can be attributed to the flexural mode of the spindle since it does not depend on the supply pressure of the air bearings. As can be seen from Figure 4, the first mode shape is conical since 1 and 2 present a 180° phase shift.   In view of these considerations, as a first approach, two single-DOF methods were employed to identify the modal parameters of the system: the logarithmic decrement  In view of these considerations, as a first approach, two single-DOF methods were employed to identify the modal parameters of the system: the logarithmic decrement method (LDM) [29] and the half-power method (HPM) [30]. Table 3 lists the mean values of the estimations of the damped frequencies, ω d , and the damping ratios, ζ, computed through HPM and LDM. The values of the damped natural frequencies increase with the supply pressure, whereas the damping ratios exhibit an opposite trend. The accuracy of these estimations was verified by comparing the experimental signals with those reconstructed through the identified modal parameters (see Figure 4) by making use of the following 1 DOF analytical formula: where constants a and b depend on the initial conditions, and the undamped natural Unfortunately, as can be seen from Figure 4, the reconstructed signals do not exhibit a satisfactory degree of accuracy if the estimated parameters are not manually corrected with a trial-and-error procedure. This mismatch was mainly due to the fact that in some cases, the analyzed modes presented a relatively high damping factor (small number of oscillations in the time signal) and a flat peak in the spectrum. In fact, despite the large amount of energy provided by hammering the spindle nose, in most cases, the rotor stopped after a few oscillations. To overcome these problems, the estimations were repeated by means of the least-squares complex exponential method (LSCEM).

Multi-DOF Modal Analysis with LSCEM
This is a multi-DOF method that works in the time domain [31]. It is based on the impulse-response function (IRF) of the MDOF system with N couples of complex and conjugate poles: where s r = ζ r ω n,r ± i ω n,r 1 − ζ 2 r and A r,ij are the r th poles and the modal constants of the system and i and j refer to the nodes where the impulse is applied and where the system response is measured, respectively. The advantage of this method is that the nonlinear solution for the eigenvalues is computationally very straightforward since it performs an exponential fitting that requires no starting values. To obtain a solution, it is only necessary to supply the algorithm with the impulse-response function, h ij (t), and the order of the model (N). A particular solution strategy was adopted in these cases. In order to provide the algorithm with the possibility to fit-noisy signals, the order of the model was chosen equal to 20, and frequency stabilization diagrams were used to distinguish the "computational modes" from the real ones. This kind of diagram was obtained by iterating the LSCEM from a minimum (3) to a maximum order (20) and superimposing the obtained frequencies on the obtained spectrum.
In order to obtain a cleaner stabilization diagram, before plotting the obtained damped frequencies, ω d , they were filtered by considering only the first 4 modes presenting the higher modal constant (ω 1,N , ω 2,N , ω 3,N and ω 4,N ). Figures 5 and 6 show the spectrum and time signal reconstructed by using the modal parameters obtained from this procedure. In Figure 5, the results correspond to the "+" symbols, and they are coloured in blue, red, cyan and black from the higher to the lower modal constant. A further filter was imposed by considering a threshold of 5% on the values of the normalized modal constants in an attempt to exclude the computational modes. The results that satisfied this threshold are circled in green (see Figure 5). Figure 5 also shows the comparison between the experimental spectrum and that reconstructed through the LSCEM. Figure 6 shows the same comparison in the time domain.
to provide the algorithm with the possibility to fit-noisy signals, the order of the m was chosen equal to 20, and frequency stabilization diagrams were used to distin the "computational modes" from the real ones. This kind of diagram was obtain iterating the LSCEM from a minimum (3) to a maximum order (20) and superimp the obtained frequencies on the obtained spectrum.
In order to obtain a cleaner stabilization diagram, before plotting the obt damped frequencies, , they were filtered by considering only the first 4 modes pr ing the higher modal constant ( 1, , 2, , 3, and 4, ). Figures 5 and 6 show the trum and time signal reconstructed by using the modal parameters obtained from procedure. In Figure 5, the results correspond to the "+" symbols, and they are col in blue, red, cyan and black from the higher to the lower modal constant. A further was imposed by considering a threshold of 5% on the values of the normalized m constants in an attempt to exclude the computational modes. The results that satisfie threshold are circled in green (see Figure 5). Figure 5 also shows the comparison bet the experimental spectrum and that reconstructed through the LSCEM. Figure 6 s the same comparison in the time domain.  As can be seen, in this case, the accuracy of the method is satisfactory. The r obtained with the LSCEM method are summed up in Table 4, where the damped na frequencies and damping ratios are reported, together with their phase, , and m As can be seen, in this case, the accuracy of the method is satisfactory. The results obtained with the LSCEM method are summed up in Table 4, where the damped natural frequencies and damping ratios are reported, together with their phase, φ, and modal constant, A. Only the first rigid mode was considered, as the second mode exhibited a natural frequency of about 220 krpm, which was almost constant with the supply pressure and had modal constants much smaller than the first mode. The evaluation of the phase shifts, ∆φ = |φ(y 1 ) − φ(y 2 )|, confirms that the mode identified is of the tilting type. Table 4. Experimental natural frequency and damping ratio of signals y 1 and y 2 at different supply pressures (only the first rigid mode is reported); the results were obtained with LSCEM.
Signal The damped frequencies obtained with the 1 DOF methods are in good accordance with those obtained with LSCEM; the error is less than 4%. The damping factors, apart from cases with p s = 0.3 MPa, differ by less than 36% for HPM and 61% for LDM.

Numerical Model
In this section, the numerical model of the spindle described in Section 2.1 is presented. The effect of the double-face thrust bearing is not taken into account, as it is negligible in the determination of radial and tilting stiffness with respect to the contribution of journal bearings.
The time-dependent Reynolds equation is considered in the fluid domain of journal bearings: where g in is the input flow per unit surface that crosses the supply orifices. The input flow is estimated using the analytical expression of the supply hole's discharge coefficient [32] without considering the Reynolds number dependence: where d s is the supply hole's diameter and h is the local air gap under the orifice. The RE is discretized with finite difference technique, and the pressure distribution at time t + 1 is obtained from the previous pressure distribution with the Euler explicit method: where f is a non-linear function that involves the pressure values in node i,j and in the adjacent nodes. The pocket downstream of each supply orifice is also considered, as it influences the dynamic characteristics of bearings. Once pressure p ij is calculated at time t + 1, this value is immediately used to calculate p i+1,j at the same time. This improves the stability of the numerical scheme, as it facilitates the transmission of the information in the adjacent nodes. The film thickness depends on the rotor center-of-mass position (x G , y G ) and on the rotor axis orientation (ϑ x , ϑ y ): Once the pressure distribution is calculated, the shear stress on the rotor surface is given by: The 4 DOF equations of motion for a rigid rotor are considered: where e and γ are the static and the dynamic rotor unbalances, respectivelye; ϕ 1 is the angle that locates the dynamic unbalance plane with respect to the static unbalance plane; F ext,x and F ext,y are the external-force components acting on the rotor at coordinates z nose ; F cx , F cy , M cx and M cy are the reaction forces and moments of the bearings, calculated by integrating the pressure distribution and the shear stress in the fluid domain.
The moments are expressed with respect to the rotor center of mass: The moments are expressed with respect to the rotor center of mass: The time integration of these equations of motion are carried out with the first order Euler method: where q is the state-space vector defined by The algorithm performed for the coupled integration of the ODE and the PDE system is described below. After assuming an initial pressure distribution, for each time iteration, the following steps are followed: 1.
The reaction forces and moments are calculated using Equations (11) and (12); 2.
The center-of-mass accelerations and the angular accelerations of the rotor are computed from the equations of motion (10); 3.
The state-space vector at time t + 1 is obtained from Equation (13); 4.
The film thickness at time t + 1 is updated from Equation (8); 5.
The pressure and shear-stress distributions at time t + 1 are updated by solving RE (2) and using Equation (9); 6.
Go back to point 1.
This method is known as the "orbit method" and has the advantage of also taking into account the non-linearities of bearings [22]. Conversely, it is more time-consuming compared to linearized methods.
For static problems, the rotor position is fixed, and only the pressure distribution is considered time-dependent. In this case, only steps 1 and 5 are considered in order to find the steady-state solution. The iterative process is stopped when the load-capacity relative variation decreases below a given tolerance:

Numerical Results
The numerical model developed was validated using both literature benchmarks and numerical results. In particular, the numerical load capacity, W, and attitude angle, Φ, were compared with the analytical values for plain dynamic journal bearings. In the comparison, the static stiffness of the spindle measured in correspondence of the nose was also considered, together with the damped natural frequencies and the damping ratio at null rotational speed.

Static Validation with Analytical Solution of Plain Dynamic Journal Bearings
The mathematical model was validated in static conditions, comparing it with the analytical solution that exists for the plain dynamic journal bearing [7]. The dimensionless load capacity is expressed with a real and imaginary parts, which represent the in-phase and the out-of-phase components with respect to the line of centers: where is the bearing number, and the eccentricity ratio. The front and rear journal bearings were simulated with no supply orifices for comparison with the dynamic-bearing analytical solution. Figure 7 shows a perfect match, both regarding the load capacity and the attitude angle, Φ.

Choice of Grid Resolution
The influence of the grid-point number on the numerical solution is evaluated both in static and dynamic conditions. A total of 31 nodes were chosen along the axial direction for each bushing. Along the circumferential direction, 48 or 96 nodes were considered. Figures 8-10 compare the pressure distribution, load capacity and air consumption obtained with these computational grids in the presence of an eccentricity of = 1 μm. As can be seen, the difference is minimal. As a result of these considerations, a 31 × 48 grid was used.

Choice of Grid Resolution
The influence of the grid-point number on the numerical solution is evaluated both in static and dynamic conditions. A total of 31 nodes were chosen along the axial direction for each bushing. Along the circumferential direction, 48 or 96 nodes were considered. Figures 8-10 compare the pressure distribution, load capacity and air consumption obtained with these computational grids in the presence of an eccentricity of e = 1 µm. As can be seen, the difference is minimal. As a result of these considerations, a 31 × 48 grid was used.
The influence of the grid-point number on the numerical solution is evaluated both in static and dynamic conditions. A total of 31 nodes were chosen along the axial direction for each bushing. Along the circumferential direction, 48 or 96 nodes were considered. Figures 8-10 compare the pressure distribution, load capacity and air consumption obtained with these computational grids in the presence of an eccentricity of = 1 μm. As can be seen, the difference is minimal. As a result of these considerations, a 31 × 48 grid was used.

Static Stiffness of Journal Bearings
The stiffness of journal bearings at zero speed was obtained from static simulations. A given displacement was imposed on the shaft (1 μm on radius), and the reaction forces in the bearings were computed in order to estimate the radial stiffness, and , of the front and rear bearings, respectively. The stiffness values are listed in Table 5, considering minimal and maximal radial clearances and taking into account the manufacturing tolerances, as shown in Section 2.1. Of course, stiffness values increase with supply pressure and decrease with increasing air gaps. The front bearing shows an increment of stiffness with the air gap only with = 0.3 MPa. The following explanation can be given after

Static Stiffness of Journal Bearings
The stiffness of journal bearings at zero speed was obtained from static simulations. A given displacement was imposed on the shaft (1 µm on radius), and the reaction forces in the bearings were computed in order to estimate the radial stiffness, k f and k r , of the front and rear bearings, respectively. The stiffness values are listed in Table 5, considering minimal and maximal radial clearances and taking into account the manufacturing tolerances, as shown in Section 2.1. Of course, stiffness values increase with supply pressure and decrease with increasing air gaps. The front bearing shows an increment of stiffness with the air gap only with p s = 0.3 MPa. The following explanation can be given after investigating the pressure distribution: the reason is that the front bearing is near the saturation condition; that is to say that the pressure drop through the input orifices is quite small (less than 0.03 MPa). The tilting stiffness of each journal bearing, calculated with respect to its center, is also evaluated in order to investigate whether its contribution to the overall tilting stiffness of the rotor bearings system is negligible or not. Table 6 compares the contribution of each bearing to the tilting stiffness obtained from the radial stiffness (left column) with the tilting stiffness of each bearing with respect to its center (right column). From this comparison, it is evident that in this case, this last bearing can be neglected so the overall tilting stiffness can be estimated from the radial stiffness of each journal bearing.

Stiffness on the Nose
The overall stiffness, k nose , measured on the nose at z = z F = −40 mm, is a function of the bearing stiffness according to: Table 7 lists the radial stiffness on the nose extrapolated from the journal bearings stiffness calculated in the centered position, with 1 µm rotor eccentricity. In case an external radial load of 30 N is applied on the nose, the eccentricities on bearings are greater and the resulting nose stiffness is lower. Table 8 lists these values.

Simplified Natural Frequencies
At null rotational speed, in case the damping effects are neglected, the cylindrical natural frequency is expressed by: The conical natural frequency is given by: where l f and l r are the distances between the rotor center of mass and the centers of the front and the rear journal bearings, respectively. In this formulation, the eventual tilting stiffness of each journal bearing is neglected. Table 9 reports the undamped natural frequencies of the rigid rotor-bearings system at null rotational speed.

Choice of the Time Step
The optimal time step, dt, must be chosen for simulations in order to guarantee a stable solution and independence of the dynamic solution on the time step. The entity of the damping factor was found to be dependent on the time step, dt, and on the grid refinement. Figure 11 compares the values obtained for the damping factor at different rotational speeds and time steps, dt, with a 31 × 96 grid. As shown, the convergence is reached with a time step of dt = 10 −8 s, and a further decrease does not significantly change the damping factor. The effect of the grid on the damping factor is visible in Figure 12. When doubling the grid nodes along the circumferential direction, it is sufficient to halve the time step, dt, in order to obtain a very similar solution. As a result of these considerations, a 31 × 48 grid with dt = 2 × 10 −8 s was chosen for dynamic simulations, as it was found to be a good compromise between accuracy and computational time.

Damped Natural Frequencies and Damping-Factor Identification
The damped natural frequencies and the damping factors were evaluated for the displacements along the x and y directions by the rotor in correspondence of planes at z = 0 (at the right edge of the front journal bearing) and at z = z L = 126 mm (at the left edge of the rear journal bearing). Moreover, to better identify both the cylindrical and the conical modes, the center-of-mass displacement, y G , and the tilting angles, ϑ x and ϑ y , were considered. Two different initial conditions were simulated since the frequency content of the spectra of the system response may depend on how the system is excited. Condition 1 is defined by: Condition 1 was adopted in an attempt to reproduce the performed experimental tests. Vertical impulses were applied to the spindle nose, and the rotor oscillations were measured on the rear and front plane. Since the results obtained with the initial conditions (condition 1) were mainly governed by the conical mode, a second condition (condition 2) was simulated in an attempt to better figure out the cylindrical mode-shape of the rotor.
ϑ y (t = 0) = 0 Condition 2 was chosen to simultaneously produce a whirl and a translation of the rotor. Due to the non-symmetry of the rotor-bearings system, both modes were visible both with initial condition 1 and with condition 2. The frequency content of rotor displacements puts in evidence the general presence of two natural frequencies, as seen in Figure 13. The prominence of one with respect to the other depends on the initial condition.
Condition 2 was chosen to simultaneously produce a whirl and a transl rotor. Due to the non-symmetry of the rotor-bearings system, both modes w both with initial condition 1 and with condition 2. The frequency content of rot ments puts in evidence the general presence of two natural frequencies, as see 13. The prominence of one with respect to the other depends on the initial con Figure 13. Examples of frequency spectra of signals , 0 , and obtained fr method.

Logarithmic Decrement Method (LDM)
Although the LSCEM proved to be the most reliable method, for the sa pleteness, the LDM method was also used to identify the main damping facto largest modal constant) and natural frequencies of the numerical signals 0 , . The LDM was used, considering time intervals Δ of 5 and 10 ms. The resu of the damping factors were compared to evaluate the sensitivity of the meth spect to the signal length. The main frequency of each signal was also evaluate uring the average period of the oscillations present in the considered time int frequency can be compared with the peaks visible in the related spectrum. Ta in Appendix A list both the natural frequencies resulting from the FFT analy frequencies and the damping factor of the signal resulting from the logarithmic method. When the oscillation was highly damped ( = 0.3 and 0.5 MPa), only window was considered. The following considerations can be made based o analysis: Figure 13. Examples of frequency spectra of signals y G , y 0 , θ x and y L obtained from the orbit method.

Logarithmic Decrement Method (LDM)
Although the LSCEM proved to be the most reliable method, for the sake of completeness, the LDM method was also used to identify the main damping factors (with the largest modal constant) and natural frequencies of the numerical signals y 0 , y L , y G and θ x . The LDM was used, considering time intervals ∆T of 5 and 10 ms. The resulting values of the damping factors were compared to evaluate the sensitivity of the method with respect to the signal length. The main frequency of each signal was also evaluated by measuring the average period of the oscillations present in the considered time intervals. This frequency can be compared with the peaks visible in the related spectrum. Tables A1-A4 in Appendix A list both the natural frequencies resulting from the FFT analysis and the frequencies and the damping factor of the signal resulting from the logarithmic decrement method. When the oscillation was highly damped (p s = 0.3 and 0.5 MPa), only a 5 ms time window was considered. The following considerations can be made based on the data analysis:

•
The frequency obtained from the average period of the oscillation is quite similar to one of the two frequencies detected in the FFT spectrum; in particular, the higher frequency is visible only on the rear plane (signal y L ), while the lower frequency is present on the other signals (y 0 , y G and ϑ x ); • in most cases, the application of LDM to signals with durations of 5 and 10 ms does not influence the results; • depending on the initial condition, one of the two modes prevails; • due to the short duration of the vibration at p s = 0.3 MPa, which is highly damped, the FFT algorithm does not provide a good estimation of the frequency spectra; • the frequency increases with the supply pressure, while the damping factor decreases; • by increasing the air gap, the frequency decreases, as well as the damping factor.

LSCE Method
The frequency content of the numerical signals y L and y 0 was investigated by adopting a procedure similar to that adopted in Section 2.4. To obtain more reliable frequencystabilization diagrams, the adopted procedure was slightly different than that used for the experimental data. In this case, to favor the exponential fitting, the model order used in the LSCEM was 80, and a random noisy signal was added to the numerical signals [31]. Figures 14 and 15 show the comparison obtained between the original and reconstructed data. The damped frequencies and damping factors evaluated on signals y 0 and y L are listed in Table 10, together with their phase shift and modal constants.
 the frequency increases with the supply pressure, while the damping fa creases;  by increasing the air gap, the frequency decreases, as well as the damping fa

LSCE Method
The frequency content of the numerical signals and 0 was investig adopting a procedure similar to that adopted in Section 2.4. To obtain more reli quency-stabilization diagrams, the adopted procedure was slightly different th used for the experimental data. In this case, to favor the exponential fitting, th order used in the LSCEM was 80, and a random noisy signal was added to the nu signals [31]. Figures 14 and 15 show the comparison obtained between the orig reconstructed data. The damped frequencies and damping factors evaluated on 0 and are listed in Table 10, together with their phase shift and modal consta

Comparison of LSCEM and LDM for Numerical Simulations
When LDM and LSCEM are compared for the numerical simulations, the damped frequencies obtained with LSCEM differ by less than 1.6% with respect the frequencies obtained with LDM, apart from the case with p s = 0.3 MPa. The estimations of damping factors differ by much more-up to 80%. This parameter is more sensitive to variations in the air gap and pressure.

Comparison between Experimental and Numerical Results
The first critical speed of the spindle was experimentally measured by analyzing the amplitude of the synchronous whirl at different rotational speeds. It was found to be 72 krpm and 94 krpm at gauge supply pressure p s = 0.5 and 0.7 MPa, respectively. If compared with the theoretical frequencies resulting from the simplified model of Section 4.5, the difference is less than 6%.
The experimental stiffness on the nose at different pressures is lower than the theoretical value calculated near the rotor-centered condition (between 27% and 39% with respect to the theoretical value). Considering that the 30 N external load imposed on the nose involves non-negligible eccentricity ratios on bearings (especially on the front bearing and with low pressure), a non-linear estimation of the stiffness is preferred (see Section 4.4). If compared with the stiffness obtained in this case, the experimental stiffness is more in accordance with the prediction (between 52% and 86% with respect to the theoretical value).
The experimental damped frequencies and the damping factors evaluated on the shaft vertical vibration at z = z L are compared with the numerical frequencies in Figure 16. The experimental damped natural frequencies are in good accordance with the numerical frequencies, as they fall in the range defined on the basis of the tolerance considered on the front (h f = 17.5 ± 1.5 µm) and rear (h r = 21.5 ± 1.5 µm) clearances. Conversely, the computed experimental damping factors are higher than the numerical simulations. The single-DOF methods were found to be less accurate than LSCEM to capture the modal parameters of experimental signals, especially the damping factors.
ical frequencies, as they fall in the range defined on the basis of the tolerance considered on the front (ℎ = 17.5 ± 1.5 μm) and rear (ℎ = 21.5 ± 1.5 μm) clearances. Conversely, the computed experimental damping factors are higher than the numerical simulations. The single-DOF methods were found to be less accurate than LSCEM to capture the modal parameters of experimental signals, especially the damping factors.

Conclusions
This paper describes the experimental and numerical investigation of an electrospindle supported by aerostatic bearings. The numerical model was validated against literature results in the case of plain aerodynamic bearings and against static and dynamic experimental tests carried out on the electro-spindle at null rotational speed. Both single-DOF methods and a multi-DOF method were proposed to estimate the damped natural frequencies and the damping factors. The accuracy of the model results with respect to the experimental data was discussed. In particular, the damping factor evaluated with LSCEM was in good agreement with the experimental data, while the single-DOF methods overestimated it. Future investigations will regard the dynamic behavior of the spindle at different rotational speeds in order to study its unbalance response and stability.
Author Contributions: Conceptualization, F.C., L.L., A.T., T.R. and V.V.; methodology, F.C. and L.L.; software, F.C. and L.L.; writing-review and editing, F.C., L.L. and A.T.; supervision, V.V.; project administration, T.R.; funding acquisition, T.R. All authors have read and agreed to the published version of the manuscript. Radial stiffness of the spindle evaluated on the nose k r Radial stiffness of the rear journal bearing k ϑ f Tilting stiffness of front journal bearing with respect to its center k ϑ r Tilting stiffness of rear journal bearing with respect to its center L Journal bearing length l f , l r Axial distance between the rotor center of mass and the centers of the front and rear journal bearings M c Reaction moment in bearings m r Mass of rotor N Order of the least squared complex exponential (LSCE) fitting p Absolute pressure p a Ambient pressure p s Supply absolute pressure r Journal bearing radius s r r th pole of the system x Generic position measured along the x-axis x G Center-of-mass position measured along the x-axis x nose Spindle nose position measured along the x-axis x 1 Position measured along the x-axis in the front plane x 2 Position measured along the x-axis in the rear plane y Generic position measured along the y-axis y G Center-of-mass position measured along the y-axis y nose Spindle nose position measured along the y-axis y rec Signal reconstructed by means of the identified modal parameters y 1 Position measured along the y-axis in the front plane y 2 Position measured along the y-axis in the rear plane z Generic position measured along the z-axis z f Front-bearing center axial coordinate z G Rotor center-of-mass axial coordinate z L Axial coordinate of the rear end of the rotor z nose Spindle nose position measured along the z-axis z r Front-bearing center axial coordinate z 0 Axial coordinate of the front end of the rotor z 1 Axial position of the front measuring plane z 2 Axial position of the rear measuring plane W Dimensionless load capacity t Time ∆T Time inteval used for LDM Eccentricity ratio φ Phase angle ϕ Angle identifing dynamic unbalance of the rotor Φ Attitude angle Λ Bearing number µ Air viscosity γ Dynamic unbalance of the rotor ϑ x , ϑ y Rotations around x and y axes τ Shear stress ζ Damping factor ω Angular speed ω con Conical mode-shape frequency of the spindle ω cyl Cylindrical mode-shape frequency of the spindle ω d Damped natural frequency ω n Undamped natural frequency

Appendix A
Tables A1-A4 report both the natural frequencies resulting from the FFT analysis, both the frequencies and the damping factor of the signal resulting from the logarithmic decrement method.