Mathematic Modelling and Numerical Analysis for a Novel Inner-Type Nutation Magnetic Drive

This paper presents an inner-type of magnetic gears for non-contact driving, based on the nutation gearing principle. An improved three-dimensional analytical model, based on the magnetic vector potential of rectangular permanent magnets, via the equivalent current, is deduced for the proposed nutation magnetic gears. The transformation of coordinates is applied for a global field solution. The analytical model and the finite element model have analyzed the output torque of the nutation magnetic drive, respectively. The results show that the values calculated by the two models are consistent. The proposed analytical model can provide a reference for the design and optimization of nutation magnetic gears and other permanent magnetic mechanisms with rectangular magnets.


Introduction
With the rapid development of industrial robots, joint reducers have been widely used in order to improve and ensure control accuracy. Robot joint reducers require characteristics such as a short transmission chain, small size, high power, light-weight and easy control. Cycloid and harmonic reducers are the two kinds of reducers that are widely used in robot joints. Meeting the same requirements, nutation reducers have attracted more and more attention recently due to its simple structure [1].
Nutation mechanical reducers were developed. Over the years, nutation mechanical gears have been widely applied in many fields, such as micro-robot reducer motors, propeller helicopter reducers, and robot wrist joint reducers [2][3][4]. Research on the nutation mechanical gear has always been studied, based on the principles of dynamics and kinematics [5,6]. However, traditional nutation drive, via the mechanical gear, generally faces some inherent issues such as friction loss, mechanical fatigue and vibration.
In order to overcome the problems caused by mechanical contact, non-contact drive has developed by applying permanent magnets. Based on the principle of magnetic field modulation, the torque density of magnetic gears can be improved [7][8][9][10]. Magnetic gears with high torque density and a high gearing ratio could be achieved to innovate the topologies. Jorgensen et al. [11] proposed a cycloidal magnetic gear, which could achieve a peak volumetric torque density of 183 kN·m/m 3 at a gear ratio of −21:1. Rens et al. [12] demonstrated a practical dual-stage implementation of the magnetic harmonic gear, which exhibits a −360:1 gear ratio and a torque density of up to 75 kN·m/m 3 . Chau et al. [13] proposed a new magnetic variable gear, which offered different gears ratios. Molokanov et al. [14] investigated a magnetic gear, which had a planetary structure with two internal engagements.
A number of mathematical models and analytical methods have been proposed to analyze the magnetic field characteristics of non-contact gears. Three-dimensional field solutions, based on the magnetic vector potential and magnetic scalar potential, were proposed for radially polarized magnetic vector potential and magnetic scalar potential, were proposed for radially polarized cylinders [15,16]. Further, analytical solutions have been extended in various magnetic gears, such as cycloid magnetic gears [11], planetary magnetic gears [17] and spur magnetic gears with a threesegment Halbach [18].
In consideration of the advantages of the magnetic gears mentioned above, the performance of traditional nutation mechanical reducers could be improved. Yao et al. [19,20] proposed a two-stage nutation magnetic drive with a transmission ratio of 105:1, and made a prototype to verify the feasibility of the proposed non-contact nutation magnetic drive. The magnetic field characteristics of the nutation magnetic drive were solved by commercial finite element software, and the magnetic field distribution and torque were analyzed. Due to the special structure, the analytical solutions for traditional magnetic gears are unavailable for nutation magnetic gears.
In this paper, an analytical solution for nutation magnetic drive is achieved. A novel inner-type of magnetic gears is first proposed. Then, the analytical solution of traditional magnetic gears is improved for the new nutation structure. The analytical expressions for a single magnet and the total magnets in three-dimensional space are deduced, respectively. The results are also compared with those solved by the finite element method.

