Experimental and Numerical Dynamic Identification of an Aerostatic Electro-Spindle

Featured Application

Aerostatic spindles are used in high-speed micromachining applications. 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.

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 non-linear 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 identification procedure and the performed experimental tests. The comparison proves satisfactory, and the possible sources of uncertainty are conjectured.

1. 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 dynamic characterization of air-bearing spindles are described in [4,5], where the frequency-response function (FRF) of the spindle is measured. Impact tests have been carried out with an impact hammer [4] or with a custom-designed impact excitation system able to provide excitations up to 20 kHz [5]. Ref. [6] describes a test rig used to measure the dynamic stiffness and damping coefficients of a hemisphere spiral groove hybrid gas bearing. In this rig, the excitation is provided by means of two electromagnetic exciters.

Mathematical models have long been used as a valuable tool to investigate the dynamic performance and stability of high-speed rotors. Early attempts at modelling this kind of system were aimed at gaining a better insight into the half-whirl instability, known today as asynchronous whirl instability [7]. At the time, this was one of the most troublesome problems in journal bearings. The first attempts to actually solve the equations of motion and the gas-lubricated equation simultaneously were reported by Sternlicht [8], Poritsky and Arwas [9]. In their solution, they linearized the equations of motion of a 2 DOF rotor, along with the Reynolds equation, to obtain lubricant pressure distribution. However, in the Reynolds equation, they omitted the time derivative of pressure, and therefore, the resulting quasi-static solutions do not adequately represent the squeeze-film effects. In the same period, the literature presented many other models based on the “$p$-$h$” linearization of the time-dependent Reynolds equation and the equations of motion of a 4 DOF rigid rotor [10,11]. However, in most cases, these kinds of simplified models cannot replace more costly trial-and-error experimentation since they do not represent a reliable theory.

The results obtained from these models have become ever more accurate with the use of finite-difference and finite-element models, where the perturbation method can be used to compute stiffness and damping coefficients for rotordynamic models. Ref. [12] investigates the impulse response of an ultra-precision raster milling spindle under the simplified hypothesis of constant coefficients. In [13], critical speeds and modal shapes of a micro-spindle machine tools were evaluated, but this analysis is limited because it only considers bearings with constant stiffnesses and neglects their damping properties. A rotordynamic model of an air spindle with a diameter of 1/8″ was considered [14], which involves coefficients that depend on the rotational speed. Here, the critical rotational speeds were calculated as a result of the rotordynamic analysis. A test rig to measure the rotordynamic response of a spindle supported on porous gas bearings is described in [15]. In [16,17,18,19], the dependence of stiffness and damping properties on both rotational speed and perturbation frequency is expressed with analytical formulations. With this technique, it is possible to calculate the critical speeds and verify the stability of rigid or flexible rotors. During the design process of a spindle supported on gas bearings, it is of great importance to estimate the conical and cylindrical natural frequencies. A clear design guideline is illustrated in [20] to avoid critical speeds within the operating speed range. The influence of temperature on the dynamic behavior of the spindle is also discussed. Conical and cylindrical modes are investigated in [21] for a non-uniform slot-clearance journal bearing, together with the onset speed of whirl instability.

The orbit method is an alternative to the linearized coefficient analysis employed in the aforementioned papers. It includes the complete nonlinear equations of the rotor and the bearings, which are integrated numerically to obtain the shaft center orbits corresponding to any set of geometrical, running, and initial conditions. This technique essentially uses the computer as an accurate experimental rig; it operates exactly in accordance with the assumed governing equations [22]. This method is currently used to calculate the unbalance response, the onset speed and the non-linear rotor behavior in case of large displacements [23]. Ref. [24] makes use of this method for a spindle with herringbone-grooved journal bearings. Ref. [25] highlights the differences between linear and non-linear analyses of a rotor supported by plain journal bearing. Past research by the authors concerns an electro-spindle of 50 mm shaft dia. supported on aerostatic bearings [26], a 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.

2. Experimental Setup

2.1. The Prototype

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. Mass properties of rotor.

Rotor mass:	$m_r$ = 382 g
Transversal inertia moment:	$I$ = 666.1 × 10−6 kg m2
Polar inertia moment:	$I_p$ = 31.1 × 10−6 kg m2

