A Study on Characteristics of Flow Field of Air Gap in Maglev Spherical Active Joints for Multi-DOF Spindle System

Fan ZHANG*, Li ZENG**, Yan LI***, Zi ZENG****, Zhida ZHU*****, Jin SUN******, Min DAI******* *School of Mechanical Engineering, Yangzhou University, Yangzhou225127, China, E-mail: fzhang@yzu.edu.cn **School of Mechanical Engineering, Yangzhou University, Yangzhou225127, China, E-mail: lizengcf@163.com ***School of Mechanical Engineering, Yangzhou University, Yangzhou225127, China, E-mail: yli@yzu.edu.cn ****School of Mechanical Engineering, Yangzhou University, Yangzhou225127, China, E-mail: zzeng@yzu.edu.cn *****School of Mechanical Engineering, Yangzhou University, Yangzhou225127, China, E-mail: zdzhu@yzu.edu.cn ******School of Mechanical Engineering, Yangzhou University, Yangzhou225127, China, E-mail: jinsun@yzu.edu.cn *******School of Mechanical Engineering, Yangzhou University, Yangzhou225127, China, E-mail: mdai@yzu.edu.cn


Introduction
Spindle tools capable of six-degree-of-freedom spatial movement in any posture is one of the key features of existing advanced CNC machine tools.The multi-degreeof-freedom spindle systems are generally classified into two types: new spindle with parallel structure and traditional spindle with cascaded structure [1].The former changes its position and posture by employing several servo-motor or hydraulically driven links with controllable lengths, whereas the latter utilizes multiple mechanical rotating structures, each of which controls the rotation for only one degree-of-freedom, thus controlling the direction of spindle tool (Fig. 1).Both types of spindles are faced with the following defects: 1. numerous joints complicates the mechanical structure and worsens the driving efficiency and reliability; 2. serious abrasion caused on the friction surfaces of joints response speed and dynamic characteristics; and 3. accumulation of mechanical errors resulting in lowering the accuracy of control system.To address these defects, quasi-human joints with better efficiency could be constructed based on simplified spherical joints, that is to say to utilize maglev spherical active joints supported and driven by motors with spherical surface in contact [2][3][4][5][6][7].Using such kind of active joints in a CNC machine will reduce the complexity of multi-DOF spindles and improve their position and posture accuracies and response speed.Under high-speed operation of spindles, rotor surface drives the air in gap to rotate along with the rotor and form a circumferential air circulation in gap, thereby engendering an air film in the flow field of air gap between the stator and the rotor to generate a levitating force or resistance [8].For the film, how its flow field pressure and dynamic characteristics are distributed and changed exerts direct influence upon the dynamic and static performance of joint.However, present studies normally ignore the effects of the levitating force or resistance generated by the air film upon rotor's movement and dynamic performance, which explains the deviation between theoretical analysis and reality in terms of their motions.A lot of scholars [9][10][11][12][13][14][15][16] have achieved some valuable results in both terms of theoretical modeling and experimental verification, which have been practiced in engineering.However, studies on the characteristics of flow field of air gap in maglev spherical active joint have been little reported in both China and other countries.This paper explored the flow field of air gap in maglev spherical active joint and its dynamic & static characteristics based on aerodynamics, established mathematical models for dynamic & static pressures and kinetic parameters of flow field of air gap, and, with simulation and experiment, investigated how the flow field of air gap impacts its static load-carrying capacity and dynamic characteristics.

