Multi-Field Coupling Dynamics Modeling of Aerostatic Spindle

The aerostatic spindle in the ultra-precision machine tool shows the complex multi-field coupling dynamics behavior under working condition. The numerical investigation helps to better understand the dynamic characteristics of the aerostatic spindle and improve its structure and performance with low cost. A multi-field coupling 5-DOF dynamics model for the aerostatic spindle is proposed in this paper, which considers the interaction between the air film, spindle shaft and the motor. The restoring force method is employed to deal with the times varying air film force, the transient Reynolds equation of the aerostatic journal bearing and the aerostatic thrust bearing is solved using ADI method and Thomas method. The transient air film pressure of aerostatic bearings is obtained which clearly presents the influence induced by the tilt motion of the spindle shaft. The motion trajectory of the spindle shaft is obtained which shows different stability of the shaft under different external forces. The dynamics model shows good performance on simulating the multi-field coupling behavior of the aerostatic spindle under external force. which is quite meaningful and useful for the further research on the dynamic characteristics of the aerostatic spindle.


Introduction
The aerostatic spindle plays a key role in an ultra-precision machine tool in the nanoprecision machining, which directly affects the machining quality [1,2]. Generally, the shaft of the aerostatic spindle is directly driven by a permanent magnet synchronous motor (PMSM) and supported by the aerostatic bearing, including the journal bearing and the thrust bearing [3].
Compared to the aerodynamic bearing, the aerostatic bearing provides higher capacity and steady force, while it contains more complex mechanism contributes to the dynamic behavior of aerostatic spindle. During the operation of the aerostatic spindle, it shows the multi-physics field coupling characteristics, in which the coupling effect of air bearing, the shaft and the motor should be considered for fully understanding its dynamics behavior [4]. However, the previous literatures commonly focused on single factor or two factors while investigating the dynamic characteristics of the aerostatic spindle [5][6][7][8]. It has been pointed out that the existence of unbalanced magnetic force (UMF) caused by rotor eccentricity at the motor can have a considerable impact on the machined surface [9][10][11], thus, the electromagnetic factor is not negligible in the dynamics model of aerostatic spindle.
Rotor dynamics modeling is the preferred method to study the dynamic characteristics of a rotor-bearing system, by which we can get detailed information of how the rotor acts under different condition. Among the previous literature, the rotor dynamics modeling method was widely used for modeling the aerodynamic bearing-rotor system with 2-DOF (Degree of Freedom) [12][13][14][15][16], the rotor trajectory and the stability problem were analyzed.
Since performance superiority of aerostatic spindle emerges in recent years, the analysis on aerostatic spindle came out frequently and have made great progress. Wang [17] implemented the rotor dynamic model of aerostatic bearing-rotor systems which applied a spherical aerostatic bearing, and the nonlinear dynamic behavior of the rotor bearing system is analyzed. Similar works were done by Zhang et al. [18], they proposed forecast orbit method to deal with the transient gas lubricated Reynold equation, and the 2D rotor trajectory under different rotor speed and imbalance mass was obtained. Two-dimensional analysis can reflect the characteristics of aerostatic bearing to some extent but not as precise as the result of the 5-DOF model. In a 5-DOF model, the tilt motion, the displacement in axial direction and the conical movement of the shaft can be acquired while the 2-DOF model cannot, and these forms of movement are also proved to have significant influence on the machined surface quality by Zhang et al. [19]. To achieve this, Li et al. [20,21] presented the 5-DOF model of an aerostatic spindle in a fly-cutting machine, in which the air pressure distribution of journal bearing and thrust bearing was calculated and the displacement of the shaft in different directions was obtained. Xu and Jiang [22,23] analyzed the 5-DOF rotor dynamics model of an aerostatic spindle, in which the stability, the unbalance response, and the forced response of the rotor-bearing system were investigated. Both Li and Xu adopted the dynamic coefficient method in their model, however it may lose the numerical accuracy when taking this method which indeed largely improved the efficiency. By contrast, the restoring force commonly used in 2-DOF models aims to obtain force directly by solving the transient Reynold equation [24,25], and it can also obtain the torque in a 5-DOF model. It can be optimized of the dynamic coefficient method, but if this measure is taken, massive extra work needs to be done. Besides, the available research on aerostatic spindle modeling fail to consider the electromagnetic factor, which has important effect on the motion error of the spindle.
In this paper, the electromagnetic factor is considered in the modeling of the aerostatic spindle, and the restoring force method is adopted. The finite difference method is adopted to discrete the transient Reynolds equation. The simulation for the motion trajectory of the spindle shaft is realized under the coupling of the restoring force of the aerostatic bearing, the UMF of PMSM and the external force on the shaft, where the shaft is regarded as a rigid rotor.