lists the mass properties of the rotor. The spindle is equipped with a water-cooling system to ensure temperature stability during operation.

2.2. Shaft-Displacement Detection

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.

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:

$$\begin{gathered} x_G=\frac{x_1\left(z_2-z_G\right) - x_2\left(z_1-z_G\right)}{z_2-z_1} \\ {y}_G=\frac{y_1\left(z_2-z_G\right) - y_2\left(z_1-z_G\right)}{z_2-z_1} \\ {\vartheta }_y=\frac{x_2-x_1}{z_2-z_1} \\ {\vartheta }_x=\frac{y_1-y_2}{z_2-z_1} \end{gathered}$$

(1)

The shaft displacement at the spindle nose ($z$ = $z_{nose}$ = −40 mm) can be extrapolated using Equation (2):

$$\begin{gathered} x_{nose}=\frac{x_1\left(z_2-z_{nose}\right) - x_2\left(z_1-z_{nose}\right)}{z_2-z_1} \\ {y}_{nose}=\frac{y_1\left(z_2-z_{nose}\right) - y_2\left(z_1-z_{nose}\right)}{z_2-z_1} \end{gathered}$$

(2)

2.3. Static Stiffness on the Nose

The static stiffness on the nose was evaluated by extrapolating the displacement produced through the application of a 30 $\pm 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. Resulting static stiffness on the nose.

${\mathit{p}}_{\mathit{s}}\left(\mathbf{MPa}\right)$	${\mathit{k}}_{\mathit{n}\mathit{o}\mathit{s}\mathit{e}}(\mathbf{N}/\mathsf{\mu }\mathbf{m})$
0.3	1.14
0.5	2.9
0.7	4.0
0.9	4.6

lists the resulting stiffness, measured at different supply pressures.

2.4. Experimental Characterization at Null Rotational Speed

The damped natural frequencies, ${\omega }_d$, and the damping factors, $\zeta$, 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.

2.4.1. Single-DOF Modal-Identification Analysis through LDM and HPM

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. Experimental natural frequency and damping ratio of signal y₂ at different supply pressures; the results were obtained with HPM and LDM.

${\mathit{p}}_{\mathit{s}}$ (MPa)	HPM		LDM
	${\mathit{\omega }}_{\mathit{d}}$ (krpm)	$\mathit{\zeta }$	${\mathit{\omega }}_{\mathit{d}}$ (krpm)	$\mathit{\zeta }$
0.3	70.23	0.2611	72.82	0.3928
0.5	87.63	0.1464	86.50	0.2347
0.7	95.42	0.1084	92.60	0.1497
0.9	100.3	0.1125	97.64	0.2325

lists the mean values of the estimations of the damped frequencies, ${\omega }_d$, and the damping ratios, $\zeta$, 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:

$$\begin{gathered} y_{rec}\left(t\right) = \left[a\cos \left({\omega }_{d}t\right) + b\sin \left({\omega }_{d}t\right)\right] · e^{-\zeta {\omega }_{n}t} \\ a=y_{rec}\left(0\right) \\ b=\frac{{\dot{y}}_{rec}\left(0\right) + \zeta {\omega }_n{y}_{rec}\left(0\right)}{{\omega }_d} \end{gathered}$$

(3)

where constants a and b depend on the initial conditions, and the undamped natural frequency is ${\omega }_n={\omega }_d/\sqrt{1-{\zeta }^2}$.

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).

2.4.2. 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:

$$h_{ij}\left(t\right) ={\displaystyle \sum }_{r=1}^{2N}{A}_{r,ij} · e^{s_{r}t}$$

(4)

where $s_r={\zeta }_r{\omega }_{n,r}\pm i {\omega }_{n,r}\sqrt{1-{\zeta }_r^2}$ 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}\left(t\right)$, 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, ${\omega }_d$, they were filtered by considering only the first 4 modes presenting the higher modal constant (${\omega }_{1,N}$, ${\omega }_{2,N}$, ${\omega }_{3,N}$ and ${\omega }_{4,N}$). Figure 5 and Figure 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.

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. Experimental natural frequency and damping ratio of signals y₁ and y₂ at different supply pressures (only the first rigid mode is reported); the results were obtained with LSCEM.