Principles of joint generating electromagnetic torque and levitating force
However, the difference from conventional motors is that, the joints should extend its rotor from uni-dimensional rotation around a fixed axis to spatial rotation around a fixed point.In this case, rotor's angular velocity vector ω has to be changed from a directed scalar to a spatial vector in an arbitrary direction.ω = ωxi + ωyj+ ωzk, where i, j and k represent unit vectors in three arbitrary axes, respectively.However, nonlinear coupling relations are shown among mechanical equilibrium equations of three orthogonal axes.Thus, to effectively and separately manipulate the three vector components of angular velocity for the rotor of joint, three axial electromagnetic torques Mx, My and Mz, which are independently controllable.In an orthogonal coordinate system, the three torque vector components are combined as M = Mxi + Myj+ Mzk, and the resultant torque is a spatial vector in an arbitrary direction.To yield electromagnetic torques in arbitrary controllable spatial directions, the maglev spherical active joint should produce a resultant magnetic field, whose amplitude and axis direction can totally be manipulated.Based on a Cartesian coordinate system, this paper prepared threephase AC windings around each coordinate axis.A two-dimensional magnetic field vector was formed by electrifying symmetrical three-phase windings with symmetrical threephase alternating currents, and the magnetic field vector rotated round its axis.These magnetic field vector could also be formed by two phase AC windings with two quadrature alternating currents.Then combining 2D magnetic field vectors of all 3 coordinate axes with each other formed a 3D magnetic field, which rotated around an arbitrary axis passing through a fixed spatial point.As shown in Fig. 2, three windings that generated 2D magnetic fields are represented by Wx, Wy and Wz, respectively.Therefore, the structure of a maglev spherical active joint can be designed according to the formation mechanisms of a spatial rotating magnetic field, as shown in Fig. 3.The core of stator also provided a spherical surface internally to accommodate to the spherical rotor and engender a uniform air gap between them.
Three windings (corresponding to three coordinate axes) with axes mutually orthogonal were also embedded into the slotted core of stator.Each winding was divided into two sections, which were symmetrically wrapped around the core of stator in two hemispherical slots.Fig. 4 shows the electromagnetic winding in z direction of the rotor of maglev spherical active joint, where Wz is divided into two symmetrical windings z W + and z W − , representing the posi- tive and negative hemispherical windings of stator, respectively, with 1 and 2 referring to the spherical rotor and stator of joint, respectively.Windings wrapped around 3 axes were energized by adjustable 3-phase alternating currents iA, iB and iC and a controller was employed to control torques and speeds (or angular displacements) for both winding sections on each axis, thereby generating parallel electromagnetic torques with same direction to drive the rotor to rotate around this axis.For example, providing the generated resultant electromagnetic torque in z-axis was =+.The levitated rotor shares the same work- ing principles with those traditional electromagnet bearings.As illustrated in Fig. 4, when the rotor in z direction deviated from its balanced position under interference, its gap in positive z-axis shrank.For such conditions, to return the rotor back to stator's centre (that's the balanced position) under closed-loop control, the currents applied onto the hemispherical windings of joint's stator below the rotor were increased by z i  , whereas the currents applied onto other side were reduced by z i  .The difference between electromagnet pulls generated by upper and lower hemispherical windings of stator enabled the rotor to return to the balanced position at the centre of stator and the resultant levitating force was z z z =−, as shown in Fig. 4. The rotation and levitation mechanisms in x-axis and y-axis can be analysed in the same way.When the rotor of maglev spherical active joint rotates at a high speed around the spindle, the air circulation in the flow field of air gap can directly affect the dynamic & static performances of joint.The characteristics of air circulation are in close relation with the variation pattern of air gap length.As shown in Fig. 5, a, a 3D coordinate system was established for the maglev spherical active joint.
The air gap is illustrated in Fig. 5, b when there is no deviation for the centre of sphere of rotor, where r is the outer radius of rotor; R is the inner radius of stator; and c is the average length of air gap when the centre of sphere of stator coincides with the centre of sphere of rotor.Fig. 5, c shows the case where the centre of sphere of rotor deviates from the balanced position, where h is the air gap length between stator and rotor at arbitrary spatial positions.

R isin cos jsin sin kcos
The eccentricity from the geometrical centre of rotor O1 to the geometrical centre of stator O is represented as: Take the geometrical centre of rotor O1 as the origin of a radius vector, which intersects with the internal surface of stator at point P.As the eccentricity e is far less than r and R,

Static load-carrying capacity model for flow field of air gap in maglev joint
3.1.Static pressure distribution model for air film in air gap of maglev joint.
When the rotor rotates at a high speed, the air in gap will form a layer of air film against the stator.The static pressure distribution on the air film can be derived from continuity equation of fluid, equation of motion, state equation of fluid and energy equation.A pressure distribution equation can also be derived for the compressive flow field in a spherical coordinate system and its dimensionless form is expressed as follows: ( ) where: Λ is the load-carrying coefficient, t is the time and ω is the circular frequency for rotor working.

Static load-carrying capacity of air film in air gap of maglev joint
As the static pressure distribution equation of air film has no relation with the time, each of following static variables is represented with a subscript 0 for the purpose of being distinguished.So according to Eq. ( 6), the static pressure distribution equation we can obtain is as follows: ( ) Then solve for the boundary conditions: 1) At the place where the air gap contacts with the atmosphere, P0=Pa; 2) The rotor has an eccentric angle We can combine Eq. ( 7) and Eq. ( 8) to get the thickness H0 and the pressure P0 of air film.Further, after we obtain the steady state, the load on the static air film and the damping moment of air film with the direction along the axis of rotation can be expressed as: where: ( )

