Evaluation of the Inﬂuence of Shaft Shape Errors on the Rotation Accuracy of Aerostatic Spindle—Part 1: Modeling

: In the research on the aerostatic spindle, for quantifying the inﬂuence of shaft shape errors on its rotation accuracy, the calculation models of unbalanced air ﬁlm force (UAFF) and centrifugal inertia force (CIF) are established. Based on the typical shape of the shaft, three composite cylinder models are constructed. The data sources are selected from the cylindricity measurement results of the shaft and substituted into the model for calculation. The results show that, for the shaft in this paper, when the roundness deviation is ≤ 1.25 µ m and the cylindricity deviation is ≤ 2 µ m, the inﬂuence of shape errors on the rotation accuracy of the shaft is small enough. The harmonic analysis of the data source shows that the magnitude and stability of the UAFF on the shaft are determined by the actual shape of the shaft’s cylindrical surface, and cannot be inferred only by the value of shape errors.


Introduction
The aerostatic spindle uses compressed air to fill the gap between the bearing and the shaft to make the shaft float inside the bearing. Compared with the traditional spindle, it has the advantages of high rotation accuracy, low noise and frictionless [1]. The rotation accuracy of the aerostatic spindle is mainly affected by the shape error of the shaft [2], which will make the air film thickness change continuously when the shaft rotates. It can be seen from the Reynolds equation that the change in air film thickness will produce pressure fluctuation and cause UAFF. It is one of the key factors influencing spindle rotation accuracy and has a significant impact on the surface processing quality of parts [3][4][5][6].
Many researchers have studied the influence of the deviations from ideal surface topography of the shaft and bearing on the performance of the aerostatic spindle and bearing. The results of a theoretical analysis by Pande and Somasundaram [7] indicated that manufacturing errors will significantly affect the performance characteristics of aerostatic journal bearings. Stout et al. [8] studied the shape error of the bearing, and the results show that the parallelism, roundness and tilt of the bearing will affect its load capacity. Bhat et al. [9] analyzed the influence of surface profiling errors on the performance of a flat pad aerostatic bearing by using the Pareto optimization method, and concluded that multi-orifice aerostatic flat pad bearings are highly sensitive to surface profile variations and the surface profiles of bearings should be accurately controlled during machining. Through numerical and experimental research, Cui et al. [10,11] found that the bearing characteristics were significantly influenced by the manufacturing error, and the radial running accuracy of journal bearings can be improved by reducing the waviness amplitude or spatial wavelength. Wang et al. [12,13] studied the axial and circumferential surface waviness of aerostatic journal bearings by using the finite difference method and linear perturbation theory, and the results show that the amplitude of surface waviness has a significant impact on the capacity, stiffness and damping of bearings. Sun et al. [14] studied the effect of shaft shape errors on dynamic characteristics of a rotor-bearing system and the quantitative relationship between shaft shape errors and bearing reaction forces of the rotor-bearing system was obtained. Dimofte [15,16] pointed out that the three-wave bearing shows a better load capacity than the three-wave-groove bearing and the three-lobe bearing at any applied load and running regime, but the three-wave bearings are mare sensitive to the direction of the applied load than the other bearings. Charitopoulos et al. [17] analyzed the effect of manufacturing errors on the load capacity and friction coefficient of the textured micro-thrust bearings and revealed that the performance of textured micro-thrust bearings could be improved by manufacturing errors. In addition to the shape error, assembly faults will also significantly affect the spindle performance, as angular misalignment deteriorates the pressure distribution between the shaft and the bearing surfaces, and an air bearing with angular misalignment is more sensitive to the dynamic motions of the shaft in the bearing [18].
It is generally believed that in order to obtain optimal performance of aerostatic bearings, the air film thickness should usually be 5-20 µm, and the shape error of the mating surface is required to be better than 1/10 of the air film thickness, which presents a challenge for the manufacturing capability [19,20], and the precision of the shaft is much more important than that of the bearing [21]. In the present research of the aerostatic spindle, the shape error of the shaft and bearing is mostly constructed by a trigonometric function, which is quite different from the actual situation. When the spindle rotates, the topography of any position on the shaft surface will affect the air film distribution, and then affect the rotation accuracy, but the shape errors, such as waviness, roundness and cylindricity, are unable reflect the topography of the shaft. In this paper, the calculation models of UAFF and CIF considering shaft shape errors are established, and then the representative shaft roundness error data selected from the measurement database of a cylindricity instrument are used for calculation. The analysis results give the recommended roundness and cylindricity deviations to ensure shaft rotation accuracy, and provide theoretical support for the influence of shape errors on the spindle rotation accuracy. The analysis and calculation method in this paper will be a useful reference in the design, manufacture and test process of the aerostatic spindle.