Mathematical Modeling of Aerostatic Spindle
The typical structure diagram of the aerostatic spindle is shown in Figure 1, the spindle shaft is directly driven by the motor, and the shaft is supported by an aerostatic journal bearing and two aerostatic thrust bearing.
Effectively modeling the dynamic behavior of the air film is quite necessary for the systematic and consistent air bearing design [26,27], thus, the modeling of the air film is of primary importance. The pressure distribution of the aerostatic bearing air film can be described by the Reynolds equation, and the Reynolds equation is obtained based on the Navier-Stokes Equations with following assumptions, the flow is isothermal, the gas viscosity is assumed to be constant, the pressure distribution in the vertical direction to the air film is assumed to be constant, the viscosity force is assumed to be much larger than inertia force, there is no velocity slip at the boundary, and the air is the ideal gas. Micromachines 2021, 12, x 3 of 23 Effectively modeling the dynamic behavior of the air film is quite necessary for the systematic and consistent air bearing design [26,27], thus, the modeling of the air film is of primary importance. The pressure distribution of the aerostatic bearing air film can be described by the Reynolds equation, and the Reynolds equation is obtained based on the Navier-Stokes Equations with following assumptions, the flow is isothermal, the gas viscosity is assumed to be constant, the pressure distribution in the vertical direction to the air film is assumed to be constant, the viscosity force is assumed to be much larger than inertia force, there is no velocity slip at the boundary, and the air is the ideal gas.

Modeling of Aerostatic Journal Bearing
The transient pressure distribution of aerostatic journal bearing can be modeled by the transient Reynolds equation in the Cartesian coordinate [28].
where x is the circumferential coordinates of journal bearing, z is the axial coordinate of journal bearing, h is the air film thickness, p is the air film pressure, μ is the dynamic viscosity of the air, 0 U is the surface velocity of the shaft, v  is the flow velocity of the supply air at the orifice, t represents the time, ρ is the air density, a P and a ρ are the pressure and density of the ambient air, The dimensionless parameters are defined as follows.