${\mathit{p}}_{\mathit{s}}$ (MPa)	Signal	${\mathit{\omega }}_{\mathit{d}}$(krpm)	$\mathit{\zeta } (-)$	$\mathit{\varphi } \left(\mathbf{deg}\right)$	${\mathit{A}}_1 (-)$
0.3	$y_2$	73.44	0.1483	3.474	3.400
	$y_1$	70.21	0.2578	169.8	7.963
0.5	$y_2$	89.50	0.1072	9.420	5.305
	$y_1$	87.23	0.1632	178.0	7.673
0.7	$y_2$	99.75	0.00762	−3.336	5.223
	$y_1$	95.46	0.1060	153.6	8.461
0.9	$y_2$	103.8	0.07330	19.69	3.150
	$y_1$	107.0	0.09230	−154.6	4.102

, where the damped natural frequencies and damping ratios are reported, together with their phase, $\varphi$, 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, $\mathsf{\Delta }\varphi = \left|\varphi \left(y_1\right) - \varphi \left(y_2\right)\right|$, confirms that the mode identified is of the tilting type.

2.4.3. LSCEM vs. LDM and HPM (Experimental)

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.

3. 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:

$$\frac{\partial }{\partial z}\left(\frac{p{h}^3}{12\mu RT}\frac{\partial p}{\partial z}\right) + \frac{\partial }{r\partial \vartheta }\left(\frac{p{h}^3}{12\mu RT}\frac{\partial p}{r\partial \vartheta }\right) + g_{in}=\frac{\omega r}{2RT}\frac{\partial \left(ph\right)}{r\partial \vartheta }+\frac{1}{RT}\frac{d\left(ph\right)}{dt}$$

(5)

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:

$$c_d=0.85\left(1-e^{-8.2\frac{h}{d_s}}\right)$$

(6)

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:

$$p_{ij}^{t+1}=p_{ij}^t+\mathsf{\Delta }t · f\left(p^t, p^{t+1}\right)$$

(7)

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 (${\vartheta }_x,{\vartheta }_y)$:

$$h\left(z,\vartheta \right) =h_0- \left[x_G+{\vartheta }_y\left(z-z_G\right)\right] \cos \vartheta - \left[y_G-{\vartheta }_x\left(z-z_G\right)\right]\sin \vartheta$$

(8)

Once the pressure distribution is calculated, the shear stress on the rotor surface is given by:

$$\tau =\frac{h}{2}\frac{\partial p}{r\partial \vartheta }+\mu \frac{\omega r}{h}$$

(9)

The 4 DOF equations of motion for a rigid rotor are considered:

$$\begin{gathered} m_r{\ddot{x}}_G-F_{cx}-F_{ext, x}=m_{r}e{\omega }^2\cos \omega t \\ {m}_r{\ddot{y}}_G-F_{cy}-F_{ext, y}=m_{r}e{\omega }^2\sin \omega t \\ I{\ddot{\vartheta }}_x+I_p\omega {\dot{\vartheta }}_y-M_{cx}-F_{ext, y}\left(z_G-z_{nose}\right) = \left(I_p-I\right)\gamma {\omega }^2\sin \left(\omega t+{\phi }_1\right) \\ I{\ddot{\vartheta }}_y-I_p\omega {\dot{\vartheta }}_x-M_{cy}+F_{ext, x}\left(z_G-z_{nose}\right) = \left(I-I_p\right){\omega }^2\gamma \cos \left(\omega t+{\phi }_1\right) \end{gathered}$$

(10)

where $e$ and $\gamma$ are the static and the dynamic rotor unbalances, respectivelye; ${\phi }_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:

$$\left[\begin{array}{c}{F}_{cx}\\ F_{cy}\end{array}\right] ={{\displaystyle \int }}_0^L{{\displaystyle \int }}_0^{2\pi }\left(-p\left[\begin{array}{c}\cos \vartheta \\ \sin \vartheta \end{array}\right] + \tau \left[\begin{array}{c}\sin \vartheta \\ -\cos \vartheta \end{array}\right]\right)rd\vartheta dz$$