Theoretical Model
According to the position of the thrust surface, the shaft of the aerostatic spindle can be divided into 'T' type, '+' type and 'H' type. The most common and stable is the 'H' type shaft, as shown in Figure 1. Thrust surfaces are distributed at both ends of the shaft, which can minimize the axial movement and tilt of the shaft during rotation. Assuming that the thrust surface of the shaft is an ideal shape, the axial movement and tilt of the shaft can be ignored.
This spindle is made for precision micro-grinding, which has high requirements for rotation accuracy. The journal bearings and thrust bearings adopt multiple orifices evenly arranged to provide sufficient load capacity. As the continuous working time of the spindle is short, it is cooled by air flow. The exhaust air at the front end of the spindle is discharged from the outer edge of the thrust bearing, and the exhaust air at the rear end of the spindle is collected to cool the motor and then discharged. The motor is an AC frameless motor, the motor rotor is installed at the rear end of the shaft and the stator is fixed in a special frame (not shown in Figure 1).  Figure 2a is a schematic diagram of the spindle considering the radial shape error of the shaft. It can be seen that the shape error of the shaft surface directly affects the uniformity of the film thickness distribution between the bearing and the shaft. A two-dimensional closed cylindrical coordinate system x o z is established on the inner surface of the bearing, where the x -axis is at the bottom of the bearing inner surface and the z -axis is parallel to the z-axis in the Cartesian coordinate system, then the film thickness and pressure at any position on the shaft surface can be indexed with h x z and p x z . The general form of the Reynolds equation applicable to this coordinate system is [22] ∂ ∂x

Model of the UAFF
(1) Some assumptions are as follows: a. The bearing is stationary, i.e., u b = 0, w b = 0. b. The shaft is rigid and without axial movement, i.e., w a = 0. c. The spindle is fully cooled and without local expansion, i.e., h(∂ρ/∂t) = 0. d. The air is ideal gas, which can be obtained according to the equation of state that is e. The shaft is in fixed eccentric motion. Based on the above assumptions, the dynamic error of the spindle is only the rotation error.
Let u = u a , add the inflow term V in at the orifice entrance [23] and Equation (1) can be changed into ∂ ∂x The first term on the right side of Equation (2) is the Couette term, which can be divided into the pressure wedge term hu(∂p/∂x ), stretch term ph(∂u/∂x ) and physical wedge term pu(∂h/∂x ), in which the stretch term of the rigid shaft is zero. When there are shape errors on the shaft surface (see Figure 2), the physical wedge term is not zero, therefore, Equation (2) can be arranged as By employing the dimensionless parameters Equation (4) can be written as: Divide the shaft into m grids along the x -direction and n grids along the z -direction. Rewrite Equation (5) by using the finite difference method (see Equations (A1)-(A10) in Appendix A) [24], and the following equation can be obtained Using Equation (6), the surface pressure P i,j at node (i, j) of the shaft can be calculated. The action range of P i,j is the area surrounded by 9 nodes in Figure 3, and the dimensionless area is ∆X ∆Z .
The resultant force of the air film force on the j-th section of the shaft in the Cartesian coordinate system is F air where The total UAFF on the shaft is