Modeling of Aerostatic Journal Bearing
The transient pressure distribution of aerostatic journal bearing can be modeled by the transient Reynolds equation in the Cartesian coordinate [28].
where x is the circumferential coordinates of journal bearing, z is the axial coordinate of journal bearing, h is the air film thickness, p is the air film pressure, µ is the dynamic viscosity of the air, U 0 is the surface velocity of the shaft, v is the flow velocity of the supply air at the orifice, t represents the time, ρ is the air density, P a and ρ a are the pressure and density of the ambient air, δ k = 1 at the orifice and δ k = 0 at other position. The dimensionless transient Reynolds equation is given by The dimensionless parameters are defined as follows.
where L j is the characteristic length of journal bearing, R is the radius of the journal bearing, C r is the radial clearance between the journal bearing and the shaft, P s is the supply pressure, Λ J is the bearing number of the journal bearing, ω is the rotating speed of the shaft. When the eccentricity or the tilt of the shaft happens as depicted in Figure 2, the thickness of the air film changes simultaneously, and the distribution of the air film thickness in the journal bearing can be expressed as where ε (i) is the eccentricity at the node i, θ (i,j) is the angle of air film at node i and j, α (i) is the eccentric angle at the node i, i b is the node of the barycenter, L s is the axial length of the shaft, M s is the number of nodes that L s divided into, ε ib is the eccentricity at the node i b , β is the tilt angle of the shaft, e ibx and e iby are the eccentricity of the node i b in x direction and y direction respectively, e ix and e iy are the eccentricity of the node i in x direction and y direction respectively. The mass flow rate can be acquired based on the isothermal assumption given by Powell [29].
where φ is the coefficient of the mass flow rate, inlet A is the minimum area of the flow The mass flow rate can be acquired based on the isothermal assumption given by Powell [29].
where φ is the coefficient of the mass flow rate, A inlet is the minimum area of the flow channel, d is the diameter of the orifice, β i = p/P s and β k = (2/(k + 1)) k/(k−1) . As shown in Figure 3, the computational domain of the journal bearing film is meshed into θ 0 : N j and Z 0 : M j , the Periodic boundary condition, the Atmospheric boundary condition and the Mass flow boundary condition are also defined.  Submitting the expression 2 S P = , the central difference method is adopted to discrete the variable in the θ and Z directions. The ADI method [18] is employed to simplify the implicit equation obtained after the discretization. The ADI format of the equation can be expressed as follows, Submitting the expression S = P 2 , the central difference method is adopted to discrete the variable in the θ and Z directions. The ADI method [18] is employed to simplify the implicit equation obtained after the discretization. The ADI format of the equation can be expressed as follows, The increment in the marching direction Z is carried out at the time step of n + 1, and the increment in the marching direction θ is carried out at the time step of n + 2. The Equation (10) and Equation (15) can be solved by the Thomas method, Then, the pressure distribution of the aerostatic journal bearing is obtained. Base on the load formulas given by Rowe [30], the air film force (i.e., the restoring force of the journal bearing) of the journal bearing in x and y direction is given by And the torque of the air film force (i.e., the restoring torque of the journal bearing) on the barycenter of the spindle shaft respect to the x axis and y axis is given by where θ is the circumferential coordinates of journal bearing, r is the radial coordinate of journal bearing, V 0 is the surface velocity of the shaft.

The dimensionless transient Reynolds equation is given by
The dimensionless parameters are defined as follows.
where L t is the characteristic length of journal bearing, and Λ T is the bearing number of the thrust bearing. The tilt motion of the shaft is considered to have influence on the air film thickness of the thrust bearing while the influence of the eccentricity is neglected. The tilt motion at the thrust bearing is shown in Figure 4, and the distribution of the air film thickness at the thrust bearing can be expressed as where θ (i,j) is the angle of air film at node i and j, γ (i) is the angle between the tile direction and the positive direction of x at the node i, i is the node number along the radial direction, R a is the inner diameter of the thrust bearing, R b is the out diameter of the thrust bearing, M t is the number of thrust bearing nodes along the radial direction, ∆z is the shaft displacement in the axial direction. L is the characteristic length of journal bearing, and T Λ is the bearing number of the thrust bearing. The tilt motion of the shaft is considered to have influence on the air film thickness of the thrust bearing while the influence of the eccentricity is neglected. The tilt motion at the thrust bearing is shown in Figure 4, and the distribution of the air film thickness at the thrust bearing can be expressed as (26) where ( )   As shown in Figure 5, the computational domain of the thrust bearing film is meshed into θ(0 : N t ) and R r (0 : M t ), the Periodic boundary condition, the Atmospheric boundary condition and the Mass flow boundary condition are also defined.  Adopt the same treatment as the journal bearing, the ADI format of the equation obtained from thrust bearing can be expressed as Adopt the same treatment as the journal bearing, the ADI format of the equation obtained from thrust bearing can be expressed as By solving the Equation (27) and (32) using Thomas method, the pressure distribution of the thrust bearing can be obtained. The air film force (i.e., the restoring force of the thrust bearing) of the thrust bearing in z direction is given by And the torque of the air film force (i.e., the restoring torque of the thrust bearing) of the thrust bearing respect to the x axis and y axis is given by