Nutation Magnetic Drive Principle
The movement principle of the nutation magnetic device is similar to the mechanical one. In this paper, an inner-type nutation magnetic drive is proposed. This sketch of the nutation magnetic drive is illustrated in Figure 1. The main parts of the nutation drive are composed of a pair of permanent magnetic gears. One magnetic gear is for the nutation labeled MG1, and the other magnetic gear is for the output labeled MG2. The magnetic gears mesh in the axial direction. The nutation sleeve is installed on the input shaft, and its axis has an angle of φ with the axis of the input shaft called the nutation angle. When the input shaft rotates, the MG1 only wobbles in the nutation sleeve and its rotation is confined by the guide pin. The gear meshing of the MG1 with the MG2 can achieve the nutation drive, which is non-contacting. Then, the magnetic torque of the MG2, driven by the MG1, is transferred to the output shaft, which is fixedly connected to MG2. The sketch of the magnetic forces between the MG1 and MG2 is shown in Figure 2. The magnets are polarized in parallel. 'N' and 'S' present the polarities of the meshing faces, which are north-polar and south-polar, respectively. Each magnetic pole on the MG2 is subjected to the push and pull force from the adjacent magnetic poles on the MG1. The MG2 will rotate when the total magnetic force is strong enough, and it corresponds to the total torque. There is angle difference α between the magnetic poles with different polarities on MG1 and MG2 when the gears are meshed, called the torque angle. The sketch of the magnetic forces between the MG1 and MG2 is shown in Figure 2. The magnets are polarized in parallel. 'N' and 'S' present the polarities of the meshing faces, which are north-polar and south-polar, respectively. Each magnetic pole on the MG2 is subjected to the push and pull force from the adjacent magnetic poles on the MG1. The MG2 will rotate when the total magnetic force is strong enough, and it corresponds to the total torque. There is angle difference α between the magnetic poles with different polarities on MG1 and MG2 when the gears are meshed, called the torque angle.
The nutation magnetic drive can achieve a high reduction ratio. The movement process of magnetic gears, with the input shaft, rotated a circle, as shown in Figure 3. The MG1 outside is the Energies 2020, 13, 1346 3 of 14 frustum of a cone shape, which is placed eccentrically relative to the MG2 with a cylindrical shape. The pole numbers of the MG1 and MG2 are 32 and 30, respectively. The clockwise movement of the nutation will result in a small change in the inner rotor rotation in an anticlockwise direction. In the first position, the reference magnetic pole 1 on the MG1 meshes with the reference magnetic pole 2 on the MG2. At that point, the center of the MG1 is on the right side of the MG2. From positions 1 to 9, the input shaft rotates a circle with a step of 45 • . Reference pole 1 is still in its previous position, while reference pole 2 moves two pole positions. It can be found that the MG1 has rotated one clockwise revolution, while the MG2 has only rotated 1/15 of a revolution anticlockwise. The reduction ratio is much higher than that of the cylindrical magnetic gears. According to the above principle, the reduction ratio R of the nutation magnetic drive can be expressed as where p 1 and p 2 are the total pole numbers of the MG1 and MG2, respectively. The reduction ratio of the inner-type nutation magnetic gear, shown in Figure 3, is −15. The nutation magnetic drive can achieve a high reduction ratio. The movement process of magnetic gears, with the input shaft, rotated a circle, as shown in Figure 3. The MG1 outside is the frustum of a cone shape, which is placed eccentrically relative to the MG2 with a cylindrical shape. The pole numbers of the MG1 and MG2 are 32 and 30, respectively. The clockwise movement of the nutation will result in a small change in the inner rotor rotation in an anticlockwise direction. In the first position, the reference magnetic pole 1 on the MG1 meshes with the reference magnetic pole 2 on the MG2. At that point, the center of the MG1 is on the right side of the MG2. From positions 1 to 9, the input shaft rotates a circle with a step of 45°. Reference pole 1 is still in its previous position, while reference pole 2 moves two pole positions. It can be found that the MG1 has rotated one clockwise revolution, while the MG2 has only rotated 1/15 of a revolution anticlockwise. The reduction ratio is much higher than that of the cylindrical magnetic gears. According to the above principle, the reduction ratio R of the nutation magnetic drive can be expressed as where p1 and p2 are the total pole numbers of the MG1 and MG2, respectively. The reduction ratio of the inner-type nutation magnetic gear, shown in Figure 3, is −15.

Analytical Model of Nutation Magnetic Gears
In the electrostatic field solving problem, the permanent magnet can be treated as the equivalent volume current density Jm and surface current density jm, according to the equivalent current model of the permanent magnet [15,16]. If the permanent magnets are high-performance permanent magnets, such as the NdFeB material, the permeability is approximately constant in the volume. The