(11)

The moments are expressed with respect to the rotor center of mass:

$$\left[\begin{array}{c}{M}_{cx}\\ M_{cy}\end{array}\right] ={{\displaystyle \int }}_0^L{{\displaystyle \int }}_0^{2\pi }\left[\begin{array}{c}\left(z-z_G\right)\left(p\sin \vartheta +\tau \cos \vartheta \right)\\ \left(z-z_G\right)\left(-p\cos \vartheta +\tau \sin \vartheta \right)\end{array}\right]rd\vartheta dz$$

(12)

The time integration of these equations of motion are carried out with the first order Euler method:

$$q_{t+1}=q_t+{\dot{q}}_t\mathsf{\Delta }t$$

(13)

where q is the state-space vector defined by

$$q= \left[\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}{x}_G& y_G\end{array}& \begin{array}{cc}{\vartheta }_x& {\vartheta }_y\end{array}\end{array}& \begin{array}{cc}\begin{array}{cc}{\dot{x}}_G& {\dot{y}}_G\end{array}& \begin{array}{cc}{\dot{\vartheta }}_x& {\dot{\vartheta }}_y\end{array}\end{array}\end{array}\right]$$

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:

• The reaction forces and moments are calculated using Equations (11) and (12);
• The center-of-mass accelerations and the angular accelerations of the rotor are computed from the equations of motion (10);
• The state-space vector at time $t+1$ is obtained from Equation (13);
• The film thickness at time $t+1$ is updated from Equation (8);
• The pressure and shear-stress distributions at time $t+1$ are updated by solving RE (2) and using Equation (9);
• 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:

$$\left|\frac{F_c^{t+1}-F_c^t}{F_c^t}\right| <{10}^{-6}$$

(14)

4. Numerical Results

The numerical model developed was validated using both literature benchmarks and numerical results. In particular, the numerical load capacity, $\overline{W}$, and attitude angle, $\mathsf{\Phi }$, 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.

4.1. 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:

$$\overline{W}=\frac{W}{ϵp_{a}LD}=\frac{i\mathsf{\Lambda }}{1+i\mathsf{\Lambda }}\frac{\pi }{2}\left[1-\frac{\tanh \left(\sqrt{1+i\mathsf{\Lambda }}\frac{L}{D}\right)}{\sqrt{1+i\mathsf{\Lambda }}\frac{L}{D}}\right]$$

(15)

where

$$\mathsf{\Lambda }=\frac{6\mu \omega }{p_a}{\left(\frac{r}{C}\right)}^2$$

(16)

is the bearing number, and

$$ϵ=\frac{e}{C}$$

(17)

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, $\mathsf{\Phi }$.

4.2. 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. Figure 8, Figure 9 and Figure 10 compare the pressure distribution, load capacity and air consumption obtained with these computational grids in the presence of an eccentricity of $e$ = 1 $\mathsf{\mu }$m. As can be seen, the difference is minimal. As a result of these considerations, a 31 × 48 grid was used.

4.3. 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. Radial stiffness of journal bearings at different supply pressures and radial clearances.

${\mathit{p}}_{\mathit{s}}$(MPa)	${\mathit{h}}_{\mathit{f}}$ = 16 μm,  ${\mathit{h}}_{\mathit{r}}$  = 20  μm		${\mathit{h}}_{\mathit{f}}$ = 19 μ$\mathbf{m},{\mathit{h}}_{\mathit{r}}=23\mathsf{\mu }\mathbf{m}$
	${\mathit{k}}_{\mathit{f}}$  (N/μm)	${\mathit{k}}_{\mathit{r}}$  (N/μm)	${\mathit{k}}_{\mathit{f}}(\mathbf{N}/\mathsf{\mu }\mathbf{m})$	${\mathit{k}}_{\mathit{r}} (\mathbf{N}/\mathsf{\mu }\mathbf{m})$
0.3	4.51	4.57	5.53	4.22
0.5	12.6	9.37	11.8	7.45
0.7	20.3	13.0	16.7	9.5
0.9	27.0	15.8	20.6	10.9

, 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. Comparison between different contributions to tilting stiffness; results refer to geometrical configuration with $h_f$ = 19 μm, $h_r$ = 23 μm and supply pressure $p_s$ = 0.7 MPa.

