# Experimental and Numerical Dynamic Identification of an Aerostatic Electro-Spindle
## 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{ph^{3}}{12\mu RT}\frac{\partial p}{ \partial z}\right)\ +\ \frac{\partial}{r\partial\theta}\left(\frac{ph^{3}}{12\mu RT}\frac{ \partial p}{r\partial\theta}\right)\ +\ g_{in}=\frac{\omega r}{2 RT}\frac{\partial(ph)}{r \partial\theta}+\frac{1}{RT}\frac{d(ph)}{dt} \tag{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\Big{(}1-e^{-8.2\frac{h}{h}}\Big{)} \tag{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}+\Delta t\cdot f\Big{(}p^{t},\ p^{t+1}\Big{)} \tag{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}\)): \begin{table}
\begin{tabular}{c c c c c c} \hline \hline \(p_{s}\) (MPa) & Signal & \(\omega_{d}\) (krpm) & \(\zeta\) (\(-\)) & \(\phi\) (deg) & \(A_{1}\) (\(-\)) \\ \hline \multirow{2}{*}{0.3} & \(y_{2}\) & 73.44 & 0.1483 & 3.474 & 3.400 \\ \cline{2-6}  & \(y_{1}\) & 70.21 & 0.2578 & 169.8 & 7.963 \\ \hline \multirow{2}{*}{0.5} & \(y_{2}\) & 89.50 & 0.1072 & 9.420 & 5.305 \\ \cline{2-6}  & \(y_{1}\) & 87.23 & 0.1632 & 178.0 & 7.673 \\ \hline \multirow{2}{*}{0.7} & \(y_{2}\) & 99.75 & 0.00762 & \(-\)3.336 & 5.223 \\ \cline{2-6}  & \(y_{1}\) & 95.46 & 0.1060 & 153.6 & 8.461 \\ \hline \multirow{2}{*}{0.9} & \(y_{2}\) & 103.8 & 0.07330 & 19.69 & 3.150 \\ \cline{2-6}  & \(y_{1}\) & 107.0 & 0.09230 & \(-\)154.6 & 4.102 \\ \hline \hline \end{tabular}
\end{table} 
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. \[h(z,\theta)\ =h_{0}-\ \left[x_{G}+\theta_{y}(z-z_{G})\right]\ \cos\theta-\ \left[y_{G}- \theta_{x}(z-z_{G})\right]\sin\theta \tag{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\dot{\theta}}+\mu\frac{\omega r}{h} \tag{9}\]  
The 4 DOF equations of motion for a rigid rotor are considered:  
\[ \begin{array}{c}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{\theta}_{x}+I_{p}\omega\dot{\vartheta}_{y}-M_{cx}-F_{ext,\ y}(z_{G}-z_{ nose})\ =\ \left(I_{p}-I\right)\gamma\omega^{2}\sin(\omega t+\varphi_{1})\\ I\ddot{\theta}_{y}-I_{p}\omega\dot{\vartheta}_{x}-M_{cy}+F_{ext,\ x}(z_{G}-z_{ nose})\ =\ \left(I-I_{p}\right)\omega^{2}\gamma\cos(\omega t+\varphi_{1})\end{array}