Multicoil torque modeling of reluctance spherical motor based on linear interpolation

Abstract

Due to the large amount of finite element calculation resources, the calculation stability is low. In order to obtain accurate torque data conveniently, a new torque modeling method is proposed in this paper. By establishing the minimum torque model, the finite element software is used to obtain elementary torque map, and the torque within the whole rotor motion range is built by polynomial fitting, and then the position torque in the process of rotor motion under different currents is obtained by multiple linear interpolation method, and converted to the global coordinate system by projection. The multicoil torque is superimposed and corrected by using the dynamic error coefficient to obtain the final synthetic torque. Finally, the effectiveness of the modeling method is verified by simulation and experiments. The results show that this method not only ensures the modeling accuracy, but also speeds up the calculation.

Keywords

Reluctance spherical motor (RSM)linear interpolation torque map multicoil dynamic error coefficient

1. Introduction

With the advancement of science and technology and the development of society, the manipulator and robot industries are gradually rising [1]. The traditional single-degree-of-freedom operation of motor no longer meets the needs of industry and life. Although the traditional motor can meet the multi-degree-of-freedom operation through multiple traditional motors, this situation will lead to a bulky system, high control requirements, mutual errors and other problems.

In order to better meet the needs of multi-degree-of-freedom operation, scholars worldwide began to study spherical motors for a long time, and put forward a variety of spherical motor structures from the earliest developed magnetoresistive stepping spherical motor to induction spherical motor, permanent magnet spherical motor, Halbach array spherical motor and magnetoresistive spherical motor [2,3]. Compared with traditional motors, these spherical motors can complete the operation of multiple degrees of freedom by one motor, which greatly simplifies the mechanical structure and reduces the manufacturing cost [4]. So far, more and more scholars have begun to study spherical motors.

Due to the strong coupling electromagnetic structure and good performance of spherical motor, the diversity of three-dimensional motion, the size and direction of electromagnetic force are different, and the time and position of each magnetic pole are different. Therefore, it is difficult to analyze the magnetic field and the torque of the motor. Generally, the finite element method is used for global analysis [5,6]. The development of permanent magnet motor is very early and fast, the reason is that there is no iron core material added to the spherical motor, which is linear. So far, the modeling methods of permanent magnet spherical motor include: fitting modeling of permanent magnet spherical motor with finite element software [7,8], three-dimensional magnetic field modeling based on spherical harmonic function and Laplace equation [9], modeling based on surface current model and Lorentz force method [10], etc.

The research of reluctance spherical motor involves the fields of electromotor, microelectronics and power electronics. In addition, its nonlinear characteristics increase the difficulty of research. However, with the progress of science and technology, scholars have also made great progress in the research of reluctance spherical motor. Lee and others designed a double excitation two degree of freedom reluctance spherical motor for surveillance camera [11], and analyzed the static torque characteristics by using the equivalent magnetic circuit method [12]. Zeng designed a spherical switched reluctance magnetic levitation motor based on the magnetic levitation theory and the principle of minimum reluctance, derived the magnetic permeability of the gap by using the magnetic segmentation method, and obtained the mathematical model of the spherical motor [13], through the established three-dimensional coordinate system, a rotor position detection method based on radial displacement sensor and angular displacement sensor is proposed [14]. In recent years, Ju proposed a three-degree-of-freedom spherical reluctance motor, established the basic mathematical model of the spherical motor according to the electromagnetic law and energy conversion principle, analyzed the motor by three-dimensional finite element method, and obtained the characteristics of magnetic field distribution and static torque during spin and pitch [15,16]. Later, Qiao based on the finite element software, the virtual displacement method was used to model the motor and establish the inductance curve and moment angle characteristics of the motor [17]. Recently, Xu proposed a sensorless rotor posture detection method based on the mutual inductance voltage of the stator coil to simplify the position detection element of a reluctance spherical motor. This method is based on the real-time voltage generated by stator coil mutual inductance, combined with motor structure constraints and intelligent algorithm to reverse the three-degree-of-freedom position angle [18].

Due to the complex magnetic field distribution in RSM and its nonlinear characteristics, it is very difficult to calculate by the analytical method, and the error after numerical differentiation is very large, so it is difficult to establish a correct model. Therefore, in order to quickly obtain an accurate torque model, modeling based on the finite element method is feasible. In this paper, the elementary model of torque is used for finite element calculation, which reduces the model calculation to the area corresponding to a single stator coil and rotor, which greatly shortens the calculation time, and can quickly calculate the torque of a single coil under the condition of withstanding current. Then, the torque in the range of rotor spin and pitch motion is built by polynomial fitting, and the torque value of each position in the process of rotor motion under single coil can be calculated by the interpolation method. The torque in the spin or pitch motion under the single coil is projected into the global coordinate system, and the dynamic error coefficient is calculated according to several important positions between the rotor and stator coil, which is fitted by binary polynomial. Finally, based on the dynamic error coefficient and multicoil torque, the synthetic torque is obtained.

2. Structure and coordinate system division

2.1. Spherical motor structure

The main parts of the motor include stator spherical shell, rotor body, coil, iron core and output shaft. The basic structure of the spherical motor is shown in Fig. 1. In order to facilitate multi-degree-of-freedom motion, the stator shell and rotor shell are designed as spherical. The rotor is divided into six poles, with an interval of 60° between two adjacent magnetic poles. The rotor is pressed by silicon steel sheet. The surface of the rotor is an arc surface in the shape of a salient pole. The stator core and housing are made of DT4 material. There are 24 energizing coils in total, which are centrally distributed and divided into three layers, with a difference of 33° for each layer. The middle layer coil is located on the equatorial plane, and the other two layers are evenly distributed on both sides of the equator, with an interval of 45° between adjacent coils on the same layer. A pole shoe is added at the end of the stator core, which is wider than the stator core. See Table 1 for the detailed data of the motor.