Tilting Stiffness of Journal Bearing with Respect to the Shaft Center of Mass		Tilting Stiffness of Journal Bearing Respect to Its Center
${\mathit{k}}_{\mathit{f}}·({\mathit{z}}_{\mathit{G}}-{\mathit{z}}_{\mathit{f}})/{\mathit{\theta }}_{\mathit{x}}(\mathbf{Nm}/\mathbf{rad})$	${\mathit{k}}_{\mathit{r}}·({\mathit{z}}_{\mathit{G}}-{\mathit{z}}_{\mathit{r}})/{\mathit{\theta }}_{\mathit{x}}(\mathbf{Nm}/\mathbf{rad})$	${\mathit{k}}_{{\mathit{\vartheta }}_{\mathit{f}}}(\mathbf{Nm}/\mathbf{rad})$	${\mathit{k}}_{{\mathit{\vartheta }}_{\mathit{r}}}(\mathbf{Nm}/\mathbf{rad})$
33507	57166	925	534

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.

4.4. 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:

$$k_{nose}=k_f\frac{z_r-z_f}{z_r-z_F}$$

(18)

Table 7. Radial stiffness on the nose at different supply pressures and radial clearances, calculated at small eccentricity (1 μm).

${\mathit{p}}_{\mathit{s}}$(MPa)	${\mathit{k}}_{\mathit{n}\mathit{o}\mathit{s}\mathit{e}}(\mathbf{N}/\mathsf{\mu }\mathbf{m})$ ${\mathit{h}}_{\mathit{f}}$ = 16 μ$\mathbf{m},{\mathit{h}}_{\mathit{r}}=20\mathsf{\mu }\mathbf{m}$	${\mathit{k}}_{\mathit{n}\mathit{o}\mathit{s}\mathit{e}}(\mathbf{N}/\mathsf{\mu }\mathbf{m})$ ${\mathit{h}}_{\mathit{f}}$ = 19 μ$\mathbf{m},{\mathit{h}}_{\mathit{r}}=23\mathsf{\mu }\mathbf{m}$
0.3	2.87	3.51
0.5	8.01	7.50
0.7	12.9	10.6
0.9	17.1	13.1

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. Radial stiffness of the nose at different supply pressures and radial clearances; these values were obtained with a 30 N load applied on the nose.

${\mathit{p}}_{\mathit{s}}$(MPa)	${\mathit{k}}_{\mathit{n}\mathit{o}\mathit{s}\mathit{e}}(\mathbf{N}/\mathsf{\mu }\mathbf{m})$ ${\mathit{h}}_{\mathit{f}}$ = 16 μ$\mathbf{m},{\mathit{h}}_{\mathit{r}}=20\mathsf{\mu }\mathbf{m}$	${\mathit{k}}_{\mathit{n}\mathit{o}\mathit{s}\mathit{e}}(\mathbf{N}/\mathsf{\mu }\mathbf{m})$ ${\mathit{h}}_{\mathit{f}}$ = 19 μ$\mathbf{m},{\mathit{h}}_{\mathit{r}}=23\mathsf{\mu }\mathbf{m}$
0.3	1.32	1.29
0.5	4.02	3.79
0.7	6.59	5.36
0.9	8.82	6.52

lists these values.

4.5. Simplified Natural Frequencies

At null rotational speed, in case the damping effects are neglected, the cylindrical natural frequency is expressed by:

$${\omega }_{cyl}=\sqrt{\frac{k_f+k_r}{m_r}}$$

(19)

The conical natural frequency is given by:

$${\omega }_{con}=\sqrt{\frac{\left(k_f{l}_f^2+k_r{l}_r^2\right)}{I}}$$

(20)

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. Undamped natural frequencies at different supply pressures and radial clearances.