Model of the CIF
Due to the shape errors of the shaft, the actual centroid of the shaft will deviate from the ideal position, and the CIF will be generated when the shaft rotates. Divide the shaft into n segments along the z-axis and, when the n is large enough, each segment can be regarded as a slice. Divide the j-th slice into m sectors, as shown in Figure 4, where the included angle of each sector is 2π/m and the radius of the i-th sector is r i,j , and if the area density of the sheet is µ density , then the mass of the i-th sector is The centroid coordinates of the i-th sector arē The total centroid coordinate of the j-th slice is The total centroid coordinate of the shaft is The mass of the shaft is Assuming that the axis of rotation coordinate of the shaft is (x rotation , y rotation ), the CIF can be calculated by the following equation: Then the resultant force on the shaft is

Construction of Shaft Cylinder
The shape of the shaft cylinder can be divided into four categories: conical (Figure 5a), concave (Figure 5b), convex ( Figure 5c) and banana-shaped (Figure 5d). The conical, concave and convex shapes are axisymmetric shapes and, when they rotate around the axis, the change in the radial film gap is only determined by the section shape, which has little impact on the stability of the shaft. When the banana-shaped shaft rotates, the centers of different sections will deviate from the axis, resulting in significant changes in the air film gap, which has a great impact on the stability of the shaft. The actual shape of the shaft cylinder is composed of one or more of the four shapes. It is necessary to construct a shaft cylinder to analyze the influence of shape errors on rotation accuracy. In this paper, the shaft cylinder is constructed by scaling and moving the datum section following the characteristic line. The characteristic line in Figure 5a is a straight line, and the characteristic lines in Figure 5b-d are curves. In order to simplify the construction process, the curve is defined as an arc. The red dotted line in Figure 5 represents the actual shape of the datum section. Since the conical, concave and convex shapes are axisymmetric shapes, the generatrix is used as the characteristic line (the red solid line in Figure 6a and the red arcs in Figure 6b,c). The horizontal distance between any point of the generatrix and the z-axis is the nominal radius of the current section (see Equation (17) (a)-(c)), and the center of all sections is on the z-axis. In Figure 5d, the connecting line of the section center is used as the characteristic line (the red arc in Figure 6d), the x-direction coordinate of any point of the characteristic line can be expressed as Equation (18) and the shape of the section at different heights is the same.
x cylinder z The radius of the characteristic arc can be obtained by using the chord length and bow altitude method: The cylindricity is determined by three factors: the section's roundness, the average size difference of sections and the center position of sections. The latter two factors correspond to ∆r a ∆r b and ∆r c in Figure 6a-c, and ∆r d in Figure 6d, respectively. The ∆r d in Figure 6d represents the deviation of the section center from the axis of the cylinder. Through the measurement of a large number of shafts, it is found that the deviation exists in all shafts, but it is generally less than 0.1 µm. Therefore, three shaft cylinders are constructed by combining banana-shaped with conical, concave and convex shapes, respectively, then the cylindricity deviation calculated by

Selection and Processing of Data Sources
In previous studies, researchers mostly used trigonometric functions to construct shaft shape errors, and the shaft surface is very smooth and the error distribution is very regular, which is quite different from the actual shaft surface. In order to study the influence of shape errors on shaft rotation accuracy, the R10 series of preferred number is used to determine five roundness deviations as the datum according to the machining requirement of aerostatic shafts, which are 0.63 µm, 0.80 µm, 1 µm, 1.25 µm and 1.60 µm, respectively. Then, five corresponding roundness measurement data are selected from the measurement database of the cylindricity instrument as the data sources. The peakto-valley roundness deviations of these data sources are 0.58 µm, 0.79 µm, 1 µm, 1.24 µm and 1.58 µm, respectively. These values are based on Gaussian filtering and least squares evaluation. To ensure the diversity of analysis, the data sources are selected from different shafts. Figure 7 shows the roundness measurement diagram (linear and polar) of each data source. In order to facilitate comparison, each linear diagram adopts the same abscissa and the same ordinate. It can be seen that the shapes of each source are quite different. The polar diagram of Figure 7a is oval, and the other diagrams are multi-lobed.
Each data source is scaled according to five datums, and a total of twenty-five sets of data are obtained for the next calculation. The scaling method is shown in Equation (21), and the data source is in polar format.