Fig. 1.

Structure diagram of the RSM.

Table 1

Main parameters of the RSM

Parameters	Values
Outer radius of stator shell	122 mm
Inner radius of stator shell	75 mm
Width of stator pole	10 mm
Length of stator pole	10 mm
Height of stator pole	15.5 mm
Outer radius of rotor	73 mm
Inner radius of rotor	30 mm
Width of rotor pole	14 mm
Height of rotor pole	12.5 mm
Air gap	1 mm
Coil turns	350 turns

The operation principle of the motor is similar to that of the traditional switched reluctance motor. The motor can rotate around the Z-axis by changing the power of the adjacent coils on the same layer; When the adjacent coils of the upper and lower layers are energized, the pitching motion around the X-axis and Y -axis can be realized. The tilt motion range of reluctance spherical motor is merely between ±33°, and the spin motion range is 0 ∼ 360°.

2.2. Global and local coordinate systems

Firstly, the coordinate system of the spherical motor is established. The spherical motor rotates with three-degrees-of-freedom, and the change is complex. In order to describe the operation and position relationship of the motor conveniently, the coordinate system is established. For the convenience of finite element calculation, the local coordinate system and global coordinate system are established based on the stator.

Fig. 2.

Stator coil marking diagram: (a) three layers; (b) middle layer.

Here, M is used to define the coils in the middle layer, B is used to define the coils in the lower layer, and T is used to define the coils in the upper layer, as shown in Fig. 2. Firstly, the local coordinate system of M1 is established, which takes the spherical center as the coordinate origin. The central direction of the stator core passing through M1 is the X-axis, the vertical upward is the Z-axis, and the Y -axis direction points to the center of the M3 stator core. Based on the local coordinate system of M1, the local coordinate system of M2 is obtained rotating 45° counterclockwise. Take the center direction of the stator core passing through M2 is the X-axis, and the vertical upward is the Z-axis, then the Y -axis direction points to the center of the M4 stator core. The local coordinate system of M3 can be obtained by rotating 90° counterclockwise on the basis of the local coordinate system of M1. Take the center direction of the stator core passing through M3 is the X-axis and the vertical upward is the Z-axis, and the Y -axis direction points to the center of the M4 stator core. Similarly, the local coordinate systems of M4, M5, M6, M7 and M8 can be established. For the uppermost local coordinate system, based on the M1 local coordinate system, the T1 local coordinate system can be obtained by pitching up 33° in the Y -axis. Based on the M1 local coordinate system, first rotate 45° counterclockwise, and then pitching up 33° in the Y -axis to obtain the T2 local coordinate system. By analogy, the local coordinate systems of T3 to T8 can be obtained. For the bottom local coordinate system, based on the M1 local coordinate system, the local coordinate system of B1 can be obtained by pitching down the Y -axis by 33°. The local coordinate system of B2 can be obtained by selecting 45° counterclockwise and then pitching down the Y -axis by 33° based on the M1 local coordinate system; By analogy, the local coordinate systems of B3 to B8 can be obtained. In order to simplify the calculation, here we take the local coordinate system of M1 as the global coordinate system. Use O − X _ji Y _ji Z _ji to represent the local coordinate system and O − XYZ to represent the global coordinate system, where $j=\{M,T,B\}$ , $i=\{1,2,\ldots ,8\}$ . The marking of local coordinate system is shown in Fig. 3.

Fig. 3.

Marking diagram of the local coordinate system.

Fig. 4.

Definition diagram of each included angle.

The included angles between the local coordinate system and the global coordinate system are defined as ∠Xx, ∠Xy, ∠Xz, ∠Yx, ∠Yy, ∠Yz, ∠Zx, ∠Zy, ∠Zz and the subdivision of the included angle is shown in Fig. 4. The upper case letter represents the name of each coordinate axis in the global coordinate system, and the lower case letter represents the name of each coordinate axis in the local coordinate system. Since the positions of the local coordinate system and the global coordinate system are determined and do not change, all included angles of each local coordinate system and the global coordinate system are determined, and 207 included angle values can be obtained directly through measurement.

It is necessary to verify the torque relationship between the local coordinate system and the global coordinate system. By applying 1 A current to M2 coil, it is calculated that under local coordinate system M2 and local coordinate system M1, the rotor teeth are at the initial position of zero degree, rotating and pitching from 0 to 45° at the same time. The nine included angle values between the two relative coordinate systems are calculated, and the torque value under local coordinate system M2 is converted to the torque value under local coordinate system M1 by transformation, as shown in Eq. (1). Comparing the calculated values with the simulation data under the local coordinate system M1, it is found that the two values are the same, as shown in Fig. 5. Therefore, it shows that the values of a torque vector in two different coordinate systems can be converted to each other; that is, the torque in any local coordinate system can also be projected into the global coordinate system. $\begin{eqnarray}\displaystyle \left[\begin{array}{@{}c@{}}T_{x1}\\ T_{y1}\\ T_{z1}\end{array}\right]=\left[\begin{array}{@{}ccc@{}}\cos \angle 45^{\circ } & \cos \angle 135^{\circ } & \cos \angle 90^{\circ }\\ \cos \angle 315^{\circ } & \cos \angle 45^{\circ } & \cos \angle 90^{\circ }\\ \cos \angle 90^{\circ } & \cos \angle 90^{\circ } & \cos \angle 0^{\circ }\end{array}\right]\left[\begin{array}{@{}c@{}}T_{x2}\\ T_{y2}\\ T_{z2}\end{array}\right]. & & \displaystyle\end{eqnarray}$ (1)

