*2.2. Stator*

The stator of the reaction sphere is fixed with the spacecraft body, which consists of the curved stators and the electromagnets. They are responsible for rotation motion and magnetic suspension respectively. The curved cores and the electromagnet cores are used to conduct the magnetic flux generated by the stator coils. To meet the requirement of generating tri-axial torques, the stator must be designed as a spherically symmetric structure. As shown in Figure 3, the stator of the reaction is formed by three circles that are orthogonal to each other. Each circle is composed of 4 curved cores and enables the rotor to rotate about one certain axis. The output torque about any axis can be realized by combining the torques generated by the three circles. On the other hand, in order to eliminate the mechanical friction and extend its working life, six electromagnets are set in the position of the six intersection of three circles to get the spherical rotor magnetically levitated in the center of the stator.

**Figure 3.** Structural design of the stator cores: (**a**) original design and (**b**) improved design.

In the original design with round-shape electromagnet cores, as shown in Figure 3a, only five complete teeth on the curve iron-core are functional, while the two teeth near the side of the curve stator are not. The notches are formed by connecting three pieces of iron-core curve stator to form one complete stator circle. It causes non-uniformity eddy current and magnetic field induced along the stator circle. To reduce such effect and get the stator magnetic fields utilized in a more sufficient way, the structure of the electromagnet cores are revised to be the cross-shape. In this way, the curved iron-core and electromagnet core can be assembled together and so that two extra teeth on the side are functional, as shown in Figure 3b.

### *2.3. Winding*

For the convenience of speed adjustment and the feature of self-starting, the stator coils are rolled as three-phase windings. In order to generate the equivalent tri-axial torques, the stator should be symmetric in both mechanical and electrical designs. Since the uniaxial plane is divided into four equal parts, as shown in Figure 4, 12 coils are missing in each circle due to the occupation of the electromagnets. All in all, 1/3 of the driving coils are replaced by electromagnet coils.

**Figure 4.** Section diagram of the reaction sphere (one-DOF): (**a**) 36 toroidal coils. (**b**) 24 toroidal coils and 4 electromagnet coils.

As can be seen in Figure 5a,c, the three-phase windings come to two types, 6-pole and 4-pole. 6-pole windings can be regarded as removing 12 coils directly from the complete set of three-phase windings, which does not make any adjustment of the rest windings. As for 4-pole windings, the phase of next coil starts right from the end of the last one so that the windings on the four curved cores can achieve consistency.

With finite element method (FEM) modeling in COMSOL Multiphysics 5.3, the fields and currents distribution on the surface of the rotor can be visually analyzed. As shown in Figure 5b,d, the gray lines stand for the magnetic line of force, the red arrows refer to the induced currents and the surface colors refer to the radius magnetic flux density. In this simulation, copper thickness, steel thickness, ampere turns, tooth thickness and tooth width are set as 0.3 mm, 3 mm, 2.4 kA, 3 mm and 25 mm respectively. In this case, the driving torques of the 6-pole arrangement and 4-pole arrangement are 33.8 mNm and 27.0 mNm, where the former is 1.25 times larger than the latter. However, it can be seen that, in 6-pole scheme, two poles are much larger than the rest four and the field distortions in the regions that face the electromagnets are serious. Also, it induces an asymmetric distribution of eddy currents by the adjacent curved stators. These features increase the difficulty of multi-DOF motion control of this 6-pole reaction sphere. Furthermore, from the radius magnetic flux density in Figure 5b,d, it can be seen that, in 6-pole arrangement it is the electromagnet cores that mainly contribute to the driving torque, while in 4-pole arrangement it is the curved cores that mainly contribute to the driving torque. The 4-pole arrangement makes the curved cores better exploited. Additionally, the case without electromagnets cores is simulated, where the torque of 4-pole scheme declines from 27 mNm to 18.4 mNm, and the torque of 6-pole scheme declines from 33.8 mNm to 9.89 mNm. This phenomenon reflects the generated torque is stronger dependent on the electromagnet cores in 6-pole scheme. The flux generated by this 6-pole arrangement mainly passes through the

inner core of the electromagnets, which can interfere the magnetic suspension and increase the risk of magnetic saturation.