${\mathit{p}}_{\mathit{s}}$(MPa)	${\mathit{h}}_{\mathit{f}}$  = 16  μm,  ${\mathit{h}}_{\mathit{r}}$  = 20  μm		${\mathit{h}}_{\mathit{f}}$ = 19 μ$\mathbf{m},{\mathit{h}}_{\mathit{r}}=23\mathsf{\mu }\mathbf{m}$
	${\mathit{\omega }}_{\mathit{c}\mathit{y}\mathit{l}}$  (krpm)	${\mathit{\omega }}_{\mathit{c}\mathit{o}\mathit{n}}$  (krpm)	${\mathit{\omega }}_{\mathit{c}\mathit{y}\mathit{l}}\left(\mathbf{krpm}\right)$	${\mathit{\omega }}_{\mathit{c}\mathit{o}\mathit{n}}\left(\mathbf{krpm}\right)$
0.3	46.6	58.9	48.2	57.7
0.5	72.4	86.2	67.8	78.0
0.7	89.2	102.8	79.1	88.9
0.9	101.1	114.3	86.7	95.9

reports the undamped natural frequencies of the rigid rotor-bearings system at null rotational speed.

4.6. 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⁻⁸ 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⁻⁸ s was chosen for dynamic simulations, as it was found to be a good compromise between accuracy and computational time.

4.7. 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, ${\vartheta }_x$ and ${\vartheta }_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:

$$\begin{gathered} x\left(z=0,t=0\right) =0 , y\left(z=0,t=0\right) =0 \\ x\left(z=z_L,t=0\right) =0, y\left(z=z_L,t=0\right) =0 \\ \dot{x}\left(z=0,t=0\right) =0 , \dot{y}\left(z=0,t=0\right) =-0.01 \mathrm{m}/\mathrm{s} \\ \dot{x}\left(z=z_L,t=0\right) =0 , \dot{y}\left(z=z_L,t=0\right) =0 \end{gathered}$$

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.

$$\begin{gathered} x\left(z=z_G,t=0\right) =1 \mathsf{\mu }\mathrm{m}, y\left(z=z_G,t=0\right) =0 \\ {\vartheta }_x\left(t=0\right) =0, {\vartheta }_y\left(t=0\right) =0 \\ \dot{x}\left(z=z_G,t=0\right) =0 , \dot{y}\left(z=z_G,t=0\right) =-0.01 \mathrm{m}/\mathrm{s} \\ {\dot{\vartheta }}_x\left(t=0\right) =0 , {\dot{\vartheta }}_y\left(t=0\right) =0 \end{gathered}$$

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.

4.7.1. 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 ${\theta }_x$. The LDM was used, considering time intervals $\mathsf{\Delta }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.

Table A1. Damped frequencies and damping factors resulting from FFT analysis and from LDM; $p_s$ = 0.3 MPa.

${\mathit{h}}_{\mathit{f}}/{\mathit{h}}_{\mathit{r}}$ (μm)	$\mathit{T}$ (ms)	Initial Condition	Signal	${\mathit{\omega }}_{\mathit{d}}\mathbf{LDM}$ (krpm)	$\mathit{\zeta }\mathbf{LDM}$ (krpm)	${\mathit{\omega }}_{\mathit{d}, 1/2}$ (krpm)
16/20	5	1	$y_0$	55.556	0.345	49.80 82.03
			$y_L$	78.125	0.219
			$y_G$	55.046	0.361
			${\theta }_x$	56.250	0.325
	5	2	$y_0$	55.556	0.346	49.80 82.03
			$y_L$	78.125	0.220
			$y_G$	55.046	0.361
			${\theta }_x$	56.426	0.326
19/23	5	1	$y_0$	49.180	0.171	48.40 68.80
			$y_L$	69.971	0.151
			$y_G$	49.180	0.174
			${\theta }_x$	49.315	0.167
	5	2	$y_0$	49.180	0.171	48.00 71.60
			$y_L$	70.175	0.151
			$y_G$	49.180	0.174
			${\theta }_x$	49.315	0.167

,

Table A2. Damped frequencies and damping factors resulting from FFT analysis and from LDM; $p_s$ = 0.5 MPa.