3. Single coil torque map and linear interpolation

3.1. Torque mathematical model and stator region

For the reluctance spherical motor, its nonlinearity greatly aggravates the complexity of the motor. Because there are ferromagnetic materials in of reluctance spherical motor, the magnetic circuit has saturation. In this paper, the virtual displacement method is preferred to analyze the electromagnetic torque between a single stator coil and rotor [17], as shown in Eq. (2). $\begin{eqnarray}\displaystyle dW_{e}=dW_{m}+dW_{mec} & & \displaystyle\end{eqnarray}$ (2) where $\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}\displaystyle dW_{e}=-eidt=id{\psi}=i\frac{\partial {\psi}(i,x)}{\partial x}dx+i\frac{\partial {\psi}(i,x)}{\partial i}di\\[10.0pt] \displaystyle dW_{m}=\frac{\partial W_{m}(i,x)}{\partial x}dx+\frac{\partial W_{m}(i,x)}{\partial i}di\\[10.0pt] \displaystyle dW_{mec}=f_{em}dx,\end{array}\right. & & \displaystyle\end{eqnarray}$ (3) where dW _e is the net electric energy input into the motor within dt time after removing the resistance loss, dW _m is the total energy absorbed by the magnetic field in dt time, and dW _mec expresses the increased mechanical energy in dt time. The spherical motor rotates at a small angle θ, as shown in Eq. (4). $\begin{eqnarray}\displaystyle T_{e}=\left.\frac{\partial W^{\prime }}{\partial {\theta}}\right|_{i_{k=\mathit{const}}}=\left.\frac{\partial W}{\partial {\theta}}\right|_{{\psi}_{k=\mathit{const}}}, & & \displaystyle\end{eqnarray}$ (4) where winding magnetic energy is $W^{\prime }=\int _{0}^{{\psi}}id{\psi}$ and winding magnetic energy is $W=\int _{0}^{i}{\psi}di$ .

Fig. 5.

Comparison of torque conversion.

Now, the finite element method is used to calculate the torque data of the spherical motor. ANSYS software is used to draw the diagram of the reluctance spherical motor and mesh it. In ANSYS, for the three-dimensional magnetic field analysis, the solution domain surrounding the whole region is established based on the air environment (vacuum), which is similar to the balloon boundary in the two-dimensional magnetic field solution. Therefore, the magnetic field boundary conditions required by RSM automatically meet the infinite field boundary conditions without additional setting. t is only necessary to set all coils of the spherical motor as insulation boundary conditions. For the grid division setting, the maximum unit side lengths of coil, rotor and stator are set to 5 mm, 8 mm and 10 mm respectively, and the maximum number of division units is limited to no more than 2000, get the grid situation in Fig. 6.

Fig. 6.

Mesh generation.

It is unrealistic to simulate the global torque of the reluctance spherical motor. Its calculation is huge and time-consuming and it also has high requirements for computers. Therefore, the method of sampling minimum model is more reasonable.

For the local coordinate system established in the previous section, the first step is to determine the area of each stator coil. Taking the coil with zero longitude and latitude as an example, the center of the coil core is taken as the regional center point, which occupies 22.5° of east longitude and west longitude respectively, and 16.5° of north latitude and south latitude respectively. The area occupies 45° longitude and 33° latitude. Therefore, this region is regarded as the elementary model domain for torque calculation.

Through the analysis of the motor structure, it is found that the shape of the coil area on the same layer of the spherical motor is the same. Through the simulation verification, it is concluded that the torque value in its local coordinate system is the same when a single coil on the same layer is energized. Therefore, as long as the torque data in the whole region under its local coordinate system can be obtained when a single coil with the same longitude on each layer of the three-layer coil is energized, the torque data in other coil coordinate systems on the same layer can be obtained directly by “analogy”. In order to simplify the amount of simulation data, through the judgment of the motor body and the verification of simulation experiments, it is found that the area and data of the upper coil and the lower coil of the motor are symmetrical about the origin of the spherical motor. Therefore, as long as the torque data in the area under the local coordinate system when a single coil of the upper or lower layer is energized can be obtained, the data of another layer can be converted, and the torque data in other coil coordinate systems on the same layer can be obtained by “analogy”, as shown in Eq. (5). In this way, the torque value in the whole spherical motor area can be obtained quickly. $\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}T_{tx}=-T_{bx}\\ T_{ty}=-T_{by}\\ T_{tz}=T_{bz},\end{array}\right. & & \displaystyle\end{eqnarray}$ (5) where T _tx, T _ty, T _tz represent the three-degrees-of-freedom torque of the upper layer; and T _bx, T _by, T _bz represent the three-degrees-of-freedom torque of the lower layer.

3.2. Torque map data acquisition

In this paper, the torque data are prestored and interpolated to calculate the torque value of the motor movement process. It is necessary to obtain a certain amount of torque data first. Through comparison, it is found that the calculation of torque data in prestored rotor coordinate system is large and time-consuming. Taking 1° as the sampling interval, the number of sampling points is 1980. The selection of torque data in prestored stator coordinate system is relatively simple, with 1485 sampling points.