Kinetic model for flow field of air gap in maglev spherical active joint
Due to effects from mass centre and external forces, the rotor rotating at a high speed could suffer turbulence at its balanced position.The dynamic response of air film affects the stability of rotor, so it is necessary to study its dynamic characteristics.Under turbulence, the rotor will deviate from its balanced point, where the three coordinate components of the turbulence are represented by Ex0, Ey0 and Ez0, respectively.Assuming the turbulence appears in a sinusoidal form, the rotor position at any time can be expressed as [18]: where: f is the dimensionless turbulence frequency, , f   = v is the turbulence frequency of rotor; i is the imaginary unit.
So we can obtain the thickness of air film for the rotor at any position: where: H0 is the thickness of static air film and Hd is the thickness of dynamic air film.
where: Hd0 is the amplitude of dynamic air film, a complex number.
Under slight turbulence, the pressure of dynamic air film can be expressed as: where: Pd is the pressure of dynamic air film and Pd0 is the pressure amplitude of dynamic air film, a complex number.Substituting Eq. ( 12) and Eq. ( 13) into Eq.( 6) gives the following dynamic equation of air film: This dynamic equation of air film (14) implies the turbulence components Ex0, Ey0 and Ez0.In order to obtain the dynamic stiffness and damping coefficient of the air film under dynamic pressure, the partial derivative method is used.Let:  Finding the derivative of Eq. ( 14) and Hd0 gives the equations for PEx, PEy, and PEz.Equations for PEx are: . We can combine with Eq. ( 15) to obtain PEx.Then, PEy and PEz can be obtained in the same way.

Dynamic characteristic parameters for flow field of air gap in maglev joint
By obtaining the pressure P0 of static air film and substituting it into the dynamic equation of air film, we can obtain the dynamic pressure distribution of air film.With PEx, PEy and PEz known, we can obtain the dimensionless dynamic stiffness and the dynamic damping coefficient, expressed as follows, respectively:

Simulation and experimental study for static & dynamic performances of flow field of air gap in maglev spherical active joint
Write a program via finite difference method based on above analysis results to calculate the static & dynamic performances of air film in the spherical joint.By investigating how the eccentric position and speed of the rotor affect the static pressure and damping moment of the air film, as well as how the turbulence frequency of rotor influences the dynamic stiffness and damping coefficient of the air film, we can obtain the parameters of spherical joint and its operation parameters, as listed in Table 1.
Table 1 Parameters of spherical joint  6 illustrates the distribution pattern of loaded pressure for the flow field of air gap in joint.The loaded pressure of the flow field of air gap would be in proportion to the load-carrying stiffness of the air gap if a controller is used to maintain the air gap.Therefore, analysing the variation pattern of load-carrying capacity for the air gap could manifest the variance tendency of the air film stiffness.Many factors, mainly the average length of air gap, the speed of rotor and the eccentricity, can impact the distribution pattern of load-carrying capacity or stiffness of air film in the flow field of air gap.Fig. 6 Loaded pressure distribution for the air film How the load-carrying capacity of air film varies with the average length of air gap is shown in Fig. 7, where the load-carrying capacity drops as the average length goes up.The reason is that the load-carrying capacity will exponentially drop as the average length goes up.
Fig. 8 shows the variation of the load-carrying capacity of air film with the angular velocity of rotor, where it can be found both the load-carrying capacity and the stiffness increase as the rotation speed increases, because the load-carrying coefficient of air film as the rotation speed increases.However, since the film is communicated to the atmosphere, the load-carrying capacity and the stiffness cannot increase endlessly but their variances tend to be gentle in the end.Fig. 7 Variation of load-carrying capacity of air film Fig. 8 Variation of load-carrying capacity of air film with angular velocity of rotor Fig. 9 exhibits the variation of the load-carrying capacity of air film with the eccentricity, where the load-carrying capacity significantly increases as the eccentricity goes up.The reason is that once the film thickness h at a position changes, the film pressure will alter correspondingly to maintain continuity.
Providing the rotation speed of joint is 100000 r/min and the turbulence displacement in positive x-axis is 1μm, three contour plots, as shown in Fig. 10, can be obtained to illustrate the dynamic pressure change of air film when the air gap in z-axis lower spherical surface of joint is set at 0.5 s, 2 s and 5 s, respectively.According to simulation results, the dynamic pressure zone of the air gap keeps changing with the rotation of rotor, and the film pressure under turbulence displacement is higher than that under static balance, which is because the turbulence displacement increases the eccentricity, enabling the film pressure to be increased as well.Based on the colour distribution in these plots, the pressure around z-axis approximately shows a sinusoidal (or cosine) distribution.Fig. 12 shows an experiment setup for controlling an inductive maglev spherical active joint studied in this paper, and this setup was consisted of a joint, an input channel for data collected in the experiment, an output channel for control, and a high-speed processor.
First, the experiment setup as shown in Fig. 12 was used to measure the resultant supporting force of joint's rotor.Then, subtract the measured electromagnetic levitation force from the resultant supporting force according to Newton's second law, and the obtained result was approximately equal to the force by which the air film in the flow field of air gap supported the rotor rotating at a high speed.we only compared the abovementioned theoretical simulation results with the experimental results in terms of variation pattern or tendency.When the rotor of the maglev spherical active joint rotated at a high speed around the z-axis, the air circulation pressure distribution pattern in the flow field of air gap, which was located in the lower hemispherical surface of joint's stator and vertical to the spindle and passed through the section with maximum load, is shown in Fig. 13.Its distribution pattern and the air film pressure simulation results are similar to those shown in Fig. 10.The air circulation pressure distribution pattern mainly was dependent on c r , the ratio between air gap length (clearance) and rotor radius.As shown in the figure, when the gap was small, the air film acted as a fluid plain bearing and the Reynolds number effect occurred.But when the ratio increased, the Bernoulli effect arose, which caused the rotor to vibrate and displace, leading to increased flow rate and decreased pressure at the side where the gap decreased.As the vibration displacement increased further, the rotor in vibrating would lose its stability, thus determining the boundary condition for the stability of vibration of rotor [19].