Analytical Model of Nutation Magnetic Gears
In the electrostatic field solving problem, the permanent magnet can be treated as the equivalent volume current density J m and surface current density j m , according to the equivalent current model of the permanent magnet [15,16]. If the permanent magnets are high-performance permanent magnets, such as the NdFeB material, the permeability is approximately constant in the volume. The permanent magnet is deemed to be a homogeneous magnetizing medium, so the volume current density J m is 0.
The torque of permanent magnet gears is generated by the Lorentz force applied to the equivalent surface current in an external magnetic field. The output torque can be expressed as where T is output torque, R is radial vector, j m2 is equivalent surface current of MG2, B ext is external magnetic flux density, and S is surface area of the magnets. In this section, the vector magnetic potential theory is first used to obtain a magnetic field solution for the source magnet MG1. Then, the equivalent current density distribution of MG2 is applied to obtain the output torque by Equation (2).

Magnetic Flux Density by a Single Pole of MG1
A single pole of MG1 in the (s p -u p, v p, w p ) coordinate system is shown in Figure 4. The subscript p presents the pole number, which is from 1 to p 1 . It is assumed that the magnetization is strictly parallel to the u p direction. Thus, the magnetization intensity M(r ) can be written as where and r is the source point coordinate and the ± term takes the alternating polarity of adjacent poles into account. Then, the equivalent surface current j m can be written as where n represents the unit external normal vector of the magnet's surface. If the magnetization direction of MG1 is positive, the results of the equivalent surface currents j m1 are shown in Table 1.
The equivalent surface currents on the other two surfaces are 0.
Energies 2020, 13, x FOR PEER REVIEW 5 of 15 (4) where n represents the unit external normal vector of the magnet's surface. If the magnetization direction of MG1 is positive, the results of the equivalent surface currents jm1 are shown in Table 1.
The equivalent surface currents on the other two surfaces are 0.

Surface
Range n jm1  Top Magnetic vector potential is used to calculate the magnetic field. For the field point outside the magnetizing medium, the magnetic vector potential generated by the equivalent magnetizing current can be expressed as [15] where µ 0 is the permittivity in vacuum, A is the vector magnetic potential and r is field point coordinates.
Taking the results in Table 1 into account, Equation (5) can be rewritten as In the rectangular coordinate system, Equation (6) can be extended as In the Cartesian coordinate system, there is Its components can be expressed as Then, the magnetic flux density, by a single pole of MG1, can be obtained.

Magnetic Flux Density by Total Poles of MG1
The coordinate transformation is applied to realize the global field solution. Due to the existence of the nutation angle, the (s p -u p , v p , w p ) coordinate system should rotate ϕ angle around the line 'AB', Energies 2020, 13, 1346 6 of 14 as shown in Figure 4, to form a new (o p -x p , y p , z p ) coordinate system. A single pole in the (s p -u p , v p , w p ) coordinate system and (o p -x p , y p , z p ) coordinate system is shown in Figure 5. In order to consider the total poles of MG1, the (o p -x p , y p , z p ) coordinate system should rotate θ p angle around axis z p to form a new (o-x, y, z) coordinate system, as shown in Figure 6. The θ p angle is defined as where p = 1, 2, . . . , p 1 .
Then, the magnetic flux density, by a single pole of MG1, can be obtained.

Magnetic Flux Density by Total Poles of MG1
The coordinate transformation is applied to realize the global field solution. Due to the existence of the nutation angle, the (sp-up, vp, wp) coordinate system should rotate φ angle around the line 'AB', as shown in Figure 4, to form a new (op-xp, yp, zp) coordinate system. A single pole in the (sp-up, vp, wp) coordinate system and (op-xp, yp, zp) coordinate system is shown in Figure 5. In order to consider the total poles of MG1, the (op-xp, yp, zp) coordinate system should rotate θp angle around axis zp to form a new (o-x, y, z) coordinate system, as shown in Figure 6. The θp angle is defined as where p = 1, 2, …, p1.  Considering the (sp-up, vp, wp) coordinate system and (o-x, y, z) coordinate system as shown in Figure 5 and Figure 6, it can obtain   Considering the (s p -u p , v p , w p ) coordinate system and (o-x, y, z) coordinate system as shown in Figures 5 and 6, it can obtain where Then, the magnetic flux density generated by a single pole in the (o-x, y, z) coordinate system can be expressed as where Therefore, the magnetic flux density generated by the total poles of MG1 in the (o-x, y, z) coordinate system can be expressed as where the (−1) p+1 term takes the alternating polarity of adjacent poles into account. In what follows, it is considered to be a source of an external field.