In the experiment, if the current is greater than 3 A, it is easy to damage the driving circuit equipment, and the current exceeds 2 A has reached the saturation state. For the middle and upper coils, a single coil M1 and T1 are powered on, and the currents are 0.5 A, 1 A, 1.5 A, 2 A, 2.5 A and 3 A respectively. Taking 1° as the sampling density in the region, each region contains 1485 sampling points, which are simulated and calculated by ANSYS. The obtained torque data are stored in 18 tables in the upper and middle layers, and the lower torque is obtained by conversion. The final torque map can be obtained by superimposing the torque under different current conditions, as shown in Fig. 7, where lon and lat are the longitude and latitude in the stator coordinate system. Figures 7(a)–7(c) show the torque in the three coordinate axes of the middle single coil, Figs 7(d)–7(f) show the torque in the three coordinate axes of the upper single coil, and Figs 7(g)–7(i) show the torque in the three coordinate axes of the lower single coil.

Fig. 7.

Three layers of single coil torque waveform diagram.

Fig. 8.

(Continued).

3.3. Coordinate transformation and torque linear interpolation

According to the verification conclusion in Section 2., the local coordinate system and the global coordinate system can be converted to each other. The corresponding torque data are also obtained in the previous section. Therefore, when any coil is energized, the torque at each position in the process of rotor motion can be calculated by linear interpolation through the stored torque value, and finally the torque in the global coordinate system can be obtained by projection conversion.

For the rotor, the spherical motor center is the coordinate origin, the Z-axis direction is the output axis direction, the X-axis points to the center of rotor pole 1, and the Y -axis points to the midpoint of rotor pole 2 and 3, where O − dqp represents the rotor coordinate system as shown in Fig. 9.

Fig. 9.

Rotor marking diagram.