In summary, while reducing the output torque (20% reduction rate in case of the Figure 5), 4-pole arrangement brings the benefits in three aspects compared with 6-pole arrangement. Firstly, it generates a more uniform distribution of magnetic field and eddy currents on the surface of the rotor. Secondly, it is beneficial for larger angular momentum storage due to its higher synchronous speed with the same power frequency. Thirdly, it reduces the risk of magnetic saturation when the curved inductors and electromagnets are working at the same time. Thus, the 4-pole scheme is adopted for the reaction sphere windings.

**Figure 5.** Winding mode and simulation of magnetic field and induced currents on the surface of the rotor: (**a**) 6-pole scheme. (**b**) Simulated distribution of 6-pole scheme. (**c**) 4-pole scheme. (**d**) Simulated distribution of 4-pole scheme.

#### **3. Analytical Model**

This section analyzes the magnetic fields on the surface of the rotor, which leads to the solution to electromagnetic driving torque. This also allows us to analyze the generated torque with respect to some of the key design parameters, such as number of poles. Several simplifications are made here. (1) The magnetic permeability of the rotor core is infinite. Thus, no magnetic saturation is encountered. (2) The length of air-gap is evenly distributed. (3) Displacement currents are neglected. Thus, the current law are simplified as Ampere's law.

Generally speaking, the magnetic flux density can be calculated by computing the curl of magnetic vector potential. In spherical coordinates (*r*, *θ*, *ϕ*), the formula can be expressed as

$$
\begin{array}{c|cccc}
\begin{array}{c|cccc}
& & & & \overline{\mathcal{e}\_r} & r\overline{\mathcal{e}\_\theta} & r\sin\theta\,\overline{\mathcal{e}\_\phi} & & \\
\hline
r^2\sin\theta & \frac{\partial}{\partial r} & \frac{\partial}{\partial \theta} & \frac{\partial}{\partial \phi} & & \\
& A\_r & rA\_\theta & r\sin\theta A\_\phi & & \\
\end{array}
\end{array}
$$

where −→*er* , −→*e<sup>θ</sup>* and −→*e<sup>ϕ</sup>* refer to the unit vectors, *Ar*, *A<sup>θ</sup>* and *A<sup>ϕ</sup>* refer to components of the vector potential in spherical coordinate system.

The general torque −→*τ* of the reaction sphere is the vector summation of 3 torques −→*τ<sup>x</sup>* , −→*τ<sup>y</sup>* , and −→*τ<sup>z</sup>* generated by 3 identical circles of the stator.

$$
\overrightarrow{\tau} = \overrightarrow{\tau\_x} + \overrightarrow{\tau\_y} + \overrightarrow{\tau\_z} \dots
$$

The following analysis will take torque about z-axis −→*τ<sup>z</sup>* for example. It is the longitudinal magnetomotive force (mmf) that induces the z-axis rotating magnetic field, in which case *A<sup>ϕ</sup>* = *Ar* = 0. Therefore, the magnetic flux density can be simplified as

$$
\overrightarrow{B}' = -\frac{\overrightarrow{\varepsilon\_r}^{\flat}}{r\sin\theta} \frac{\partial A\_\theta}{\partial \rho} + \frac{\overrightarrow{\varepsilon\_\theta}^{\flat}}{r} \frac{\partial r A\_\theta}{\partial r} \equiv B\_r \overrightarrow{\varepsilon\_r}^{\flat} + B\_\phi \overrightarrow{\varepsilon\_\phi}.\tag{2}
$$

In spherical coordinate system, by Ampere's law −→*<sup>J</sup>* <sup>=</sup> <sup>Δ</sup> <sup>×</sup> −→*<sup>H</sup>* ,

$$
\begin{array}{c|ccc}
\hline
\overline{J} & \frac{1}{r^2 \sin \theta} & \frac{\overline{\partial}}{\partial r} & r \sin \theta \, \overline{\mathcal{E}}\_{\overline{\theta}}^{\flat} \\
\hline
\overline{J} & \frac{\partial}{\partial r} & \frac{\partial}{\partial \theta} & \frac{\partial}{\partial \overline{\theta}} \\
H\_r & rH\_\theta & r \sin \theta H\_\varphi & \\
\hline
\end{array}
$$

In this case, only those currents flowing in *θ*-direction induce torque about z-axis,