Modeling of PMSM
Ideally, the magnetic force of PMSM is symmetric. However, the UMF will come out when the rotor eccentricity happens or the structure of the PMSM is not symmetric. Kawase et al. [31] used 3-D finite element method to analyze the UMF of the PMSM, and it is found that the axial component of the UMF is relatively small compared to other two components. Here, we assume the structure of the PMSM is symmetric and the axial component of the UMF is negligible. According to Maxwell stress tensor method, the 2D magnetic forces in the radial and tangential direction can be expressed as follows [32].
where B r is the radial flux density, B θ is the tangential flux density, µ 0 is the permeability of the air. Figure 6 shows the schematic structure of motor with eccentric rotor [4], and the 2-D magnetic force in cartesian coordinates can be expressed as To simulate the 3-D state of the PMSM, the multiple slice method is employed. As shown in Figure 7, the magnetic force of the 3-D PMSM can be expressed as And the magnetic torque of the 3-D PMSM is given by   To simulate the 3-D state of the PMSM, the multiple slice method is employed. As shown in Figure 7, the magnetic force of the 3-D PMSM can be expressed as

Dynamics Modeling of ABMS
The spindle shaft is regarded as a rigid rotor in the model, and the spindle shaft is considered to move with 5-DOF. The dimensionless acceleration, velocity and the attitude in the space can be calculated by

Dynamics Modeling of ABMS
The spindle shaft is regarded as a rigid rotor in the model, and the spindle shaft is considered to move with 5-DOF. The dimensionless acceleration, velocity and the attitude in the space can be calculated by where the dimensionless parameters are given as follows.
where F ex , F ey , F ez are the dimensionless external force applied on the spindle shaft in x, y, z direction, F bx , F by , F bz are the dimensionless air film force in x, y, z direction, F mx , F my are the dimensionless UMF in x, y direction, T ex , T ey are the dimensionless external torque with respect to x, y axis, T bx , T by are the dimensionless air film torque with respect to x, y axis, T mx , T my are the dimensionless magnetic torque with respect to x, y axis, M is the dimensionless mass of the spindle shaft, I x , I y , I z are the rotational inertia of the spindle shaft with respect to x, y, z axis, X, Y, Z are the dimensionless displacement in x, y, z direction, Θ X , Θ Y are the rotation angle of the spindle shaft with respect to x, y axis.
The spindle shaft is still at the initial condition, thus . δ = {0} at the initial time. The flow chat for calculating the 5-DOF dynamics model of the ABMS is shown in Figure 8.

Detailed Parameter
A case study is conducted to verify the effectiveness of the dynamics model of the ABMS. Table 1 lists the detailed parameters of the simulated model.

Detailed Parameter
A case study is conducted to verify the effectiveness of the dynamics model of the ABMS. Table 1 lists the detailed parameters of the simulated model.