Parameters of Cylinder Construction
Firstly, a cylinder with the same deviation of roundness and cylindricity is constructed to study the influence of roundness on shaft rotation accuracy separately, that is ∆r i = 0(i = a, b, c, d), which is constructed by using the five datums in the previous section.
Secondly, the cylinder with compound deviations is constructed. According to the machining requirements of aerostatic shafts, four cylindricity deviations are determined, which are 1.5 µm, 2 µm, 2.5 µm and 3 µm, respectively. Only three sources are selected to construct the cylinder in order to control the amount of calculation, which are Source-B, Source-C and Source-D, respectively. The value of each parameter can be calculated by using Equation (21), where ∆r d = 0.1 µm. See Table 1 for other parameters.

Settings of Numerical Calculation
The spindle has two rows of orifices, with 8 orifices in each row, see Table 2 for the specific parameters. To balance the calculation efficiency and accuracy, set 360 grids in the in the circumferential direction (x -direction) and 100 grids in the axial direction (z -direction), and the speed range is [1000, 15,000] r/min.

Effect of the Roundness
The twenty-five sets of roundness data are filtered [25,26] and imported into the simulation model for calculation, and the calculation results of UAFF (see Figure 8a-e) and CIF (see Figures 8f) on the shaft are obtained. It can be seen that for the shaft constructed based on the same data source, in general, the smaller the roundness deviation, the smaller and more stable the force on the shaft at different speeds. Comparing Figure 8a-e horizontally, it can be found that even if the roundness deviation values are the same, the force on the shaft constructed based on different sources will vary greatly. The peak-to-valley roundness deviation of Source-B (0.79 µm) is smaller than Source-C (1 µm) and Source-D (1.24 µm), but the UAFF is larger. The reason for these results is that the roundness deviation only represents the difference between the maximum and minimum values in the roundness measurement data, and cannot reflect the actual shape of the shaft, while the shaft shape directly affects the film thickness and then affects the UAFF. When the roundness error value is large, for example, 1.6 µm, the UAFF of the shaft changes greatly at different speeds, and the force will be very small at some speeds, such as 10,000 r/min and 12,000 r/min in Figure 8b and 9000 r/min and 14,000 r/min in Figure 8e. This is because the air in some parts cannot reach a stable state in time and cannot produce enough air film force at higher speeds [2]. The CIF is directly proportional to the square of the angular velocity (see Equation (15)), therefore, only the CIF at the speed of 15,000 r/min is calculated (see Figures 8f). It can be seen that the value of the CIF is very small under all settings, and the impact on the shaft rotation accuracy can be ignored.
By using Equation (22), the stiffness of the shaft can be calculated to evaluate the influence of the UAFF on the rotation accuracy of the shaft [27].
The main stiffness of the shaft is about 40 N/µm. If the shaft eccentricity caused by UAFF is limited to 0.05 µm, the UAFF shall be less than 2 N. As can be seen from Figure 8a-e, the UAFF on every shaft is less than 2 N. However, when the roundness deviation is 1.6 µm, the UAFF fluctuates greatly with the speed; when the roundness deviation is 0.63 µm, the UAFF on the shaft is the smallest and changes little with the speed, but it presents a challenge for the manufacturing capability. Therefore, the roundness deviation of the shaft can be selected as 0.8 µm, 1 µm and 1.25 µm.