$$J\_{\theta} = -\frac{1}{r\sin\theta} \left( \frac{\partial r\sin\theta H\_{\varphi}}{\partial r} - \frac{\partial H\_{r}}{\partial \varphi} \right) = -\frac{1}{\mu\_{0}r\sin\theta} \left( \frac{\partial r\sin\theta B\_{\varphi}}{\partial r} - \frac{\partial B\_{r}}{\partial \varphi} \right). \tag{3}$$

Combine (3) and (2),

$$J\_{\theta} = -\frac{1}{\mu\_0 r} \frac{\partial^2 r A\_{\theta}}{\partial r^2} - \frac{1}{\mu\_0 r^2 \sin^2 \theta} \frac{\partial^2 A\_{\theta}}{\partial \phi^2}. \tag{4}$$

Express the magnetic potential *A<sup>θ</sup>* written in complex form *ejω<sup>t</sup>* [27], where *ω* is the angular frequency of input currents, *j* is the imaginary unit. Then,

$$J\_{\theta} = \sigma E = -\sigma \frac{\partial A}{\partial t} = -j\sigma \omega A\_{\theta}. \tag{5}$$

Combine (4) and (5), the differential equation for magnetic field distribution on the surface of the rotor is obtained as

$$
\sigma j \omega A\_{\theta} = \frac{1}{\mu\_0 r} \frac{\partial^2 r A\_{\theta}}{\partial r^2} + \frac{1}{\mu\_0 r^2 \sin^2 \theta} \frac{\partial^2 A\_{\theta}}{\partial \phi^2}. \tag{6}
$$

To make (6) solvable, apply the the method of separation of variables as

$$A\_{\theta}(r,\theta,\varphi) = X(r)\,\chi(\theta,\varphi). \tag{7}$$

Substitute (7) into (6),

$$
\sigma j \omega = \frac{1}{\mu\_0 r X} \frac{\partial^2 r X}{\partial r^2} + \frac{1}{\mu\_0 r^2 \sin^2 \theta Y} \frac{\partial^2 Y}{\partial \varphi^2}. \tag{8}
$$

With high-order time and space harmonics neglected, it is assumed that *ejω<sup>t</sup>* only has the first mmf space harmonic. In this way,

1 *Y* sin2 *θ ∂*2*θY ∂ϕ*<sup>2</sup> <sup>=</sup> <sup>−</sup>*p*2, (9)

where *p* is the number of pair pole. The analytical solution to (9) is given as

$$Y = \mathbb{C}\_1 e^{j p \cdot \rho \sin \theta} + \mathbb{C}\_2 e^{-j p \cdot \rho \sin \theta}. \tag{10}$$

From (8) and (9), it yields

$$r^2 \frac{\partial^2 X}{\partial r^2} + 2r \frac{\partial X}{\partial r} - X(r^2 \mu \sigma j \omega + p^2) = 0. \tag{11}$$

Set X(*r*,*θ*) = <sup>√</sup><sup>1</sup> *r* X1(r,*θ*). Then, (11) are transformed to modified Bessel's Equation (12) [28],

$$r\_1 \frac{\partial^2 X\_1}{\partial r\_1^2} + r\_1 \frac{\partial X\_1}{\partial r\_1} - X\_1 (r\_1^2 + p^2 + \frac{1}{4}) = 0,\tag{12}$$

where *r*<sup>1</sup> = *r* !*μσjω*. The analytical solution to (12) is as follows:

$$X(r) = \frac{1}{\sqrt{r}} X\_1(r, \theta) = \frac{1}{\sqrt{r}} \mathcal{C}\_3 I \frac{1}{\sqrt{r^2 + \frac{1}{4}}} (r \sqrt{\mu r j \omega}) + \frac{1}{\sqrt{r}} \mathcal{C}\_4 K \frac{1}{\sqrt{r^2 + \frac{1}{4}}} (r \sqrt{\mu r j \omega}), \tag{13}$$

where *Ip*(*r*1) and *Kp*(*r*1) are the two linearly independent solutions to the modified Bessel's Equation [28]:

$$\begin{split} I\_a(\mathbf{x}) &= \sum\_{m=0}^{\infty} \frac{1}{m!\Gamma(m+\mathfrak{a}+1)} (\frac{\chi}{2})^{a+2m} , \\ K\_a(\mathbf{x}) &= \frac{\pi}{2} \frac{I\_{-a}(\mathbf{x}) - I\_a(\mathbf{x})}{\sin a\pi} . \end{split} \tag{14}$$

Subsequently, the magnetic flux density can be calculated through (15), with boundary conditions of *Hϕ*|*r*=*Rs* = 0 and *Hϕ*copper|*r*=*Rc* = *Hϕ*air|*r*=*Rc* , *Br*copper|*r*=*Rc* = *Br*air|*r*=*Rc* for continuity, where *Rc* is the radius of rotor outer surface and *Rs* is the radius of rotor inner surface. The nickel layer is neglected because it is much thinner than copper layer and steel layer.

$$\begin{split} B\_{\mathsf{T}} &= -\frac{\mathsf{X}(r)}{r \sin \theta} \frac{\partial \mathsf{Y}(\theta, \varphi)}{\partial \varphi}, \\ B\_{\mathsf{F}} &= \frac{\mathsf{X}(r) \mathsf{Y}(\theta, \varphi)}{r} + \mathsf{Y}(\theta, \varphi) \frac{\partial \mathsf{X}(r)}{\partial r}. \end{split} \tag{15}$$

Based on the flux density distribution, the driving torque can be evaluated by Lorentz forces induced by eddy currents on the copper layer of the rotor. The electromagnetic torque arising from the rotor steel is negligible [27].

$$T\_{eL} = -\int\_0^{2\pi} \int\_{R\_s}^{R\_c} \int\_0^{\pi} f\_\theta B\_r r^3 \sin^2\theta \,d\theta dr d\varphi. \tag{16}$$

In case of proposed design, 16 of the 36 teeth in circle are occupied by 4 electromagnets (see Figure 4a). After structural optimization (see Figure 3b), each side tooth of these curved cores can be approximated as a half tooth concerning its combination with the electromagnet core. This helps form 4 equivalent teeth. So the output torque generated by the four curved inductors is 2/3 of the complete electromagnetic torque. Then, substitute (5), (7) and (10) into (16),

$$T\_{\varepsilon LRS} = \frac{2}{3} T\_{\varepsilon L} \tag{17}$$

$$I\_{\theta} = \frac{2}{3} \int\_{0}^{2\pi} \int\_{R\_{s}}^{R\_{c}} \int\_{0}^{\pi} J\_{\theta} \frac{\partial A\_{\theta}}{\partial \varphi} r^{2} \sin \theta \, d\theta dr d\varphi \tag{18}$$

$$\mathcal{L} = \frac{2}{3} \sigma p \omega \int\_{R\_4}^{R\_c} X^2(r) r^2 dr \int\_0^{2\pi} \int\_0^{\pi} \chi^2(\theta, \varphi) \sin^2 \theta d\theta d\varphi \tag{19}$$

$$=\frac{\pi^2}{3}\sigma p\omega \int\_{R\_\delta}^{R\_\epsilon} X^2(r)r^2 dr. \tag{20}$$

From (20), it can be seen that higher electrical conductivity and pole numbers can enhance the output torque. So it is reasonable to adopt high conductivity materials like copper as the outer layer of the rotor. Also, the torque is expected to be reduced when choosing the 4-pole winding instead of 6-pole winding. To raise the output torque, the design parameters like tooth thickness and tooth width are required to be optimized in developing the prototypes. However, the analytical torque model (20) does not include these parameters. To authors' best knowledge, all the proposed analytical fields and torque models of inductive reaction spheres are developed based on the assumptions that the stator is slotless [27,29–31]. Otherwise, the differential equation will be too complex to find an analytical solution. Due to the nonlinear dependency of design parameters and the 3-D fields distribution, it is too complicated to derive a purely analytical torque model if the stator contains slots. Hence, data-based modeling technique is used to predict the output torque with design parameters considered to be optimized in the following section.

#### **4. Torque Density Optimization**

#### *4.1. Data-Based Torque Model*

In this optimization, by dimension constraints, the air gap length is fixed as 1 mm, the outside diameter of the rotor is fixed as 99 mm and the outside diameter of the stator is fixed as 123 mm. The goal is to raise the output torque under the constraints of mass and MOI by searching the optimal tooth thickness (*S*), tooth width (*B*), ampere turns (*N I*), copper thickness (*Tco*) and steel thickness (*Tst*) as shown in Figure 6. Due to the time consuming of FEM calculation, the comprehensive design is not practical. Thus, Taguchi design of five factors five levels is adopted to select the representative samples data necessary for building learning models as shown in Table 1 [32]. Consider the overall size and the process feature, the ranges of *Tco*, *Tst*, *N I*, *S* and *B* are set as 0.1–0.5 mm, 1–5 mm, 1.5–2.7 kA, 2–4 mm and 20–30 mm respectively.

Eight groups of the Taguchi design form the 200 trials with no repeated experiments. Each group consists of 25 samples, which selects the representative samples from the whole parameter space by following the two rules. (1) The value in column has the same number of occurrences. (2) Any pair of values between two arbitrary columns have the same number of occurrences. Subsequently, of all the eight groups, seven groups are used to train the learning model and the rest one is used to test the prediction performance of the model. The output torque is obtained by FEM calculation. To reveal the uneven distribution of the air gap length and raise the accuracy of the calculation, 3D FEM is adopted instead of 2D one. The whole mass is estimated by the CAD model volume and its density. The materials used in FEM simulation are listed in Table 2. The rated current is set to 1.4 A. The power source is 50 Hz.

**Figure 6.** Design parameters to be optimized.

**Table 1.** The parameters of sample space.


**Table 2.** Material components of FEM and lab prototype.


The various parameters of the inductive RS have high nonlinear relationships due to the 3-D field distribution and electromagnetic conversion. Fitting neural network (FNN) has the strength of realizing the complex non-linear mapping between the input and the output. So it is suitable for solving the problem of inductive RS torque prediction with complicated internal mechanism. Due to the large calculation and time consumption of 3D finite element simulation, obtaining a large number of training samples are not practical. Support vector machines (SVM) fit well for the small or medium size of training samples. It determines the final results with only a few support vectors. The major redundant samples are eliminated after training. Random forest (RF) can handle data with many features and have high training speed with independent decision trees. This fits the efficient multi-variable optimization of the inductive RS. Also, it can reduce the risk of over-fitting by increasing tree numbers. It is hard to give strict conclusion of application domains for the three methods. In practice, we choose the appropriate method according to the specific situation.

To enhance the prediction performance, the tree numbers of RF, the hidden layer size of the FNN and parameters of the SVM are optimized. As shown in Figure 7, increasing the number of trees generally improves the accuracy of RF and prediction tends to be stable. It reaches the highest average accuracy of 0.0022 Nm with 44 trees. As for fitting neural networks, it takes a row vector of N hidden layer sizes, and a BP training function, and returns a feed-forward neural network with N + 1 layers. Theoretically, multi-layer fitting neural networks can fit any finite input-output relationship well given enough hidden neurons. However, the risk of over fitting increases with the increase of

hidden neurons, where it may fit well for the training data but badly for the test data. Thus, the hidden layers should be chosen properly. As shown in the second figure of Figure 7, it reaches the highest average accuracy of 0.0011 Nm with 10 hidden layers. As for SVM, the radial basis function is set as the kernel function considering the nonlinearity of the torque model. The SVM parameter cost (*c*), radial basis function parameter (*g*) and epsilon in loss function (*p*) are chosen to be optimized. This optimization is carried out by genetic algorithms. To avoid over fitting, *c*, *g* and *p* are constrained within 100, 100 and 1 respectively. It shows that the SVM prediction model reaches the highest average accuracy of 0.00093 Nm when *c p* and *g* are 70.716, 0.019 and 0.066 respectively.

**Figure 7.** Data-based models parameters optimization: (**a**) RF regression with increasing trees, best root-mean-square error (RMSE) of 0.0022 Nm with 44 trees. (**b**) NN regression with increasing hidden layers, best RMSE of 0.0011 Nm with 10 hidden layers. (**c**) SVM regression with genetic algorithm optimizing parameters, best RMSE of 0.0009 Nm with *c* (70.716), *p* (0.019) and *g* (0.066).

All in all, the regression results of the test set show that the SVM model achieves the least RMSE over the other two in the case of small-sample database. The optimal predictions by RF, FNN and SVM are shown in Figure 8a, with these test samples arranged by Taguchi design of five factors five levels. Figure 8b shows the prediction errors of RF, FNN and SVM for each sample. It can be seen that the

SVM torque model has a stable prediction, with all errors restricted in 2 mNm. Thus, the SVM model is chosen to be the objective function of the optimization in next subsection.

**Figure 8.** RF, FNN and SVM models prediction performance with test data: (**a**) torque predictions and (**b**) absolute errors.

#### *4.2. Design Parameters Optimization*

The reaction sphere design parameters optimization is a multi-parameter optimization problem with constraints. The objective function is the SVM torque prediction model established above. The constraints include the box constraints, linear inequality constraints and nonlinear inequality constraints. The box constraints (22)–(26) are set to be consistent with the range of the sample spaces shown in Table 1. The linear inequality constraints are considered between coil turns, tooth thickness and tooth width: The coil turns and tooth thickness are limited by the tooth pitch (27); The coil turns and tooth width are limited by the fixed stator arc length (28). The parameters of the linear inequality constraints are validated by the lab possessed parts shown in Figure 9. The nonlinear inequality constraints consider the application in micro spacecrafts: the MOI of the rotor is larger than 0.001 kg·m2 (29) and the whole mass of the actuator is no more than 3 kg (30).

**Figure 9.** Parts of the reaction sphere prototype for constraints and masses validation.

The mass function M(*Tco*, *Tst*, *N I*, *S*, *B*) evaluates mass of the whole actuator, consisting of rotor, stator core, toroidal winding and electromagnets. The mass of the rotor and stator core is estimated by their measured volume of the digital model and their corresponding densities. The mass of the toroidal winding is calculated by the product of turn number and mass of each turn. The mass of the electromagnets is a constant, which is weighed with prototypes. Several sample components are selected to validate its accuracy as shown in Table 3. Due to the insulating coating, the mass test results of the curved cores are slightly larger than the estimated ones.

**Table 3.** Mass estimation and experimental validation.


Based on the analysis above, the mathematical description of the optimization problem can be expressed as:

	- *s*.*t*. 0.1 mm ≤ *Tco* ≤ 0.5 mm, (22)
		- 1 mm ≤ *Tst* ≤ 5 mm, (23)

$$1.5\,\text{kA} \le N I \le 2.7\,\text{kA},\tag{24}$$

$$2\text{ mm} \le S \le 4\text{ mm},\tag{25}$$

$$\text{20 mm} \le B \le \text{30 mm},\tag{26}$$

$$1\,\mathrm{NI}/33.6 \le 140-20\,\mathrm{S},\tag{27}$$

$$1\text{ NI}/33.6 \le 340 - 10S - 8B,\tag{28}$$

*Ico* <sup>+</sup> *Ist* <sup>≥</sup> 0.001 kg · <sup>m</sup>2, (29)

$$\mathbf{M}(T\_{\varepsilon\alpha\prime}T\_{st\prime}\,\mathrm{NI},\mathrm{S},\mathrm{B}) \le \mathfrak{Z}\mathrm{kg}.\tag{30}$$

This constrained optimization is carried out by interior point methods on MATLAB2016a optimization toolbox, with parameter setting as follows: function tolerance and constraints tolerance: 10−6, max iterations: 1000. Inner max iterations: 100, relative tolerance: 10−2. The optimization results are listed in Table 4. The iteration process is shown in Figure 10. The final *Tco*, *Tst*, *N I*, *S*, *B* for developed 4-pole prototype are rounded as 0.5 mm, 2.4 mm, 2.688 kA, 3 mm and 27.5 mm respectively. The driving torque is 0.0371 Nm (FEM) with the whole mass of 3 kg and MOI of 0.001 kg·m2.

**Figure 10.** Design parameters optimization by interior point method for better output torque.



#### **5. Feasibility Test with Lab Prototypes**

To validate the feasibility of the proposed design, the experimental setup is developed to test the speed performance of the lab prototypes, which is composed of converters, DC power supply, reaction sphere actuators, controllers and amplifiers as shown in Figure 11. The three converters are responsible for driving the reaction sphere actuators in the three orthogonal axis respectively. The DC power supply provides direct currents for amplifiers and controllers for magnetic levitation control.

In this magnetic levitation system, to get the feedback position signal in the closed-loop control, an eddy current displacement sensor made of tiny coils is placed at the center of each electromagnet. These sensors provide the signal of air gap distance to the controller. From here, a negative feedback control loop is set up by two serial-link PD controller to regulate the levitation distance. Specifically, the controller uses STM32F051K6T6 (Suzhou Joint-Tech Mechanical and Electrical Technology Co., LTD, Suzhou, China), which collects the signal of eddy current displacement sensors and generates PWM signal to drive the power amplifier. Each driver amplifies the control signal and drives the two electromagnetic coils on the same axis. What should be mentioned is that the electromagnets about Z-axis need to overcome the gravity of the rotor. They require higher currents compared with the counterparts about the other two axes. A 4-pole reaction sphere prototype and a 6-pole reaction sphere prototype are developed to test the speed performance with the same controllers and drivers. A binary painting strategy for speed measurement is presented as shown in Figure 12. The four silvery white circles on the surface of the rotor are located on the four vertexes of a tetrahedron. In case of measuring the X, Y and Z rotational speed, the laser focus will come across the high reflective areas, which is the nickel protection layer shown in Figure 2, 3 times when the rotor turn around once. Therefore, the output speed can be estimated as the one-third of counts per minute. The frequency of AC is set to 60 Hz in X, Y and Z axis respectively of the 4-pole reaction sphere prototype and in X and Y axis respectively of the 6-pole reaction sphere prototype. The Z axis of the 6-pole reaction sphere is set to 54 Hz, for higher frequency cannot generate a stable rotation with its significantly increasing vibration. This results from the parts manufacturing and assembling error, as well as the dynamical stability of the sphere about Z-axis rotation under gravity force for 6-pole RS, probably due to its unequal pole. The speed test results are listed in Table 5. It is shown that the 4-pole reaction sphere prototype has the maximum of 775 rpm, 791 rpm and 795 rpm in X, Y and Z axis respectively, which obtains a better speed performance over the 6-pole one as predicted. As for speed uncertainties, the 4-pole prototype shows a relatively smaller fluctuations, which may benefit from the uniform fields distributions.

**Figure 11.** Experimental setup for driving the 4-pole reaction sphere and 6-pole reaction sphere under magnetic levitation.

**Figure 12.** Spherical rotor with black painting for speed counts: (**a**) diagram and (**b**) prototype.

**Table 5.** Maximum speed of the four-pole reaction sphere prototype and the six-pole reaction sphere prototype.


The speed control experiment is not so successful. The tested output speed is much less than the synchronization speed (1800 rpm @ 60 Hz). When the frequency is over 60 Hz, the increasing vibration leads to the mechanical collision between the rotor and the stator. One reason lies in the quality of the prototypes. The magnetic levitation demands for the high precision of parts processing and

assembling, which is not well guaranteed in the developed prototypes. This results in the non-unique performance in different axes, which is particularly evident between Z axis and the other two axes in the 6-pole prototype.

A more essential reason lies in a drawback of the proposed design: When producing the driving torque, the rotation windings on the curved inductor will generate an extra thrust in the meantime. This extra thrust will push the electromagnets to generate the flux and attract the rotor to be back in the center of the stator. Then, additional damping torques will be induced by the interaction between the rotating rotor and the electromagnets flux, which leads to speed reduction.

#### **6. Conclusions**

This paper presents the new design, modeling and parameter optimization methods of a magnetically levitated inductive reaction sphere. It provides a simple, low-cost and miniaturized solution to spacecraft attitude control via momentum exchange. The compact design of electromagnets and curved inductors has the dimensional ratio of 1:1.4 between the rotor and the stator assembly. The simplified, analytical driving torque model is derived based on Ampere's law and separation of variables, showing that the conductivity of the rotor outer layer and pole number of the stator fields contribute to the driving torque. To raise the agility of attitude adjustment of satellites, it is crucial to improve the torque density of this 4-pole inductive RS besides introducing slots in the iron-core curve stator. Here, its curved inductors and electromagnets are redesigned as the cross-shape, allowing two more slots near the notches of the curve stator to be functional. This also improves the uniformity of the induced eddy current and magnetic field along the stator circle. As the full analytical torque model is not available for inductive RS with slotted curve stators, a constrained optimization problem is formulated to maximize the output torque. In addition, the design parameters are efficiently solved by SVM with only small volume of training data. The lab prototype is developed for feasibility test. It shows that the optimized RS can rotate over 700 rpm in X, Y and Z axis with the angular momentum of 0.08 kg·m2/s.

Small volume and small mass are two essential issues in reaction spheres design. The proposed compact design enables the developed prototype to have the rotor mass of 0.66 kg and the whole assembly of 3 kg under the dimension constraints of 136 mm × 136 mm × 136 mm. The necessary angular momentum of the actuator is expected to be within 0.1–1 kg·m2/s for actual micro-spacecraft (10–100 kg) attitude control [33,34], and the practical design target is set as 0.4 kg·m2/s. If we maintain the design parameters in this paper, the rotor has to spin about 3800 rpm to achieve this, demanding a high level of dynamic balance. Due to the current restraint of motor drivers and the error during part fabrication and assembly, it is not practical to realize such speed with our prototypes. A more practical way is to increase both the output speed and the mass of the steel layer of the rotor appropriately. For example, we can achieve the angular momentum of 0.4 kg·m2/s by increasing the speed to 1500 rpm and the steel thickness of the hollow sphere to 9.5 mm respectively. Remarkably, the design and optimization methods proposed in this paper are still applicable in the practical reaction sphere development.

**Author Contributions:** L.Y. derived the electromagnetic modeling, optimization method and prepared the manuscript; J.Z. and J.C. proposed the original inductive RS design; S.-L.C. advised the theoretical study, designed the experiments and refined the manuscript; Y.L. and L.Y. assembled the prototypes and carried out experimental test; C.Z. and G.Y. supervised the research throughout.

**Funding:** This research was funded by Zhejiang provincial public welfare research program(Y19E070003, LGG18E070007) and the innovation team of key components and key technology for the new generation robot(2016B10016).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Proof: (A1) is a Monotonically Increasing Function**

The MOI of a hollow and uniform sphere is

$$J = J(R\_0) = \frac{2}{5} M \frac{R^5 - R\_0^5}{R^3 - R\_0^3}, \quad 0 \le R\_0 \le R,\tag{A1}$$

where *<sup>R</sup>*<sup>0</sup> is a variable, *<sup>R</sup>* and *<sup>M</sup>* are constants, *<sup>R</sup>*0, *<sup>R</sup>*, *<sup>M</sup>* <sup>∈</sup> <sup>R</sup>+. Take logarithms on both sides of (A1),

$$
\ln(f) = \ln\left(\frac{2}{5}M\right) + \ln(R^5 - R\_0^5) + \ln(R^3 - R\_0^3). \tag{A2}
$$

Take the derivative of (A2) with respect to *R*0, it yields

$$\frac{J'}{J} = \frac{R\_0^2 \left(-5R\_0^2 R^3 + 2R\_0^5 + 3R^5\right)}{(R^5 - R\_0^5)(R^3 - R\_0^3)}.\tag{A3}$$

Set *<sup>H</sup>*(*R*0) = −5*R*<sup>2</sup> <sup>0</sup>*R*<sup>3</sup> + 2*R*<sup>5</sup> <sup>0</sup> <sup>+</sup> <sup>3</sup>*R*5, so *<sup>H</sup>*(*R*0)|*R*0=*<sup>R</sup>* <sup>=</sup> 0. In addition, *dH dR*<sup>0</sup> <sup>=</sup> <sup>10</sup>*R*0(*R*<sup>3</sup> <sup>0</sup> − *<sup>R</sup>*3) < 0. In this way, for *<sup>R</sup>*<sup>0</sup> < *<sup>R</sup>*, *<sup>H</sup>*(*R*0) > 0. This yields *<sup>J</sup>* = *JH*(*R*0)/[(*R*<sup>5</sup> − *<sup>R</sup>*<sup>5</sup> <sup>0</sup>)(*R*<sup>3</sup> − *<sup>R</sup>*<sup>3</sup> <sup>0</sup>)] > 0.

This concludes *J*(*R*0) is monotonically increasing when *R*<sup>0</sup> ∈ [0, *R*). Additionally, *J*min = *J*(0) = 2 <sup>5</sup> *MR*<sup>2</sup> <sup>0</sup>, *J* < *J*(*R*) = lim *R*0→*R* 2 <sup>5</sup> *<sup>M</sup> <sup>R</sup>*5−*R*<sup>5</sup> 0 *<sup>R</sup>*3−*R*<sup>3</sup> 0 = <sup>2</sup> <sup>3</sup> *MR*<sup>2</sup> <sup>0</sup>, *J*(0) = 60%*J*(*R*).

#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