Numerical Result
The model selects the same motor parameter as that in reference [4], and the variation trend of UMF has been given as follows.
f mx = −4.4ε 2 + 7.4ε + 0.3 · sin 2πω 60 t + 9.4ε 2 − 3.1ε + 3.8 · sin θ r (ε ≤ 0.5) −4.4ε 2 + 7.4ε + 0.3 · sin 2πω 60 t + (7.0ε + 1.1) · sin θ r (ε > 0.5) (53) Two cases are calculated, i.e., the dynamics model with 400 N and 50 N external force. In the first case, 400 N external force is applied on the shaft end in x direction, the default dimensionless time is set to be 100. Figure 9a shows the air pressure distribution of the aerostatic journal bearing, Figure 9b and c show the air pressure distribution of the front and rear aerostatic journal bearing respectively, And Figure 9d shows the integrated air pressure distribution vision of aerostatic bearing with 400 N external force at dimensionless time 100.  The result shows that when applies the external load on the shaft end, the spindle shaft can have a certain tilt angle and it affects the air pressure distribution of the aerostatic bearings obviously, which means the restoring force and torque of the air film applied on the spindle shaft are also changed. The tilt angle of the spindle shaft versus x axis and y axis is shown in Figure 10a,b, the tilt angle is converged and varies periodically with a certain amplitude, the extent of the tilt angle versus y axis is larger than the tilt angle versus x axis, this is mainly because the external force is applied in the x direction, thus the external torque is applied on y axis. Figure 11a,b shows the motion trajectory of the shaft barycenter and the shaft end respectively. The result shows that the shaft acts stable with 400 N external force, its trajectory in x-y plane converged to a certain region. From a spatial perspective, the shaft end has left its initial position towards its balance position as shown in Figure 11c, and the displacement of the shaft end in z direction is shown in Figure 11d. The displacement in z direction is also converged and it varies periodically versus time.
(c) (d)   Figure 11a,b shows the motion trajectory of the shaft barycenter and the shaft end respectively. The result shows that the shaft acts stable with 400 N external force, its trajectory in x-y plane converged to a certain region. From a spatial perspective, the shaft end has left its initial position towards its balance position as shown in Figure 11c, and the displacement of the shaft end in z direction is shown in Figure 11d. The displacement in z direction is also converged and it varies periodically versus time. In the second case, 50 N external force is applied on the shaft end in x direction. The result is quite different with the first case, Figure 12a,b shows the tilt angle of the shaft versus x axis and y axis which are both diverged, which means the shaft is unstable under the current condition. In the second case, 50 N external force is applied on the shaft end in x direction. The result is quite different with the first case, Figure 12a,b shows the tilt angle of the shaft versus x axis and y axis which are both diverged, which means the shaft is unstable under the current condition.  Figure 13a shows the corresponding air pressure distribution of the aerostatic journal bearing, Figure 13 b,c show the air pressure distribution of the front and rear aerostatic journal bearing respectively, Figure 13d shows the integrated air film pressure distribution vision of aero-static bearing with 50 N external force at dimensionless time 100. According to the result, the air film pressure distribution of the aerostatic bearing becomes quite uneven, both the aerostatic journal bearing and the aerostatic thrust bearing show severe aerodynamic effect. Figure 14a,b shows the motion trajectory of the shaft barycenter and the shaft end respectively. The result shows that the trajectory of shaft barycenter does not show the divergent trend while the trajectory of shaft end does, which means that the spindle shaft is at conical unstable state. It can be observed clearer in a spatial perspective as shown in Figure 14c, the trajectory of the shaft end moves as a spiral path. The displacement of the shaft end in the axial direction is also divergent as shown in Figure 14d.  Figure 13a shows the corresponding air pressure distribution of the aerostatic journal bearing, Figure 13b,c show the air pressure distribution of the front and rear aerostatic journal bearing respectively, Figure 13d shows the integrated air film pressure distribution vision of aero-static bearing with 50 N external force at dimensionless time 100. According to the result, the air film pressure distribution of the aerostatic bearing becomes quite uneven, both the aerostatic journal bearing and the aerostatic thrust bearing show severe aerodynamic effect. Figure 14a,b shows the motion trajectory of the shaft barycenter and the shaft end respectively. The result shows that the trajectory of shaft barycenter does not show the divergent trend while the trajectory of shaft end does, which means that the spindle shaft is at conical unstable state. It can be observed clearer in a spatial perspective as shown in Figure 14c, the trajectory of the shaft end moves as a spiral path. The displacement of the shaft end in the axial direction is also divergent as shown in Figure 14d.