Firstly, the initial position information of six rotor poles shall be determined, and the longitude and latitude position information of the rotor pole shall be converted into three-degrees-of-freedom coordinates to facilitate the quaternion calculation of the rotor coordinates. For any point P (x, y, z) between the spherical surface and the spherical center L, as shown in Fig. 10, it can be transformed through the relationship between angles, and its orthogonal coordinate system and longitude and latitude conversion formulas are as shown in Eqs (6) and (7). $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \left\{\begin{array}{@{}l@{}}\displaystyle {\theta}_{lon}=\arctan \left(\frac{z}{\sqrt{x^{2}+y^{2}}}\right)\\[15.0pt] \displaystyle {\theta}_{lat}=\arctan \frac{y}{x}.\end{array}\right.\end{eqnarray}$ (6) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \left\{\begin{array}{@{}l@{}}x=l\times \cos {\theta}_{lat}\times \cos {\theta}_{lon}\\ y=l\times \cos {\theta}_{lat}\times \sin {\theta}_{lon}\\ z=l\times \sin {\theta}_{lat}\end{array}\right.\end{eqnarray}$ (7) where θ_lat, θ_lon represent the angle between longitude and latitude of rotor pole.

Fig. 10.

Coordinate transformation and Euler angle.

After obtaining the rotor coordinates of three-degrees-of-freedom, when the rotor rotates to a new position, the coordinates can be calculated by quaternion [19,20]; that is, the new coordinates rotated in the global coordinate system can be obtained by Eqs (8) and (9), and then converted into longitude and latitude coordinates by (6). After calculating the longitude and latitude coordinates of the rotor position during rotation, it is also necessary to judge the position. $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle Q_{\mathit{ZYX}}=[x,y,z,{\omega}]^{T}=\left[\begin{array}{@{}c@{}}\displaystyle \cos \frac{{\alpha}}{2}\cos \frac{{\beta}}{2}\sin \frac{{\gamma}}{2}+\sin \frac{{\alpha}}{2}\sin \frac{{\beta}}{2}\cos \frac{{\gamma}}{2}\\[10.0pt] \displaystyle \cos \frac{{\alpha}}{2}\sin \frac{{\beta}}{2}\cos \frac{{\gamma}}{2}-\sin \frac{{\alpha}}{2}\cos \frac{{\beta}}{2}\sin \frac{{\gamma}}{2}\\[10.0pt] \displaystyle \sin \frac{{\alpha}}{2}\cos \frac{{\beta}}{2}\cos \frac{{\gamma}}{2}+\cos \frac{{\alpha}}{2}\sin \frac{{\beta}}{2}\sin \frac{{\gamma}}{2}\\[10.0pt] \displaystyle \cos \frac{{\alpha}}{2}\cos \frac{{\beta}}{2}\cos \frac{{\gamma}}{2}-\sin \frac{{\alpha}}{2}\sin \frac{{\beta}}{2}\sin \frac{{\gamma}}{2}\end{array}\right].\end{eqnarray}$ (8) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \left[\begin{array}{@{}c@{}}x_{2}\\ y_{2}\\ z_{2}\end{array}\right]=\left[\begin{array}{@{}ccc@{}}{\omega}^{2}+x^{2}-y^{2}-z^{2} & 2(xy-{\omega}z) & 2(xz+{\omega}y)\\ 2(xy+{\omega}z) & {\omega}^{2}-x^{2}+y^{2}-z^{2} & 2(yz-{\omega}x)\\ 2(xz-{\omega}y) & 2(yz+{\omega}x) & {\omega}^{2}-x^{2}-y^{2}+z^{2}\end{array}\right]\left[\begin{array}{@{}c@{}}x_{1}\\ y_{1}\\ z_{1}\end{array}\right],\end{eqnarray}$ (9) where 𝛼, 𝛽, 𝛾 represent three Euler angles and ω is the constraint parameter.

The area of the coil and the range of motion of the rotor are different. For the middle-level coil, the longitude direction of the rotor motion area is 360°, and the dimension is 33° north and south respectively. When the upper and lower coils are powered on, the motion range is 360° longitude direction, and the dimension direction is 33° South or North relative to the center position of the upper or lower coils. The motion area is beyond the coil area. Through analysis, it is found that the torque its corresponding motion range is zero at longitude 30° and the maximum at longitude 22.5°. When the dimension direction is 16.5°, the torque is the maximum, when the dimension is 33°, the torque is the minimum, close to 0. When the dimension exceeds 33°, the torque is 0. When the rotor moves to this area, the torque data at the required position can be obtained by multiple polynomial fitting through the prestored torque data, as shown in Eq. (10). $\begin{eqnarray}\displaystyle y_{I}(x,{\omega})={\omega}_{0}+{\omega}_{1}x+{\omega}_{2}x^{2}+\cdots +{\omega}_{n}x^{n}, & & \displaystyle\end{eqnarray}$ (10) where I represents the current value, x is the angle, y is the torque, and ω is the coefficient at this current value.

When the rotor rotates to any position within the operation range of the spherical motor, the new longitude and latitude coordinates after the rotor rotation are obtained through the combination of multiple coordinate changes and quaternion calculation. First, it is necessary to determine whether the rotor teeth are located in the coil area by Eq. (11) below $\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}|{\theta}_{lon_{i}}-{\theta}_{lon_{j}}|\leq 22.5\\ |{\theta}_{lat_{i}}-{\theta}_{lat_{j}}|\leq 16.5,\end{array}\right. & & \displaystyle\end{eqnarray}$ (11) where θ_{lon
_i} and θ_{lat
_i} represent the longitude and latitude of the rotor pole, θ_{lon
_j} and θ_{lat
_j} represent the longitude and latitude of the energized coil.

Fig. 11.

Linear interpolation graph.

If the rotor is not in a single stator area, the polynomial fitting method described above needs to be carried out to calculate the data of 8 similar points Q₁₁, Q₁₂, Q₁₃, Q₁₄, Q₂₁, Q₂₂, Q₂₃, Q₂₄ around the rotor; Where Q₁₁, Q₁₂, Q₁₃, Q₁₄ are the values closest to the lower boundary between the actual current and the storage current, and Q₁₄, Q₂₁, Q₂₂, Q₂₃, Q₂₄ are the values closest to the upper boundary between the actual current and the storage current. At the same time, the longitude and latitude positions of the rotor pole under the two currents are rounded up and down to obtain the position information of these eight points. If it is located in the area, the sampling data of 8 similar points around the rotor can be obtained directly from the prestored torque data, the position of each point is shown in Fig. 11. T₁ can be obtained by interpolation of Q₁₁ and Q₁₃, T₂ can be obtained by interpolation of Q₁₂ and Q₁₄, B₁ can be obtained by interpolation of Q₂₁ and Q₂₃, and B₂ can be obtained by interpolation of Q₂₂ and Q₂₄. Then P₁ can be obtained by interpolation of T₁ and T₂, P₂ can be obtained by interpolation of B₁ and B₂, and finally O can be obtained by interpolation of P₁ and P₂. P₁ and P₂ are the torque data under two different currents. Their calculation formula is shown in Eq. (12) to obtain the torque value under the final position and current, and the schematic diagram is shown in Fig. 11. $\begin{eqnarray}\displaystyle T(O) & {\approx} & \displaystyle \left\{\frac{I-I_{1}}{I_{2}-I_{1}}\left(\begin{array}{@{}l@{}}\displaystyle \frac{({\theta}_{lon4}-{\theta}_{2lon})}{({\theta}_{lon4}-{\theta}_{lon3})}\left(\frac{({\theta}_{lat4}-{\theta}_{2lat})}{({\theta}_{lat4}-{\theta}_{lat3})}T(Q_{21})+\frac{({\theta}_{2lat}-{\theta}_{lat3})}{({\theta}_{lat4}-{\theta}_{lat3})}T(Q_{23})\right)\\[10.0pt] \displaystyle +\frac{({\theta}_{2lon}-{\theta}_{lon1})}{({\theta}_{lon2}-{\theta}_{lon1})}\left(\frac{{\theta}_{lat4}-{\theta}_{2lat}}{{\theta}_{lat4}-{\theta}_{lat3}}T(Q_{22})+\frac{({\theta}_{2lat}-{\theta}_{lat3})}{({\theta}_{lat4}-{\theta}_{lat3})}T(Q_{24})\right)\end{array}\right)\right.\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle +∼\left.\frac{I_{2}-I}{I_{2}-I_{1}}\left(\begin{array}{@{}l@{}}\displaystyle \frac{({\theta}_{lon2}-{\theta}_{1lon})}{({\theta}_{lon2}-{\theta}_{lon1})}\left(\frac{({\theta}_{lat2}-{\theta}_{1lat})}{({\theta}_{lat2}-{\theta}_{lat1})}T(Q_{11})+\frac{({\theta}_{1lat}-{\theta}_{lat1})}{({\theta}_{lat2}-{\theta}_{lat1})}T(Q_{13})\right)\\[10.0pt] \displaystyle +\frac{({\theta}_{1lon}-{\theta}_{lon1})}{({\theta}_{lon2}-{\theta}_{lon1})}\left(\frac{{\theta}_{lat2}-{\theta}_{1lat}}{{\theta}_{lat2}-{\theta}_{lat1}}T(Q_{12})+\frac{({\theta}_{1lat}-{\theta}_{lat1})}{({\theta}_{lat2}-{\theta}_{lat1})}T(Q_{14})\right)\end{array}\right)\right\},\end{eqnarray}$ (12) where θ represent longitude and latitude of different points, T is the torque under different positions and I, I ₁, I ₂ are three different current values.

The interpolated torque value also needs to be projected. Each coil has its own local coordinate system, and the interpolated data are also in the local coordinate system. Therefore, it is necessary to project all the data in the local coordinate system to the global coordinate system, so the calculated torque data are the torque data in the global coordinate system. The projection formula is shown in Eq. (13) below $\begin{eqnarray}\displaystyle \left[\begin{array}{@{}c@{}}T_{x}\\ T_{y}\\ T_{z}\end{array}\right]=\left[\begin{array}{@{}ccc@{}}\cos \angle Xx & \cos \angle Xy & \cos \angle Xz\\ \cos \angle Yx & \cos \angle Yy & \cos \angle Yz\\ \cos \angle Zx & \cos \angle Zy & \cos \angle Zz\end{array}\right]\left[\begin{array}{@{}c@{}}T_{xi}\\ T_{yi}\\ T_{zi}\end{array}\right], & & \displaystyle\end{eqnarray}$ (13) where T _x, T _y and T _z represent the torque in the global coordinate system; T _xi, T _yi and T _zi represent the torque in the local coordinate system of the i the coil.

3.4. Finite element simulation verification

In order to verify the correctness of the model under the single coil, the rotor is located at the initial position with longitude and latitude of 0°. 1.2 A current is applied to M2 coil, and the rotor will pitch from −16° to 16° around Y -axis; B5 coil is energized with a current of 2.2 A, and the rotor spins from −30° to 30° around the Z-axis.

Fig. 12.

Comparison curves of the torque motion states: (a) spinning motion; (b) pitching motion.

The torque value calculated by single coil interpolation during motion and the comparison results of ANSYS simulation torque are shown in Fig. 12, respectively. The two results are highly consistent and the torque waveform is highly consistent. Because the friction of reluctance motor itself is small, the model is effective.

4. Multi coil torque superposition

4.1. Obtain dynamic error coefficient

The operation of the motor requires the adjacent coils to be energized. It is found that for the 24 stator coils of the reluctance spherical motor, the nonlinear interference between indirect adjacent coils is small, and the farther the distance is, the smaller the error is. For a pair of coils with a difference of 180°, the minimal influence can be ignored and superimposed directly. The formula is shown in Eq. (14). For the reluctance spherical motor, the continuous commutation of the adjacent coils makes the motor run. Therefore, it is necessary to focus on the coupling relationship between adjacent coils and iron core saturation.

The motion of the reluctance spherical motor is a combination of spin and pitch motion in order to shorten the calculation time and reduce the amount of data storage. Through the analysis of the magnetic field and simulation data of the motor, the torque of the spin motion and pitch motion of the motor in different longitude and latitude has the same regular distribution, but there are different values. The torque of the same movement under different currents also has the same regular distribution and is only different in value. However, due to the saturation of the stator core, it is verified that the motor reaches the saturation state when the current reaches 2 A in order to shorten the torque calculation time under different currents and ensure the accuracy of calculation data. Considering the piecewise linearization of the current, 0 ∼ 2 A is a section and 2 ∼ 3 A is another section. In the corresponding areas of the two sections, the difference of dynamic error parameters is very small. Therefore, it is only necessary to consider the dynamic error parameters under each current in the two sections. It is found that under the influence of nonlinearity, the higher the coincidence degree between rotor teeth and stator, the greater the superposition error.

The rotor is located at the position with zero longitude and latitude. Through single coil interpolation is used to calculate that M1 and M2 coils and their diagonal coils are energized with 1 A current respectively. Select four latitude positions, calculate the torque in the global coordinate system during the rotor motion from −25° to 25° spin motion, and then calculate the torque under the same motion when M1 and M2 are energized with 1 A current at the same time. Single coil interpolation is also used to calculate the 1 A current of M1 and T1 coils and their diagonal coils respectively. Select five longitude positions, calculate the torque of the rotor in the global coordinate system during the pitch motion from −16° to 17°, and calculate the torque under the same motion when 1 A current is applied to M1 and T1 at the same time. Equation (15) is used to calculate the dynamic error coefficient, and the results are shown in Fig. 13. When the rotor is far away from the coil, the error coefficient is very large. When the rotor moves to the center of the coil, the value of the error coefficient is very small. The reason is that the torque value is very small at this time. The instability of finite element simulation leads to obvious torque error.

Fig. 13.

Comparison curves of the error coefficient motion states: (a) spinning motion; (b) pitching motion.

Fig. 14.

Several important positions: (a) just contact; (b) 5° position; (c) alignment.

For any group of adjacent coils and their diagonal coils, select 1 A and 3 A current respectively, and calculate the torque in three important cases of spin or pitch motion, as shown in Fig. 14. The single coil calculation torque method in the previous section is used to calculate the torque of two adjacent single coils and their diagonal coils in the global coordinate system, and the rotor is in the same position for spin or pitch motion. After the superposition of two single coil torque data and the comparison with the synthetic torque of adjacent coils at the corresponding position of the simulation, the value of dynamic error coefficient under different sectional currents is calculated, and the formula is expressed as Eq. (15). $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \left\{\begin{array}{@{}l@{}}T_{that}=T_{i-that}+T_{j-that}\\ T_{\mathit{fy}}=T_{i-fy}+T_{j-fy}.\end{array}\right.\end{eqnarray}$ (14) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \left\{\begin{array}{@{}l@{}}\displaystyle x_{{\theta}1}=\frac{T_{\mathit{FEM}∼1}}{T_{\sum 1}}(0<i<2A)\\[13.0pt] \displaystyle x_{{\theta}2}=\frac{T_{\mathit{FEM}∼2}}{T_{\sum 2}}(2A\leq i\leq 3A),\end{array}\right.\end{eqnarray}$ (15) where T _i−that, T _j−that, T _i−fy and T _j−fy are the torque value of spin and pitch motion under the i or j coil; x _θ represents the error coefficient; T _∑ and T _FEM represent the synthetic torque directly superimposed by multiple coils and the torque calculated by the finite element method respectively.

4.2. Multicoil torque

During the operation of the rotor of the spherical motor, it is necessary to simultaneously energize the adjacent coils for commutation, so as to make the reluctance spherical motor spin or pitch. For a pair of adjacent coils, in order to obtain an accurate model and reduce the error. The dynamic error coefficients of spin and pitch motion obtained in the previous section are fitted by bivariate polynomial respectively, and the formulas are shown in Eqs (16) and (17). The error coefficient is “symmetrical” about the center line of the stator coil, and the corresponding error coefficient can be obtained according to the corresponding position by using polynomials. $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \left\{\begin{array}{@{}l@{}}\displaystyle {\rho}_{that}(x,y)=\mathop{\sum }_{i}^{m}\mathop{\sum }_{j=0}^{i}c_{ij}x^{j}y^{i-j}(0<i<2A)\\ \displaystyle {\rho}_{that}(x,y)=\mathop{\sum }_{i}^{m}\mathop{\sum }_{j=0}^{i}c_{ij}x^{j}y^{i-j}(2A\leq i\leq 3A).\end{array}\right.\end{eqnarray}$ (16) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \left\{\begin{array}{@{}l@{}}\displaystyle {\rho}_{\mathit{fy}}(x,y)=\mathop{\sum }_{i}^{m}\mathop{\sum }_{j=0}^{i}c_{ij}x^{j}y^{i-j}(0<i<2A)\\ \displaystyle {\rho}_{\mathit{fy}}(x,y)=\mathop{\sum }_{i}^{m}\mathop{\sum }_{j=0}^{i}c_{ij}x^{j}y^{i-j}(2A\leq i\leq 3A),\end{array}\right.\end{eqnarray}$ (17) where i and j are coefficients, m is the highest order, x and y are the spin angle and pitch angle, 𝜌_that and 𝜌_fy are the values of spin and pitch error coefficients, respectively.

According to the operation requirements of motor rotation and pitch, the torque in the global coordinate system calculated by single coil is used, and a pair of coil torque is obtained by linear superposition of diagonal coils. Then the torque is multiplied by the error coefficient corresponding to the angle and current, and finally superimposed under adjacent coils to obtain the multi coil synthetic torque. The formulas are shown in Eq. (18). $\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}\displaystyle T_{\sum \mathit{that}}=\sum {\rho}_{that}\times T_{that}\\[10.0pt] \displaystyle T_{\sum \mathit{fy}}=\sum {\rho}_{\mathit{fy}}\times T_{\mathit{fy}},\end{array}\right. & & \displaystyle\end{eqnarray}$ (18) where T _that and T _fy are the synthetic torque of a diagonal coil, T _∑that and T _∑fy are multicoil synthetic torque.

4.3. Finite element simulation verification

In order to verify the correctness of the model, M2 and M3 coils and their corresponding diagonal coils are energized with 1 A and 2 A currents, the longitude and latitude of the rotor is set to zero, and the rotor spins from 0° to 30° around the Z-axis. T5 and M5 coils and their corresponding diagonal coils are energized with 1 A current at the same time, and the rotor spins 180° horizontally and pitches from 0° to 30° around the Y -axis. The torque value calculated by multiple coils during movement and the comparison results of ANSYS simulation torque are shown in Fig. 15.

Fig. 15.

Three dimensional torque comparison of multicoil: (a) spinning motion; (b) pitching motion.

Fig. 16.

The experiment platform.

The calculated torque value is basically consistent with the output torque obtained by ANSYS simulation, but there is a certain error. The main reasons for the error between them are described below. Firstly, the simulation process of ANSYS software is unstable under the influence of computer, and there will be a certain fluctuation of torque. When the torque is very small, the influence of this fluctuation is more obvious. Secondly, in the process of motor torque superposition, there is a certain difference between the fitting error coefficient and the actual error coefficient, resulting in a small amount of error in the superposition.

4.4. Experimental verification

The experimental platform of the torque model is shown in Fig. 16, including computer, DC power supply, driving circuit, motor prototype and torque detection bench. The torque detection bench is composed of two 57 stepper motors, two torque sensors, 57 motor controller, planetary reducer, rotating support and stepper motor driver. Since the motor is symmetrical about the X- and Y -axis structure, the horizontal equipment can be directly used to control the operation of the motor around the X- and Y -axis directions. The vertical direction equipment is used to control the movement of the motor around theZ-axis direction. The 24 stator coils of the motor are independently powered by 24 driving units through STM32 chip, and each current can be controlled by the computer. In order to verify the feasibility of the torque model, the experimental data are compared with the calculated values of the model. In order to reduce the measurement error and friction error of the motor in the experiment and facilitate control, the rotor offset is manually adjusted in the spin motion experiment. Two pairs of adjacent coils are selected with longitude and latitude of 0° and longitude 45° and latitude 0° and their diagonal coils are applied with 1 A and 1.5 A current respectively. The spin motion torque is measured from 0° to 20° around the Z-axis.

The comparison of the experimental results is shown in Fig. 17. Through the comparison, it is found that the torque trend obtained by the torque model and the experiment is the same, and there is a certain error in the value, of which the maximum error is 15 m⋅Nm. Through the error analysis of the data, the root mean square error is 6.33 m⋅Nm. The main causes of experimental error are as follows: 1. The torque sensor is manually installed on the rotor output shaft, which may not be firmly installed; 2. There is also a certain friction between the bench for testing torque and the spherical motor during operation; and 3. During the experiment, the initial position of the rotor is manually calibrated, so the angle will deviate during the spin movement of the rotor.

Fig. 17.

Comparison of experimental results.

5. Conclusion

The modeling of the reluctance spherical motor is still in the initial exploration stage. The difficulty of the motor modeling lies in the magnetic saturation at the stator core and stator yoke. Due to the existence of nonlinearity, the reluctance spherical motor can not be linearly superimposed with most permanent magnet spherical motors to obtain the global torque under different currents. In this paper, the reluctance spherical motor is transformed into a single coil to model the elementary torque between the corresponding rotors. By establishing the local coordinate system of the stator, the elementary torque in the model under single coil is calculated by ANSYS software, the torque value within the motion range of the whole rotor is built by polynomial fitting, and the calculated position torque is calculated by using the derived multiple linear interpolation formula, and finally projected into the global coordinate system. The dynamic error coefficients are obtained according to several important positions and fitted by bivariate polynomials. The torques of adjacent coils are corrected and superimposed based on the dynamic error parameters to obtain the multicoil synthetic torque.

This method uses the elementary model of torque to avoid the consumption of a lot of computing resources and greatly reduces the amount of computing. At the same time, the accuracy of the modeling method is also verified by simulation and experiment.

Footnotes

Acknowledgement

This work is supported in part by the Key Project of the China National Natural Science Foundation (51637001) and the Nature Science Research Project of Anhui province (1908085QE236).

References

Kim

C.K.

, Three-degrees-of-freedom passive gravity compensation mechanism applicable to robotic arm with remote center of motion for minimally invasive surgery, IEEE Robotics and Automation Letters 4(4) (2019), 3473–3480.

Chai

Gan

and Chen

, A novel tiered type permanent magnet spherical motor and its rotor orientation measurement principle, IEEE Access 8 (2020), 15303–15312.

H.-F.

, Calculation of magnetic field for cylindrical stator coils in permanent magnet spherical motor, Journal of Electrical Engineering & Technology 13(6) (2018), 2158–2167.

Kasashima

, Torque control method of an electromagnetic spherical motor using torque map, IEEE-ASME Transactions on Mechatronics 21(4) (2016), 2050–2060.

Zhao

Modeling and simulation of a permanent magnetic spherical motor based on solidworks and ADAMS, in: 2015 IEEE 10th Conference on Industrial Electronics and Applications (ICIEA), IEEE, 2015.

Research on torque optimization of the spherical motor based on SVM, in: 2016 IEEE 11th Conference on Industrial Electronics and Applications (ICIEA), IEEE, 2016.

Zhang

B.-J.

3D magnetic field finite element analysis of dual-stator PM spherical motor, in: 2019 31st Chinese Control And Decision Conference (CCDC), IEEE, 2019.

Liu

J.M.

, Magnetic field modeling and validation for a spherical actuator with cylindrical permanent magnets, Simulation Modelling Practice and Theory 98 (2020).

Xia

C.-L.

, Design and analysis of a variable arc permanent magnet array for spherical motor, IEEE Transactions on Magnetics 49(4) (2013), 1470–1478.

10.

Zhou

, Torque calculation of permanent-magnet spherical motor based on permanent-magnet surface current and Lorentz force, IEEE Transactions on Magnetics 56(5) (2020).

11.

Lee

, Newly structured double excited two-degree-of-freedom motor for security camera, IEEE Transactions on Magnetics 44(11) (2008), 4041–4044.

12.

Kim

Y.-B.

Lee

J.-S.

and Kwon

B.-I.

, Experiment and torque modeling of double-excited, two-degree-of-freedom motor based on magnetic equivalent circuit analysis, Journal of Electrical Engineering & Technology 8(1) (2013), 130–136.

13.

Zhang

Mathematical modeling and simulation of spherical switched reluctance motor with magnetic levitation, in: 2012 International Conference on Frontiers of Mechanical Engineering and Materials Engineering (MEME), 2012.

14.

Z.-J.

Description of Rotor’s position of magnetic suspension spherical reluctance motor and detection, in: 2012 International Conference on Frontiers of Mechanical Engineering and Materials Engineering (MEME), 2012.

15.

Tao

W.-Q.

Design and analysis of a novel spherical motor based on the principle of reluctance, in: 2018 IEEE International Power Electronics and Application Conference and Exposition (PEAC), IEEE, 2018.

16.

L.-F.

, Nonlinear modeling analysis of a reluctance spherical motor based on torque-angle characteristics, International Journal of Applied Electromagnetics and Mechanics 62(3) (2020), 577–593.

17.

Wang

, Structural design of salient pole type reluctance spherical motor, Electric Machines and Control 25(1) (2021), 90–102.

18.

J.-Z.

, Sensorless posture detection of reluctance spherical motor based on mutual inductance voltage, Applied Science-Basel 11(8) (2021).

19.

Chen

, Mapping relationship between unit quaternion, rodrigues transform model and Euler angle, Beijing Surveying and Mapping 34(1) (2020), 44–50.

20.

Xia

X.-W.

, Discussion on the transformation between attitude quaternion and Euler angles in large-angle maneuver scenario, Missiles and Space Vehicles 321(5) (2012), 47–53.