${\mathit{h}}_{\mathit{f}}/{\mathit{h}}_{\mathit{r}}$ (μm)	$\mathit{T}$ (ms)	Initial Condition	Signal	${\mathit{\omega }}_{\mathit{d}}\mathbf{LDM}$ (krpm)	$\mathit{\zeta }\mathbf{LDM}$ (krpm)	${\mathit{\omega }}_{\mathit{d}, 1/2}$ (krpm)
16/20	5	1	$y_0$	74.534	0.119	73.60 100.8
			$y_L$	99.448	0.100
			$y_G$	73.846	0.115
			${\theta }_x$	75.710	0.125
	5	2	$y_0$	74.257	0.143	73.20 101.1
			$y_L$	99.291	0.078
			$y_G$	73.892	0.142
			${\theta }_x$	75.567	0.145
19/23	5	1	$y_0$	65.934	0.063	66.00 86.60
			$y_L$	86.331	0.027
			$y_G$	65.934	0.065
			${\theta }_x$	65.753	0.059
	5	2	$y_0$	65.934	0.063	66.00 86.60
			$y_L$	86.331	0.027
			$y_G$	65.934	0.065
			${\theta }_x$	65.753	0.059
	10	1	$y_0$	65.934	0.059	66.00 86.60
			$y_L$	80.488	0.073
			$y_G$	65.854	0.059
			${\theta }_x$	65.854	0.060
	10	2	$y_0$	65.934	0.059	66.00 86.60
			$y_L$	80.685	0.071
			$y_G$	65.854	0.059
			${\theta }_x$	65.854	0.060

,

Table A3. Damped frequencies and damping factors resulting from FFT analysis and from LDM; $p_s$ = 0.7 MPa.

${\mathit{h}}_{\mathit{f}}/{\mathit{h}}_{\mathit{r}}$ (μm)	$\mathit{T}$ (ms)	Initial Condition	Signal	${\mathit{\omega }}_{\mathit{d}}\mathbf{LDM}$ (krpm)	$\mathit{\zeta }\mathbf{LDM}$ (krpm)	${\mathit{\omega }}_{\mathit{d}, 1/2}$ (krpm)
16/20	5	1	$y_0$	87.72	0.053	87.6 113.4
			$y_L$	112.9	0.052
			$y_G$	87.21	0.052
			${\theta }_x$	88.50	0.057
	5	2	$y_0$	87.60	0.047	87.6 113.4
			$y_L$	113.2	0.042
			$y_G$	87.27	0.046
			${\theta }_x$	88.13	0.049
	10	1	$y_0$	87.464	0.064	87.6 113.4
			$y_L$	113.6	0.041
			$y_G$	87.71	0.034
			${\theta }_x$	114.13	0.053
	10	2	$y_0$	86.44	0.049	87.6 113.4
			$y_L$	113.3	0.047
			$y_G$	95.00	0.038
			${\theta }_x$	113.5	0.035
19/23	5	1	$y_0$	75.71	0.019	75.60 96.00
			$y_L$	95.24	0.010
			$y_G$	75.71	0.019
			${\theta }_x$	75.95	0.019
	5	2	$y_0$	75.76	0.019	75.60 96.00
			$y_L$	84.41	0.026
			$y_G$	75.76	0.019
			${\theta }_x$	75.85	0.019
	10	1	$y_0$	75.95	0.021	75.60 96.00
			$y_L$	95.74	0.031
			$y_G$	75.71	0.014
			${\theta }_x$	76.19	0.028
	10	2	$y_0$	75.76	0.020	75.60 96.00
			$y_L$	96.33	0.031
			$y_G$	75.85	0.015
			${\theta }_x$	75.76	0.028

and

Table A4. Damped frequencies and damping factors resulting from FFT analysis and from LDM; $p_s$ = 0.9 MPa.

${\mathit{h}}_{\mathit{f}}/{\mathit{h}}_{\mathit{r}}$ (μm)	$\mathit{T}$ (ms)	Initial Condition	Signal	${\mathit{\omega }}_{\mathit{d}}\mathbf{LDM}$ (krpm)	$\mathit{\zeta }\mathbf{LDM}$ (krpm)	${\mathit{\omega }}_{\mathit{d}, 1/2}$ (krpm)
16/20	5	1	$y_0$	97.297	0.024	97.40 122.0
			$y_L$	121.390	0.036
			$y_G$	97.297	0.024
			${\theta }_x$	97.826	0.025
	5	2	$y_0$	97.297	0.031	97.40 122.0
			$y_L$	122.245	0.023
			$y_G$	97.035	0.013
			${\theta }_x$	124.030	0.040
	10	1	$y_0$	97.335	0.025	97.40 122.0
			$y_L$	122.160	0.025
			$y_G$	97.222	0.024
			${\theta }_x$	97.561	0.025
	10	2	$y_0$	97.110	0.028	97.40 122.0
			$y_L$	122.450	0.019
			$y_G$	97.335	0.019
			${\theta }_x$	123.010	0.024