Discussion
In the motion process of the spindle shaft, the external force applied on the shaft end produce the torque on the shaft, resulting in the tilt motion of the shaft, then the distribution of the motor flux and the air film pressure is changed. Under the multi-field coupling effect of the external force, the air film force and the magnetic force, the shaft may finally

Discussion
In the motion process of the spindle shaft, the external force applied on the shaft end produce the torque on the shaft, resulting in the tilt motion of the shaft, then the distribution of the motor flux and the air film pressure is changed. Under the multi-field coupling effect of the external force, the air film force and the magnetic force, the shaft may finally stabilize in a certain spatial attitude or become unstable and even hit the bearing sleeve.
In this paper, a 5-DOF dynamics model of the aerostatic spindle was implemented using the restoring force method, while previous literatures mainly used the dynamic coefficient method [23,33]. The dynamic coefficient method improves the calculation speed compared to the restoring force method, but the dynamic behavior of the shaft can be obtained directly by employing the restoring force method. The proposed 5-DOF dynamic model of aerostatic spindle has considered the influence of the tilt motion and axial displacement of the shaft, which is neglected in the 2-DOF model. In the 2-DOF model of aerostatic spindle only the external force, the restoring force of the air film and the unbalanced magnetic force of the motor are considered, while the external torque, the restoring torque of the air film and the magnetic torque of the motor are not considered. However, they are all considered and integrated in the 5-DOF model. The motion trajectory of the shaft center in the 2D plane characterizes the dynamic behavior of the shaft, but it cannot show the conical motion of the shaft with the external force and external torque, which is quite important because the shaft has the spatial motion with 5-DOF.
In the simulation result, the dynamic behavior of the spindle shaft is stable with 400 N force and unstable with 50 N force. Similar phenomenon has been observed in the research of journal bearings, such as the research done by Khonsari et al. [34]. The stability boundary of the journal bearing varies with the change of the Sommerfeld Number, and the load on the bearing is one of the parameters that determine the value of the Sommerfeld Number. The stability boundary expands with the increase of the external load. To our knowledge, we think the reason behind dynamic behavior is probably the same as that of the journal bearings. That is, when the external load increase, the initial position locates in the unstable region under the low external load may locate in the stable region due to the expanding of the stability boundary. The Sommerfeld Number is also related to the rotation velocity [34]. Thus, the rotation velocity also affects the stability of the shaft, which is probably because that the change of the rotation velocity results in different aerodynamic effect and affects the stability. The critical value for distinguishing the stability of the shaft can be obtained in the 2-DOF model as the works done by Yang et al. [24]. It is sure that the critical value for the conical stability of the spindle shaft in the 5-DOF model also exists. It is meaningful to study the stability criterion of the shaft with 5-DOF. However, it will be much difficult to obtain the critical value in the 5-DOF model compared to that in the 2-DOF model. The 5-DOF model is much time costly to solve due to its high complexity and difficulty, thus, it will be further explored in the future research.
The proposed model can be improved in somehow. The shaft is regarded is a rigid body rather than the flexible body, and the structural deformation is neglected, however, it may have potential influence on the dynamic behavior of the shaft. Laha et al. [35] investigated the stability of the rotor supported by air journal bearing considering the rotor flexibility, in which the object was 2-DOF model. If the shaft is considered as a flexible rotor in a 5-DOF model, extra calculations need to be done, and the improvement of the solving code is also required. Due to the limitation of the current experimental condition, the validation of the theoretical results will be conducted in the future work.

Conclusions
In this study, a 5-DOF dynamics model of the aerostatic spindle is established. The modeling method and the numerical solution is given. Through the numerical model, the transient pressure distribution of the air film as well as the motion trajectory of the spindle shaft is obtained as expected. The result shows the significant meaning of considering the influence of the shaft tilt on the air film pressure distribution, which further determines the restoring force and torque of the aerostatic bearing. Besides, the different conical stability behavior of the spindle shaft is shown under different external load through the model.