2.
c n   Where: ωc is the critical speed of the rotor losing stability; and ωn is the intrinsic frequency of rotor.14, a shows how the loadcarrying capacity of air film changes with the rotation speed in general, the capacity increases as the speed increases, with the variation tendency and the simulation analysis results as shown in Fig. 7.However, due to restricted practical test conditions and test errors, there was little consistency between test data and theoretical data.Fig. 14, b represents how the load-carrying capacity of air film changes with the air gap under different rotation speeds of the joint, where it can be found that, the shown tendency of the radial loadcarrying capacity changing with the air gap under different rotation speeds is basically consistent with the theoretical simulation analysis results shown in Fig. 7 i.e. the capacity drops as the air gap increases.Fig. 14, c shows how the load-carrying capacity of air film changes with the eccentricity under different rotation speeds of joint.It can be found in this figure that regardless of the speed of joint, all curves show basically the same variation pattern, which is generally consistent with the theoretical simulation analysis results.As eccentricity goes up, the load-carrying capacity air film increases and its hydrodynamic effect turns more significant, which can benefit the operation performance of joint.

Conclusions
Since a machine spindle requires multi-DOF movement, this paper proposes a multi-DOF maglev spherical active joint, and analyses the three-dimensional rotation of active joint rotor around any axis of a fixed spatial point, as well as the principle of maglev fields.It also established the mathematical models for static & dynamic pressures and dynamic parameters of air gap flow field in the maglev spherical active joint, based on Reynold's equation of dynamic pressure of flow field air.In addition, simulation and experimental study are performed on air film's load-carrying capacity in the air gap flow field, as well as its dynamic characteristics.The experimental and simulation results are qualitatively compared to validate the effectiveness of air gap flow field model, with discussions on the main factors affecting the joint's load-carrying capacity.The research shows that in the drive application of multi-DOF active joint, the influence of air gap flow field is ignored, but its effect grows more significant with the increase of spindle rotor speed.Therefore, the research in this paper provides a new idea and method to make better use of 6-DOF maglev spherical active joint in spindle research and development, and to explore the dynamic characteristics of high-speed spindle and the use of pneumatic structure.

Fig. 1
Fig. 1 Multi-DOF CNC machines: a -spindle with cascaded structure, b -spindle with parallel structure

2. 2 .Fig. 5 A
Fig. 5 A 3D spherical coordinate system for flow field of air gap in a maglev spherical active joint: a -3D coordinate system of active joint, b -the air gap without deviation of the ball center, c -the air gap with deviation of the ball center let P0 = Pa at the place where ( )

4. 1 .
Dynamic pressure distribution model for air circulation in air gap of maglev joint

1 .
At the place where the air gap contacts with the atmosphere,

1 .
Fig.6illustrates the distribution pattern of loaded pressure for the flow field of air gap in joint.The loaded pressure of the flow field of air gap would be in proportion to the load-carrying stiffness of the air gap if a controller is used to maintain the air gap.Therefore, analysing the variation pattern of load-carrying capacity for the air gap could manifest the variance tendency of the air film stiffness.Many factors, mainly the average length of air gap, the speed of rotor and the eccentricity, can impact the distribution pattern of load-carrying capacity or stiffness of air film in the flow field of air gap.

Fig. 9 Fig. 11
Fig. 9 Variation of load-carrying capacity of air film with eccentricity Fig. 11 shows the dynamic stiffness diagram and the dynamic damping coefficient diagram of air film when the eccentricity is set to εx = 0.4, εy = 0.4 and εz = -0.6.We can see from Fig. 11 that the dynamic stiffness and the damping coefficient show close relations with the turbulence frequency of rotor.Generally, when f < 2, as the turbulence turns more frequent, dynamic stiffness Kxx, Kyy and Kzz increase, while the absolute values of Kxy, Kxz, Kyx, Kyz, Kzx and Kzy decrease.All absolute values of dynamic

Fig. 12 Fig. 14 Factors
Fig.12 Control system for joint prototype