19/23	5	1	$y_0$	82.645	−0.001	82.80 102.6
			$y_L$	87.209	−0.019
			$y_G$	82.645	−0.001
			${\theta }_x$	82.873	−0.002
	5	2	$y_0$	83.102	0.001	82.80 102.6
			$y_L$	102.560	0.022
			$y_G$	82.418	−0.005
			${\theta }_x$	83.565	0.009
	10	1	$y_0$	82.759	−0.002	82.80 102.6
			$y_L$	85.511	−0.007
			$y_G$	82.759	0.002
			${\theta }_x$	82.854	−0.002
	10	2	$y_0$	82.854	0.001	82.80 102.6
			$y_L$	103.330	0.035
			$y_G$	82.664	−0.004
			${\theta }_x$	82.949	0.004

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 ${\vartheta }_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.

4.7.2. 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 frequency-stabilization 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]. Figure 14 and Figure 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. Modal parameters resulting from LSCEM.

${\mathit{p}}_{\mathit{s}}\left(\mathbf{MPa}\right)$	${\mathit{h}}_{\mathit{f}}/{\mathit{h}}_{\mathit{r}}$	Signal	${\mathit{\omega }}_{\mathit{d},1}$ (krpm)	${\mathit{\zeta }}_1$ (-)	${\mathit{\phi }}_1$ (deg)	${\mathit{A}}_1$ (-)	${\mathit{\omega }}_{\mathit{d},2}$ (krpm)	${\mathit{\zeta }}_2$ (-)	${\mathit{\phi }}_2$ (deg)	${\mathit{A}}_2$ [-]
0.3	16/20	$y_L$	80.51	0.1675	−19.07	0.0737	n.a.	n.a.	n.a.	n.a.
		$y_0$	58.28	0.3427	−54.14	0.9738	n.a.	n.a.	n.a.	n.a.
0.3	19/23	$y_L$	70.00	0.1556	−18.80	0.1224	n.a.	n.a.	n.a.	n.a.
		$y_0$	50.00	0.2044	−70.50	1.115	n.a.	n.a.	n.a.	n.a.
0.5	16/20	$y_L$	100.6	0.0804	−21.40	0.0618	n.a.	n.a.	n.a.	n.a.
		$y_0$	74.82	0.1394	−76.58	0.7109	n.a.	n.a.	n.a.	n.a.
0.5	19/23	$y_L$	86.97	0.0620	−21.69	0.0710	66.22	0.0667	100.1	0.0377
		$y_0$	65.86	0.0629	82.46	0.7587	n.a.	n.a.	n.a.	n.a.
0.7	16/20	$y_L$	113.3	0.0365	21.8	0.029	87.83	0.0604	−50.2	0.0135
		$y_0$	87.69	0.0652	115.17	0.3245	n.a.	n.a.	n.a.	n.a.
0.7	19/23	$y_L$	96.09	0.0356	−23.30	0.0496	75.80	0.0193	79.62	0.0218
		$y_0$	75.81	0.0166	−85.73	0.6339	n.a.	n.a.	n.a.	n.a.
0.9	16/20	$y_L$	122.3	0.0124	4.855	0.0251	97.46	0.0281	−14.7	0.0104
		$y_0$	97.30	0.0265	157.25	0.3423	n.a.	n.a.	n.a.	n.a.
0.9	19/23	$y_L$	102.5	unst.	unst.	unst.	82.03	unst.	unst.	unst.
		$y_0$	82.03	unst.	unst.	unst.	unst.	unst.	unst.	unst.

, together with their phase shift and modal constants.

4.7.3. 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.

5. 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.

6. Conclusions

This paper describes the experimental and numerical investigation of an electro-spindle 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.