External Magnetic Flux Density of MG2
The B ext on MG2 is generated by MG1. The MG1 is in the (o-x, y, z) coordinate system, while the MG2 is in the (o -x , y , z ) coordinate system, as shown in Figure 7. The relationship can be expressed as

Torque Expression
The (o'-x', y', z') rectangular coordinate system should be transferred into the (o'-r', θ', z') cylindrical coordinate system to calculate the output torque. Only the radial component of Bext is needed, which can be written as Only the currents on the left and right surfaces are needed in order to calculate the output torque. If the magnetization direction of MG2 is positive, the results of the equivalent surface current jm2 are shown in Table 2. Taking the slip angle α into account, the output torque by the total poles of MG2 can be expressed as  Here, x 1 is the inner radius of MG1, x 2 is outer radius of MG1, and g is the minimum gap between MG1 and MG2.
Therefore, the external magnetic flux density on MG2, generated by the total poles of MG1 in the (o -x , y , z ) coordinate system, can be expressed as

Torque Expression
The (o -x , y , z ) rectangular coordinate system should be transferred into the (o -r , θ , z ) cylindrical coordinate system to calculate the output torque. Only the radial component of B ext is needed, which can be written as Only the currents on the left and right surfaces are needed in order to calculate the output torque. If the magnetization direction of MG2 is positive, the results of the equivalent surface current j m2 are Energies 2020, 13, 1346 9 of 14 shown in Table 2. Taking the slip angle α into account, the output torque by the total poles of MG2 can be expressed as (36) Left

Application of Green Formula
In order to solve Equations (11)-(13), the Green formula is introduced to calculate the magnetic flux density, which is given as In the rectangular coordinate (s p -u p , v p , w p ), it can be written as It can obtain where Energies 2020, 13, 1346 10 of 14 Then, Equations (11)-(13) can be organized as

Application of Simpson Method
The Simpson method is applied for further calculation of the integral equations [16]. The integral line is uniformly divided. The meshing is defined by where a 1 is the lower limit of the integration, a 2 is the upper limit of the integration, and i m is the total number of grids along the integral line. Additionally, the integration coefficients of the Simpson formula are defined by Applying Equation (49) along the three axes of the (s p -u p , v p , w p ) coordinate system, the a(i) becomes u(q), v(m) and w(n), respectively. The S a (i) becomes S u (q), S v (m) and S w (n), respectively. Then, the Equations (46)-(48) can be written as Energies 2020, 13, 1346 11 of 14 In the same way, applying Eq. (44) along the three axes of the (o -x , y , z ) coordinate system, the a(i) becomes x (s) and z (p), respectively. The S a (i) becomes S x (s) and S z (p), respectively. Then, Equation (35) can be written as Where

Model Verification
The finite element method (FEM) is one of the effective methods for electromagnetic calculation [22]. In order to verify the accuracy of the three-dimensional analytical model, a commercial finite element software ANSYS Maxwell was used to calculate the output torque for comparison. The solution type of the magneto static was selected to solve the problem. The computational domain was set to 800 × 800 × 250 mm 3 , which was about 10 times that of the magnetic gears. The material of the magnetic gears was set to NdFeB35, and the material of the space was set as a vacuum. In order to set the magnetization direction, the sub-coordinate system was used for every magnet. Default boundary conditions were used. The boundary conditions of the object interfaces were natural boundaries, and the outer boundaries were Neumann boundaries. The adaptive mesh refinement was used, with the limitation of the maximum length of elements being 2 mm for all magnets. The maximum energy error was set as 0.5%. As a final step, an example of the mesh distribution is shown in Figure 8. It can be found that the sizes of the mesh grids on the magnets were very small after the computational convergence. All the parameter values were used for finite element analysis and analytical solutions are given in Table 3.
Energies 2020, 13, x FOR PEER REVIEW The finite element method (FEM) is one of the effective methods for electromagnetic calculation [22]. In order to verify the accuracy of the three-dimensional analytical model, a commercial finite element software ANSYS Maxwell was used to calculate the output torque for comparison. The solution type of the magneto static was selected to solve the problem. The computational domain was set to 800 × 800 × 250 mm 3 , which was about 10 times that of the magnetic gears. The material of the magnetic gears was set to NdFeB35, and the material of the space was set as a vacuum. In order to set the magnetization direction, the sub-coordinate system was used for every magnet. Default boundary conditions were used. The boundary conditions of the object interfaces were natural boundaries, and the outer boundaries were Neumann boundaries. The adaptive mesh refinement was used, with the limitation of the maximum length of elements being 2 mm for all magnets. The maximum energy error was set as 0.5%. As a final step, an example of the mesh distribution is shown in Figure 8. It can be found that the sizes of the mesh grids on the magnets were very small after the computational convergence. All the parameter values were used for finite element analysis and analytical solutions are given in Table 3.  The calculated output torques, with the torque-angle by the analytical model and the finite element model, are compared in Figure 9. Good agreement can be seen between the two calculation methods. The relative errors are all within 4%. The maximum torque of the gear is 18.75 N·m, which occurs at a torque angle of 6 • .  The calculated output torques, with the torque-angle by the analytical model and the finite element model, are compared in Figure 9. Good agreement can be seen between the two calculation methods. The relative errors are all within 4%. The maximum torque of the gear is 18.75 N•m, which occurs at a torque angle of 6°.  The numerical results of the output torque with different nutation angles under the torque angle of 6 • are shown in Figure 10. Good agreement can also be seen between the two calculation methods. When the nutation angle is increased from 1 • to 10 • , the drive torque is decreased from 19.44 N·m to 14.93 N·m. When the nutation angle is 0 • , the nutation magnetic gear is akin to a one-wave harmonic magnetic gear [10]. It can be found that the maximum torque is 19.44 N·m. Accordingly, the maximum torque density is 128.78 kN·m/m 3 . The numerical results of the output torque with different nutation angles under the torque angle of 6° are shown in Figure 10. Good agreement can also be seen between the two calculation methods. When the nutation angle is increased from 1° to 10°, the drive torque is decreased from 19.44 N•m to 14.93 N•m. When the nutation angle is 0°, the nutation magnetic gear is akin to a one-wave harmonic magnetic gear [10]. It can be found that the maximum torque is 19.44 N•m. Accordingly, the maximum torque density is 128.78 kN•m/m 3 . Compared to the FEM, this analytical method has many advantages, such as a clear concept, a clear physical model and fast calculation speed. However, there are always assumptions in the calculation, which may cause small deviations, as shown in Figure 9 and Figure 10. One cause may Compared to the FEM, this analytical method has many advantages, such as a clear concept, a clear physical model and fast calculation speed. However, there are always assumptions in the calculation, which may cause small deviations, as shown in Figures 9 and 10. One cause may be the assumption of a unity-relative permeability of the permanent magnets in the derivation of the analytical model. The finite element model used the relative permeability of 1.09978 for the NdFeB35 magnets. A disadvantage of the proposed analytical method is limited to regular geometries of magnets and arrangement.

Conclusions
A novel inner-type of nutation magnetic gears with rectangular permanent magnets was proposed. A three-dimensional analytical model was deduced on the basis of the magnetic vector potential and equivalent current of rectangular permanent magnets. The analytical expression of the out torque is expressed by Equation (54). The calculated output torques by the analytical model were compared to those of the finite element model, and the results by the two methods are consistent. It was shown that, at gear ratios of −15:1 and small gear radius of 43 mm, a torque density of up to 128.78 kN·m/m 3 can be achieved. The analytical method can also be applied for other permanent magnetic mechanisms with rectangular magnets.