Harmonic Analysis of Data Source
The radial measurement accuracy of the Taylor Hobson 585LT cylindricity instrument is ±0.01 µm, and the default number of sampling points is 18,000 points per revolution, so the data measured by the cylindricity instrument can reflect the surface topography of the shaft well. The topography can be decomposed into several fluctuations by harmonic analysis, and Figure 9 is the harmonic analysis diagram of Source-B under 1-50 undulations per revolution (UPR). It can be seen that the amplitude of the 1st harmonic is the largest, which is produced by the eccentricity of the shaft to the datum axis and has nothing to do with the shape of the shaft. The 1st harmonic will be suppressed during roundness analysis, as shown in Figure 7, which will not affect the measurement results. There are two reasons for the 2nd harmonic, one is that the axis of the shaft is not parallel to the cylindricity instrument datum axis, and the other is that the shaft is oval. If the shaft is machined by centreless grinding or held by an overtightened three jaw chuck, the 3rd harmonic will be produced [28]. Higher order harmonics can be induced by the effects of the tool on the part, such as machine chatter. Harmonic analysis is carried out on all sources in turn. After suppressing the first harmonic, draw a 100% stacked bar (see Figure 10) for the wave of 2-50 UPR according to its amplitude proportion. The numbers in the stacked bar are the amplitude proportion of the current wave, those whose amplitude proportion is less than 2% are not displayed. It can be observed from Figure 10 that the amplitude proportion of low-order harmonics is large and the high-order harmonics is small. The amplitude of the 2nd and 3rd harmonic of Source-A and Source-C accounts for a large proportion, while the shapes of Source-A and Source-C are relatively smooth as Figure 7a,c shows, so the change in air film gap will be relatively small when the shaft rotates, and the UAFF should be small, which is consistent with the calculation results (see Figure 8a,c).
The details above show that through the harmonic analysis of the measured data of shaft roundness, the UAFF can be reasonably predicted.

Effect of the Cylindricity
Source-B, Source-C and Source-D are used to construct the cylinder and calculate the UAFF according to the parameter settings in Table 1. Figures 11-13 are the calculation results of different cylinders. For ease of comparison, the three subgraphs of each cylinder adopt the same abscissa and the same ordinate. The UAFF on conical shafts constructed by different data sources and cylindricity are relatively close, and the force is smaller as the speed increases, as shown in Figure 11. As can be seen from Figure 12, for a concave shaft, the UAFF on the shaft is directly proportional to the cylindricity deviation. When the cylindricity deviation is less than 2 µm, the UAFF at all speeds is not more than 2 N. The UAFF on the convex shaft is inversely proportional to the cylindricity error value, and the UAFF on this shaft is the smallest of the three, as shown in Figure 13. In general, the UAFF on the three shafts decreases with the increase in rotating speed. When the cylindricity deviation is less than 2 µm, the resulting shaft eccentricity shall not exceed 0.05 µm.

Conclusions
When the spindle is running, the rotation accuracy is affected by the coupling of manufacturing errors, assembly faults, unbalanced magnetic pull in the motor, thermal deformation and other factors. It is difficult to determine the influence of a certain factor on the final rotation accuracy. Based on the setting and calculation of this paper, the following conclusions can be drawn:

1.
The shape error of the cylindrical surface of the shaft will make it subject to UAFF during rotation. In general, the smaller the error, the smaller the UAFF. When the roundness deviation is ≤1.25 µm, and the cylindricity deviation is ≤2 µm, the influence of the shape error on the rotation accuracy of shaft is small enough.

2.
The CIF caused by the shape error is very small (<1 ×10 −5 N), and the influence of this force on the rotation accuracy can be ignored. The magnitude of the CIF is determined by the shaft mass, rotation speed and the eccentricity of the center of mass A i,j P 2 i+1,j + B i,j P 2 i−1,j + C i,j P 2 i,j+1 + D i,j P 2 i,j−1 − E i,j P 2 i,j = F i,j P i,j + G i,j P i+1,j − P i−1,j . (A7) Set Equation (A7) is further simplified as E i,j P 2 i,j + F i,j P i,j − I i,j = 0. (A9) By using the solution of one variable quadratic equation, the following equation can be obtained: