Loss Estimation and Thermal Analysis of a Magnetic Levitation Reaction Flywheel with PMB and AMB for Satellite Application

Abstract

The magnetic levitation reaction flywheel (MLRW) is a novel actuator of spacecraft attitude control because of its significant advantages, including lack of friction and active suppression of vibration. However, in a vacuum environment, the poor heat dissipation conditions make it more sensitive to various losses and rises in temperature. Therefore, increasing temperature is the key issue for components used in space. In this study, the losses of the three kinds of heat-generating areas in the MLRW, namely, the passive magnetic bearing (PMB), the active magnetic bearing (AMB) and brushless DC motor (BLDCM), were analyzed and calculated. Based on the electromagnetic field theory, the loss model of PMB was proposed. Based on the finite element method (FEM) and Bertotti model, the loss power of the AMB and the BLDCM was obtained. The calculated loss values were brought into the FEM to calculate the temperature field distribution of the MLRW system. Then, the key factors affecting the heat dissipation of the flywheel were obtained by combining thermal network analysis with the temperature field distribution. Finally, a prototype was fabricated. The maximum estimated and experimental temperatures were 34.8 °C and 36.8 °C, respectively, both at the BLDCM stator. The maximum error was 5.4%, which validates the calculated model.

1. Introduction

Reaction flywheels supported by AMBs and PMBs are very important apparatus in spacecraft attitude adjustment [1,2]. Compared with the traditional reaction flywheel supported by mechanical bearings, the MLRW has the advantages of not making contact, not requiring lubrication, and active vibration suppression [3,4,5]. In essence, it solves the problems of contact friction, lubrication, oil pollution, and vibration of the mechanical bearings, which have negative effects on the precision and life of the satellite platform. The magnetically suspended reaction flywheel is driven by a BLDCM and supported by an AMB in two radial-direction degrees-of-freedom (DOFs) and a PMB in the other three DOFs. Due to working in a high vacuum environment, there is no heat conduction by cross-ventilation. Simultaneously, there is no contact between the rotor and stator in the MLRW. Therefore, the heat is transferred between the stator and rotor of the MLRW only by the thermal radiation. Poor heat transfer efficiency causes the temperature to increase, and potentially overheating, in the MLRW system, which is the key issue for the flywheel application [6,7].

In an MLRW, the main losses are the copper and iron losses in the AMB, PMB, and BLDCM. Therefore, in order to reduce the copper losses and iron losses of the AMB, a novel structure of a radial hybrid magnetic bearing with permanent magnet (PM) bias was proposed [8]. A function of the MB power consumption was established, and was simulated and verified using the FEM [9]. Regarding the PMB, a novel radial PMB [10] and axial-radial PMB [11] were presented to increase the stability and accuracy, and decrease the losses. The loss of the flywheel system was calculated to determine the distribution of the main heat source of the system, and an equivalent thermal network model was established based on the whole mechanical topology structure [12]. Regarding the losses of the motor, a zero cogging brushless DC motor [13] and coreless stator brushless DC motor, which have high torque density, high positional stability, and low iron loss, were designed for spacecraft applications. In ref. [14], an analytical model for predicting the iron losses in high-speed shotless PM machines is presented, and was verified by 2-D FEM. Ref. [15] presents the analytical method of calculating losses of AMBs based on the reluctance network method. The proposed method of thermal analysis of the system was compared with other popular loss estimation methods. A revised method based on the Epstein Frame Method with rationale steps for deduction of the core loss in high-speed electric machines was proposed [16]. An extended survey on the evolution and the modern approaches in thermal analysis of electrical machines is presented in [17]. The 3-D FEM is also used to analyze the temperature distribution in high-speed PM electrical machines [18].

The purpose of this study was to obtain the thermal distribution of the high-speed MLRW to avoid deterioration of the electromagnetic performance of the PM, the soft magnetic alloy of the BLDCM, and the AMB under high temperature conditions. Therefore, it is necessary to study the temperature rise of MLRW systematically. The main contributions of this study are the systematic derivation of the PMB loss model, the loss calculation of the components of the MLRW, thermal analysis, and structural optimization. The copper losses of AMB and BLDCM were calculated by combining system structure analysis with Ohm’s law. The iron losses of AMB and BLDCM are composed of eddy current loss, hysteresis loss, and excess loss. The calculation method combines electromagnetic field simulation and Bertotti’s model [19]. The temperature field distribution of the MLRW system was predicted by 3-D FEM simulation. A thermal network model was set up. The structural optimization based on 3D finite element simulation and a thermal network was predicted. Finally, a prototype was fabricated to verify the calculated model.

2. Structure Scheme and Force Analysis of the MLRW

2.1. The MLRW Structure Scheme

A schematic section view of the active-passive MLRW is shown in Figure 1. The rotor driven by the BLDCM is supported by a radial AMB and an axial PMB in five DOFs. The BLDCM is proposed because of its high efficiency, long lifetime, and low power consumption. The PM biased hybrid radial magnetic bearing is used to reduce power losses while providing high precision control, achieving a greater bearing capacity and stiffness. The PMB is proposed to reduce the axial size of the shaft, reduce the power losses, and reduce the complexity of the control system in the MLRW. The reason for choosing this structure design is to reduce the coupling degree of the magnetic circuit of the PMB, AMB and BLDCM, improve the system control accuracy, and reduce the system size and power loss.

2.2. Force Analysis of the MLRW

As shown in Figure 2, the centroid coordinate system (oxyz) of the rotor is established. Since the radial AMB controls the two DOFs of the translational motions of the flywheel’s rotor, and the other three DOFs are controlled by the axial PMB, the kinetic equations of the rotor are shown below:

(1)$\{\begin{array}{l}m\ddot{x}=f_x=(k_{ix}{i}_x+k_{x}x)\\ m\ddot{y}=f_y=(k_{iy}{i}_y+k_{y}y)\\ J_z\dot{\mathsf{\Omega }}=T_e-T_d=\left(e_a{i}_a+e_b{i}_b+e_c{i}_c\right)/\mathsf{\Omega }-T_d\end{array}$

where m is the mass of the rotor; f_x and f_y are the bearing force of the AMB in the x-axis and y-axis, respectively; k_ix, k_iy, k_x and k_y are the current stiffness and displacement stiffness of the AMB, respectively; $\dot{\mathsf{\Omega }}$ is the angular velocity of the MLRW; J_z is the moment of inertia of the rotor around the z-axis; T_e is the torque produced by BLDCM; T_d is the disturbance torque; e_a, e_b and e_c are the induced back-EMF in the stator windings of the phases of the BLDCM; i_a, i_b and i_c are the phase currents of the BLDCM.

The losses in the MLRW, including the stator winding loss and iron core losses of the AMB and BLDCM, and the iron core losses of the PMB, should be calculated for analysis of the thermal field and the overheating issue of the system.

3. Loss Estimation of the PMB

As there is no hysteresis loss in the PM, the excess current loss is ignored. In addition, because the other structures in the PMB are non-magnetic materials, the PMB loss is mainly due to the eddy current loss in the PM.

It can be seen from the references that the eddy current loss in the PM of the PMB is caused by the gyro effect and the PM’s non-uniform magnetization. Because the gyro effect caused by the magnetic bearing rotor radial periodic offset is controlled by the radial AMB, for the axial PMB structure used in this paper, only the eddy current loss caused by the non-uniform magnetization of the PM is considered.

The PMB structure scheme is shown in Figure 3a. Table 1 shows the design parameters for the PMB. In order to analyze the loss using the electromagnetic field theory, the polar coordinates are established. The bearing configuration schematic is shown in Figure 3b, where ${\mathsf{\Sigma }}_1$ and ${\mathsf{\Sigma }}_2$ represent the surfaces of the inner and outer diameters of the stator, respectively, and ${\mathsf{\Sigma }}_3$ and ${\mathsf{\Sigma }}_4$ represent the surfaces of the inner and outer diameters of the rotor, respectively.

Using the relationship between magnetization M and current density J, the relationship between current density J and magnetic vector potential A, the relationship between magnetic vector potential A and electric field intensity E, and the relationship between electric field intensity E and eddy current loss, we can obtain the eddy current loss of the rotor PM:

(2)${\mathrm{P}}_{eddy-stator}=\sigma {\displaystyle \sum _{i=1}^{10}{\displaystyle {\int }_{V_R}d{\tau }_S\left[{\left(e_{11}+e_{21}+e_{31}\right)}^2+{\left(e_{12}+e_{22}+e_{32}\right)}^2\right]}}$

where $\sigma$ is the conductivity of the PM; ${d\tau }_S=R_{S}d{\mathsf{\Phi }}_{S}d{z}_{S}d{R}_S$ is the micro-element of the stator volume; and $e_{11}$, $e_{21}$, $e_{31}$, $e_{12}$, $e_{22}$ and $e_{32}$ represent the micro-elements of the electric field intensity generated by the rotor magnetizing current in the stator PM:

(3)$\{\begin{array}{l}{e}_{11}=\frac{{\mu }_0\omega }{4\pi }\left[{\displaystyle {\int }_{\mathsf{\Sigma }3}{M}_i\cos {\mathsf{\Phi }}_3\cdot \frac{d{S}_3}{r_{3S}}+}{\displaystyle {\int }_{\mathsf{\Sigma }3}{M}_i\sin {\mathsf{\Phi }}_3\frac{\omega R_3{R}_S\sin \left({\mathsf{\Phi }}_S-{\mathsf{\Phi }}_3\right)}{r_{3S}^3}d{S}_3}\right]\\ e_{12}=-\frac{{\mu }_0\omega }{4\pi }{\displaystyle {\int }_{\mathsf{\Sigma }3}{M}_i\sin {\mathsf{\Phi }}_3}\frac{d{S}_3}{r_{3S}}\\ e_{21}=-\frac{{\mu }_0\omega }{4\pi }\left[{\displaystyle {\int }_{\mathsf{\Sigma }4}{M}_i\cos {\mathsf{\Phi }}_4\cdot \frac{d{S}_4}{r_{4S}}+}{\displaystyle {\int }_{\mathsf{\Sigma }4}{M}_i\sin {\mathsf{\Phi }}_4\frac{\omega R_4{R}_S\sin \left({\mathsf{\Phi }}_S-{\mathsf{\Phi }}_4\right)}{r_{4S}^3}d{S}_4}\right]\\ e_{22}=\frac{{\mu }_0\omega }{4\pi }{\displaystyle {\int }_{\mathsf{\Sigma }4}{M}_i\sin {\mathsf{\Phi }}_4}\frac{d{S}_4}{r_{4S}}\\ e_{31}=-\frac{{\mu }_0\omega }{4\pi }{\displaystyle {\int }_{V_R}{M}_i\cos {\mathsf{\Phi }}_R\cdot \frac{d{\tau }_R}{r_{RS}^2}}\\ e_{32}=\frac{{\mu }_0\omega }{4\pi }({\displaystyle {\int }_{V_R}{M}_i\sin {\mathsf{\Phi }}_R\cdot \frac{d{\tau }_R}{r_{RS}^2}}-{\displaystyle {\int }_{V_R}{M}_i\cos {\mathsf{\Phi }}_R\cdot \frac{2R_S{R}_R\sin \left({\mathsf{\Phi }}_S-{\mathsf{\Phi }}_R\right)}{r_{RS}^4}}d{\tau }_R)\end{array}$

where $M_i$ represents the $i\mathrm{th}$ harmonic amplitude of the magnetization of the rotor PM; $r_{\alpha S}$ is the distance between the stator and surface ${\mathsf{\Sigma }}_{\alpha }(\alpha =3,4)$, the micro-element of the surface of the inner and outer diameters of the rotor $d{S}_{\alpha }=R_{\alpha }d{\mathsf{\Phi }}_{\alpha }d{Z}_{\alpha }(\alpha =3,4)$, the micro-element of the rotor volume ${d\tau }_R=R_{R}d{\mathsf{\Phi }}_{R}d{Z}_{R}d{R}_R$, and the rotor angular velocity $\omega =2\pi \times n/60$.

A similar expression can be written if we replace the subscripts as follows:

$\begin{array}{l}{V}_S\to V_R\\ V_R\to V_S\\ r_{3S}\to r_{1R}\\ r_{4S}\to r_{2R}\\ {\mathsf{\Sigma }}_4\to {\mathsf{\Sigma }}_2\\ {\mathsf{\Sigma }}_3\to {\mathsf{\Sigma }}_1\\ \omega \to -\omega \end{array}$

Through the experiment, the axial magnetization, measured in the circumferential direction, of the PM in the PMB is shown in Figure 4.

Through the FFT transformation, the magnetization magnitude of the harmonic of the stator and rotor magnets is obtained (Figure 5).

Bringing the results into Equations (2) and (3), the eddy current loss of the rotor is P_rotor-eddy = 8.463 mW, and the eddy current loss of the stator is P_stator-eddy = 6.982 mW.

4. Losses Estimation of the AMB

4.1. Stator Winding Losses

The AMB’s losses consist of stator winding losses, hysteresis loss, eddy current loss, and the excess loss of the radial magnetic bearing. The scheme of AMB is shown in Figure 6a. The entity stator of the AMB is shown in Figure 6b.

The gravitational system, composed of the AMB stator and rotor, is essentially unstable. In order to overcome this nonlinear relationship between the electromagnetic force and current, in this paper, an AMB biased PM is used to linearize it. The magnetic circuit produced by the PM is shown in Figure 6c. The bias magnetic field produced by the PM provides the main bearing capacity, whereas the control current applied to the windings on the stator X+, X−, Y+, Y− is only used as a magnetic control force. The magnetic circuits produced by X+, X− stator windings with the positive control current are shown in Figure 6d,e.

We can see that, in the same direction as the control current, due to the different direction of the winding, the excitation of the magnetic field enhances the side of the air gap while weakening the other side of the air gap flux. Table 2 shows the design parameters for the AMB.

The AMB’s winding loss can be calculated by:

(4)$p_{Cu\_coil}=n_{mb}R{i}^2=n_{mb}\frac{N{p}_e}{A_e}{i}^2{\rho }_{cu}\left[1+\delta (T_1-20)\right]$

where n is the number of windings, R is the winding resistance, i is the winding current, N is the number of coil winding turns, p_e is the coil winding perimeter of an AMB’s pole, A_e is the effective cross-section area of coil winding, ρ_cu is the resistivity of the coil material (ρ_cu = 0.01851 Ω·mm²/m), δ is temperature coefficient of copper (δ = 0.004/°C), T₁ is the test temperature. Because of the PM biased AMB, no biased current exists in the windings. In the MLRW, winding losses and currents of the AMB are calculated by Equation (4) and shown in Table 3, with a rotation speed of 5000 r/min and extreme vacuum of 1.8 × 10⁻⁷ Pa. The total winding loss of the AMBs is 0.074 W at 20 °C. The winding loss is affected by temperature, and the resultant curves are shown in Figure 7.

4.2. Iron Loss Estimation of AMB

In this study, the iron losses in the AMB can be divided into the eddy current loss, the hysteresis loss, and the excess loss. All of these are affected by the core material, the magnetic field changes in frequency, the magnetic flux density, and the volume of the core. The iron losses of the AMB occur in the stator and rotor.

The stator core and rotor core use lamination material 1J50, which is used to reduce losses in the eddy currents. The eddy current loss may appear in the case of a shaft position fluctuation. The eddy current loss can be estimated as:

(5)$P_{R\_eddy}=k_c{f}_R^2{B}_R^2{V}_R$

where k_c is the coefficient of eddy current loss; f_R is the operating frequency of the magnetic field changes in the frequency of the AMB, which is twice as large as the mechanical frequency because two flux density periods are contained in a geometric period of the shaft; B_R is the maximum flux density; V_R is the affected volume.

The hysteresis loss can be calculated according to the analytical equation:

(6)$P_{R\_hys}=k_h{f}_R{B}_m^{\alpha }{V}_R$

where k_h is the factor of hysteresis loss and α is the Steinmetz constant.

The excess loss can be expressed as:

(7)$P_{R\_ex}=k_e{f}_R^{1.5}{B}_R^{1.5}{V}_R$

where k_e is the coefficient of the anomalous eddy current loss.

The combined estimation equation for the iron core losses can be also presented in the form:

(8)$P_{R\_core}=\left(k_h{f}_R{B}_R^{\alpha }+k_c{f}_R^2{B}_R^2+k_e{f}_R^{1.5}{B}_R^{1.5}\right)V_R$

The iron loss depends on the loss factor of the iron core material, the magnetic field changes in the frequency, the flux density amplitude, and the volume of the core. In this study, the AMB was made of a silicon steel sheet with a thickness of 0.2 mm. According to the loss curves of the silicon steel sheets with a sinusoid supply (Figure 8), the loss factors and Steinmetz constant calculated by the direct fitting computed method are shown in Table 4.

The main parameters and loss calculation results are shown in Table 5, which shows that the total loss is 0.786 W, and the percentages of different loss types are shown in Table 6. The eddy current and hysteresis losses are the major losses, whereas the excess loss is the smallest loss (0.13% of the total iron loss), and can be ignored in the iron loss. The hysteresis loss, classical eddy current loss, and excess eddy current loss versus magnetic magnitude, while the AMB is working with a rotation speed of 5000 r/min, are shown in Figure 9.

In addition to the stator and rotor core of the AMB, the iron loss also exists in the magnet ring, the PM, and the return ring. As solid materials are used in the above elements, the iron loss calculation formula is slightly different from that in Equation (8). Ignoring the excess loss, the iron loss can be expressed as:

(9)$P_{S\_iron}=P_{S\_h}+P_{S\_e}=\left(k_{S\_h}+k_{S\_e}\right){\sigma }_S{\delta }_S^2{f}_S^2{B}_S^2{V}_S$

where k_{S_h} is the hysteresis loss factor, k_{S_e} is the eddy current loss factor, σ_S is the conductivity, δ_S is the thickness, B_S is the maximum flux density, f_S is the operating frequency of the AMB, V_S is the effective volume. The iron cobalt vanadium soft magnetic alloy was adopted in this prototype. The main parameters and loss calculation results are shown in Table 7.

5. Loss Estimation of the BLDCM

Space applications generally require low-power reaction wheel systems. Thus, a three-phase BLDCM with ironless and slotless stator was used in the MLRW for driving the rotor to rotate at a high speed [20,21]. The structure and parameters of the BLDCM are shown in Figure 10 and Table 8, respectively. The stator winding coil and Hall sensor are fixed on the stator frame constructed of polyimide. The rotor consists of the PM poles comprising the SmCo magnet, outer rotor core, and inner rotor core made of the 1J50 material (perm alloy with 50% nickel). The range of the speed is from −5000 r/min to 5000 r/min, which can output torques in both directions. Therefore, the loss type in the BLDCM consists of the iron core loss from the rotor and the coil copper loss from the stator.

5.1. Stator Winding Losses

Because the current passes through the windings and the polyimide stator frame, the stator loss of the BLDCM consists mainly of the copper loss. Considering the influence of the temperature on the armature windings, the coil copper loss can be expressed as:

(10)$P_{coil}^M=n_M{R}_c{i}_c^2=n_M\frac{N_c{p}_c}{A_c}{\rho }_{cu}\left[1+\delta \left(T_c-20\right)\right]i_c^2$

where n_M is number of the phases and T_c is the test temperature. R_c, i_c, N_c, p_c, A_c, ρ_cu and δ represent the same parameters as in Equation (4) for the BLDCM. The estimated results of the coil losses have a positive correlation with the temperature of the BLDCM, as shown in Figure 11. At room temperature, the total copper losses are 2.91 W.

5.2. Rotor Losses Estimation

In the BLDCM, the rotor loss consists of the inner rotor core loss, outer rotor core loss, and PM loss. Similarly to the AMB, the iron losses are separated into the hysteresis loss, eddy current loss, and excess loss. PM losses are caused by the induced eddy current.

The iron core is made up of thin laminations in order to reduce iron core loss. The iron core loss depends on the loss factor of iron core material, flux density amplitude, and magnetic field changes in frequency. From refs. [22,23], the iron loss of the 1J50 material consists of the eddy current loss, hysteresis loss, and excess loss, and can be calculated as:

(11)$\begin{array}{l}{P}_M=P_h^M+P_c^M+P_e^M\\ \{\begin{array}{l}{P}_h^M=K_{h}f{\displaystyle \sum _{k=0}^{\infty }k\left(B_{k\max }^{\alpha }+B_{k\min }^{\alpha }\right)}{V}_M\\ P_c^M=K_c{f}^2{\displaystyle \sum _{k=0}^{\infty }{k}^2\left(B_{k\max }^2+B_{k\min }^2\right)V_M}\\ P_e^M=\frac{K_e}{{\left(2\pi \right)}^{\frac{3}{2}}}\frac{1}{T}{\displaystyle {\int }_0^T\left({\left|\frac{d{B}_r\left(t\right)}{dt}\right|}^{1.5}+{\left|\frac{d{B}_{\theta }\left(t\right)}{dt}\right|}^{1.5}\right)}dt{V}_M\end{array}\end{array}$

where K_h, K_c, K_e and α are shown in Table 4 for the same material; B_kmax and B_kmin are the maximum and minimum values of the elliptical harmonic of each order, respectively [24]; B_r and B_θ are the radial and circumferential components, respectively. f is the operating frequency of the magnetic field, T = 60/(p $\times$ n) = 0.0015 s, and V_pm is the effective volume.

5.2.1. Loss Estimation of the Inner Rotor Core

When the MLRW works with the rotation speed of 5000 r/min, the maximum and the minimum FFT values of the harmonics in the inner rotor core are calculated (Figure 12).

According to Equation (11), since the amplitudes of the harmonics over an order of 17 are smaller and negligible, only the maximum 17-order harmonic magnetic flux density is considered. It can be calculated that the hysteresis loss and eddy current loss in the inner rotor core are 0.472 and 0.00168 W, respectively.

(12)$\{\begin{array}{l}{P}_{in\_h}^M=K_h{f}_M{\displaystyle \sum _{k=0}^{17}k\left(B_{k\max }^{\alpha }+B_{k\min }^{\alpha }\right)V_{in}^M}\\ P_{in\_c}^M=K_c{f}_M^2{\displaystyle \sum _{k=0}^{17}{k}^2\left(B_{k\max }^2+B_{k\min }^2\right)}{V}_{in}^M\end{array}$

where f_M = $p_M\times n_m/60$ Hz, $p_M$ is the pole pair number of the BLDCM, $n_M$ is the rated speed, $V_{in}^M$ is the effective volume of the inner rotor core.

The first derivative of the radial and axial flux density in the inner rotor core is calculated in Figure 13a,b, respectively.

According to the excess loss formula given in Equation (11) and the computed data in Figure 13a,b, it can be calculated that the excess loss in the inner rotor core is 0.0009 W.

(13)$P_{in\_e}^M=\frac{K_e}{{\left(2\pi \right)}^{\frac{3}{2}}}\frac{1}{T}{\displaystyle {\int }_0^T\left({\left|\frac{d{B}_r\left(t\right)}{dt}\right|}^{1.5}+{\left|\frac{d{B}_{\theta }\left(t\right)}{dt}\right|}^{1.5}\right)}dt{V}_{in}^M$

5.2.2. Loss Estimation of the Outer Rotor Core

Similarly, the maximum and the minimum FFT values of the harmonics in the outer rotor core are calculated (Figure 14).

It can be calculated that the hysteresis loss and eddy current loss in the outer rotor core are 0.626 and 0.00364 W, respectively.

(14)$\{\begin{array}{l}{P}_{ou\_h}^M=K_h{f}_m{\displaystyle \sum _{k=0}^{17}k\left(B_{k\max }^{\alpha }+B_{k\min }^{\alpha }\right)}{V}_{ou}^M\\ P_{ou\_c}^M=K_c{f}_m^2{\displaystyle \sum _{k=0}^{17}{k}^2\left(B_{k\max }^2+B_{k\min }^2\right)}{V}_{ou}^M\end{array}$

where $V_{ou}^M$ is the effective volume of the outer rotor core.

The first derivative of the radial and axial flux density in the outer rotor core is calculated in Figure 15a,b, respectively.

According to the excess loss formula given in Equation (11) and the computed data in Figure 15a,b, it can be calculated that the excess loss in the inner rotor core is 0.00207 W.

(15)$P_{ou\_e}^M=\frac{K_e}{{\left(2\pi \right)}^{\frac{3}{2}}}\frac{1}{T}{\displaystyle {\int }_0^T\left({\left|\frac{d{B}_r\left(t\right)}{dt}\right|}^{1.5}+{\left|\frac{d{B}_{\theta }\left(t\right)}{dt}\right|}^{1.5}\right)}dt{V}_{ou}^M$

5.2.3. Loss Estimation of the PM

The calculation results of the harmonic amplitude in the PM of the BLDCM at the rated speed are shown in Figure 16.

There is no hysteresis loss in the PM, and the eddy current loss produced by the harmonic flux density can be expressed as:

(16)$P_{ie}=\frac{1}{T}{\displaystyle {\int }_0^T{\displaystyle {\int }_0^{\frac{d_p}{2}}{\displaystyle {\int }_0^L{\displaystyle {\int }_{-\frac{{\tau }_i}{2}}^{\frac{{\tau }_i}{2}}\frac{J_i^2\left(x,t\right)}{{\sigma }_i}dxdydzdt}}}}$

Since the other higher harmonics magnetic flux density are small and negligible, only the first 17 orders of the harmonics are calculated in the eddy current loss.

(17)$P_{PM}^M={\displaystyle \sum _{n=1}^{17}\frac{{\alpha }_{PM}{\sigma }_{PM}{L}_{PM}{d}_{PM}{\tau }_n^3{B}_n^2{\omega }_{en}^2}{16}\left(\frac{{\delta }_n}{3{\tau }_n}\frac{\sinh \frac{{\delta }_n}{{\tau }_n}-\sin \frac{{\delta }_n}{{\tau }_n}}{cosh\frac{{\delta }_n}{{\tau }_n}-\cos \frac{{\delta }_n}{{\tau }_n}}\right)}$

where τ_n is the polar distance of the space harmonics (${\tau }_n=\frac{\pi D_r}{2np}$, D_r is the inner diameter of the PM), ω_en is the electrical angular velocity, ${\sigma }_{PM}$ = 1.17 × 10⁶ S/m is the conductivity of the PM, ${\alpha }_{PM}$ = 0.67 is the pole-arc coefficient. Substituting parameters and values into Equation (17), the result of the eddy current loss of the PM of the rotor is $P_{PM}^M$ = 3.42 W.

6. Thermal Field Analysis and Measurement for the Prototype of MLRW

6.1. Calculation for Heat Generation Rate

The various losses in the MLRW can eventually be converted to a temperature increase. As mentioned above, these losses consist of the copper losses and iron core losses in the AMB and the BLDCM, and they are regarded in terms of the heat generation rate per unit volume. The models of energy transfer, including the heat conduction and thermal radiation, play a major role in cooling MLRW when working in a highly vacuum state.

The heat generation rates of the copper loss of the coils in the MBs and the BLDCM are calculated by Equation (18) and the results are shown in Table 9. The total copper loss of the MLRW is 2.98 W.

(18)$q_{loss}=\frac{P_{loss}}{V_{EV}}$

The heat generation rates of the iron core loss in the MBs and the BLDCM are calculated by Equation (18) and shown in Table 10. The loss percentages of the loss types and the different parts are shown in Table 11 and Table 12, respectively. The total loss, including the copper loss and the iron core loss, is 8.75 W. The total iron core loss of the MLRW is 5.77 W. The estimated results show the BLDCM loss is the main loss, accounting for 84.98% of the total loss in the MLRW.

6.2. Convection Heat Transfer Coefficient

The thermal analysis of the system conforms to the law of conservation of energy, namely, for any closed system:

(19)$Q-W=\Delta U+\Delta KE+\Delta PE$

where Q is the heat in the system; W is the external work; $\Delta U$ is the change in the internal energy of the system; $\Delta KE$ is the kinetic energy of the system; $\Delta PE$ is the potential energy of the system. For most engineering heat conduction problems, $\Delta KE=\Delta PE=0$. For the flywheel thermal analysis studied in this paper, when the flywheel is working in a stable situation, the heat generated by the flywheel is equal to the heat dissipated into the space; that is, $Q-W=0$.

Because of the high-vacuum working situation, the effect of heat convection caused by air can be ignored. The heat transfer formula of the heat conduction due to the contact between different temperature objects is:

(20)$\raisebox{1ex}{$Q_i$}\!\left/ \!\raisebox{-1ex}{$t$}\right.=KA\left(T_{hot}-T_{cold}\right)/d$

where $Q_i$ is the heat transferred in time $t$; $K$ is the heat transfer coefficient; $A$ is the effective area between the two contact objects; $T$ is the temperature of the object; and $d$ is the distance between the two planes.

The heat transfer formula of heat radiation is:

(21)${\mathsf{\Phi }}_{1,2}=\frac{{\sigma }_b\left(T_1^4-T_2^4\right)}{\frac{1-{\epsilon }_1}{{\epsilon }_1{A}_1}+\frac{1}{A_1{X}_{1,2}}+\frac{1-{\epsilon }_2}{{\epsilon }_2{A}_2}}$

where ${\mathsf{\Phi }}_{1,2}$ is the heat flux between surfaces 1 and 2; ${\sigma }_b=5.67\times {10}^{-8} \mathrm{W}/\left({\mathrm{m}}^2{\times \mathrm{K}}^4\right)$ is the radiation constant for the blackbody; ${\epsilon }_1$, ${\epsilon }_2$ are the emissivities of the radiation of surfaces 1 and 2; A₁, A₂ are the effective areas of radiating faces 1 and 2, respectively; T₁, T₂ are the surface temperatures of the radiation surfaces 1 and 2, respectively; $X_{1,2}$ is the coefficient between the radiating surfaces 1 and 2 (representing the ratio of the total radiation emitted by surface 2 to the amount of the radiation absorbed by surface 1).

6.3. Thermal Field Analysis of the MLRW

The thermal behavior of the MLRW depends on its cooling capability and losses in the system. The 3-D FEM of the MLRW was built based on the ANSYS software. Its mesh is shown as Figure 17a, and the total number of nodes was 426,034. The ambient temperature was set to 20 °C. Conduction and radiation are the heat transfer modes for the internal cooling of MLRW. The heat convection can be ignored on account of the high-vacuum environment. Due to the tight assembly between the components, the heat transfer mode is mainly heat conduction. Although the rotor is supported by magnetic levitation in the vacuum environment, the cooling of the rotor is only by means of thermal radiation. Based on steady-state thermal analysis, thermal parameters such as thermal conductivity, specific heat capacity, emissivity, and loss-load per unit volume of each component are assigned to the FEM of thermal distribution. The estimated temperature result is shown in Figure 17b. (In order to facilitate observation, a portion of the lower temperature components are hidden). The estimated maximum temperature is 54.2 °C, which is located at the BLDCM stator.

6.4. Thermal Optimization Design Based on Thermal Network Model

When the heat dissipation method and the heat transfer mode of the system are analyzed by the thermal network method, the components of the MLRW are equivalent to the individual node units. Each node contains the size of its own heat generation value and various heat-related parameters, such as the thermal conductivity and specific heat capacity. The heat transfer path between the components is expressed by the thermal conductivity or thermal resistance between the node units. The equivalent node unit and heat dissipation path is shown in Figure 18.

Then, the heat transfer relationship of each node unit can be expressed as the equivalent thermal network model, as shown in Figure 19.

The heat change formula on any node is:

(22)$c_i\frac{d{T}_i}{dt}={\displaystyle \sum }_j{E}_{i,j}\left(T_j^4-T_i^4\right)+{\displaystyle \sum }_j{D}_{i,j}\left(T_j-T_i\right)+q_i$

where $c_i$ represents the specific heat capacity of the node i; $T_i$ represents the temperature of node i; $q_i$ is the amount of heat generated internally by the corresponding node i per unit time; $E_{i,j}$ is the thermal radiation coefficient between node i and node j; $D_{i,j}$ is the thermal conductivity between node i and node j.

According to the Equations (20) and (21), the thermal conductance formula caused by the radiation and conductivity is:

(23)$\{\begin{array}{l}{E}_{i,j}=\frac{{\sigma }_b}{\frac{1-{\epsilon }_i}{{\epsilon }_i{A}_i}+\frac{1}{A_i{X}_{i,j}}+\frac{1-{\epsilon }_j}{{\epsilon }_j{A}_j}}\\ D_{i,j}=\frac{{\lambda }_i{A}_{i,j}}{{\delta }_{i,j}}\end{array}$

where ${\delta }_{i,j}$ is the distance of the thermal conductivity and ${\lambda }_i$ is the conductivity coefficient of node i.

According to the thermal field diagram of the MLRW shown in Figure 17b, the temperature of the PM and stator windings of the BLDCM correspond to node i₇ and i₁₀ in the equivalent thermal network model. The thermal parameters, such as the thermal radiation coefficient, thermal conductivity, and heat rate of different components in the MLRW for the thermal distribution, are based on the steady-state thermal analysis [25,26]. Therefore, the optimization target is to reduce the temperature of node i₇ and i₁₀. Then, in the thermal network equations of the whole system, there are:

(24)$\{\begin{array}{l}{E}_{7,8}\left(T_8^4-T_7^4\right)+E_{7,10}\left(T_{10}^4-T_7^4\right)+E_{7,13}\left(T_{13}^4-T_7^4\right)+D_{7,14}\left(T_{14}-T_7\right)+q_7=0\\ ⋮\\ E_{7,10}\left(T_{10}^4-T_7^4\right)+D_{9,10}\left(T_{10}-T_9\right)+D_{10,13}\left(T_{13}-T_{10}\right)+q_{10}=0\end{array}$

By calculation and comparison of the values of the thermal conductance of each heat transfer path, it was found that the thermal conductance of node i₁₄ is too small, which affects the heat dissipation efficiency of nodes i₇ and i₁₀. Thus, the polyimide of the fixed BLDCM stator windings was replaced by 1060 alloy, which has better thermal conductivity and hear emissivity. The results of the 3-D FEM verification are shown in Figure 20, and the maximum temperature located at the BLDCM stator is reduced from 54.2 to 34.8 °C.

7. Experimental Test

The prototype of the MLRW with the rated speed of 5000 r/min was manufactured in this study (Figure 21a). The MLRW was measured to confirm the loss estimation and the thermal field analysis. Four temperature measurement points were selected in the prototype, as shown in Figure 21b, where the thermistors were located at the PMB stator, base, BLDCM stator, and AMB stator.

The temperature rise curve when the MLRW reaches the rated rotation speed of 5000 r/min is shown in Figure 21c. It is shown that, when the MLRW thermal distribution is stable, the highest temperature is located at the end of the BLDCM stator windings, and the value is 36.8 °C. The maximum temperatures of the AMB stator and PMB stator are 30.1 and 24.5 °C respectively. The maximum temperature increase in the MLRW satisfies the safety margin of the magnetic and metal materials. The maximum temperature of the base is 24.5 °C. The errors in the calculated and measured values are shown in Table 13, with the maximum error of less than 10%. Thus, the experimental results verified the accuracy of the MLRW loss estimation and thermal field analysis.

8. Conclusions

In this article, the temperature increase due to the copper and iron losses in a MLRW was calculated. The copper and iron core losses in the PMB, AMB, and BLDCM were all predicted by their analytical equations. The thermal field of the prototype was analyzed by the 3-D FEM based on the loss values and the heat generation rates. The estimated maximum temperature was found to be located at the BLDCM stator, and the temperature was 34.8 °C. A prototype was fabricated, and the maximum measured temperature, also located at the BLDCM stator, was 36.8 °C. The maximum error between the calculated and measured temperatures was 5.4%, which verifies the loss estimation model and the thermal field analysis and optimization.

Author Contributions

Methodology and Writing—original draft, Z.H.; Project administration, T.W.; Supervision, X.L.; Validation, Y.S. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the National Natural Science Foundation of China (Grant Nos. 51925501, 82027802, 62073010 and 61773038), and in part by the supported by Ningbo Natural Science Foundation (No. 2021J011).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Structure of the MLRW.

Figure 2. Force analysis and coordinate system of the MLRW.

Figure 3. PMB scheme. (a) The structure of the PMB. (b) Bearing configuration schematic.

Figure 4. Axial magnetization of the PM in the PMB.

Figure 5. The magnetization magnitude of the harmonic of the PMB.

Figure 6. Structure of the AMB. (a) Scheme of the AMB structure. (b) The stator of the AMB. (c) Scheme of the PM magnetic circuit. (d) Scheme of the magnetic circuit of X+ winding with positive excitation current (upper stator). (e) Scheme of the magnetic circuit of X− winding with positive excitation current (upper stator).

Figure 7. Winding loss versus temperature for the AMB.

Figure 8. Loss curves of the silicon sheet with a sinusoid supply.

Figure 9. Hysteresis loss curve, classical eddy current loss curve, and excess eddy current loss curve in the AMB versus magnetic magnitude.

Figure 10. Structure of the BLDCM.

Figure 11. Winding loss of the BLDCM affected by temperature.

Figure 12. The maximum and minimum FFT values of the harmonics in the inner rotor core.

Figure 13. The first derivative of the radial and axial flux density in the inner rotor core. (a) Radial. (b) Axial.

Figure 14. The maximum and minimum FFT values of the harmonics in the outer rotor core.

Figure 15. The first derivative of the radial and axial flux density in the outer rotor core. (a) Radial. (b) Axial.

Figure 16. The harmonic amplitude in the PM of the BLDCM.

Figure 17. Finite element analysis of the MLRW. (a) 3-D finite element mesh of the MLRW. (b) The predicted thermal field of the MLRW by FEM.

Figure 18. Equivalent node unit and heat dissipation path of MLRW.

Figure 19. Equivalent thermal network model of the MLRW.

Figure 20. Optimized thermal field of MLRW by FEM.

Figure 21. The testing setup of the MLRW prototype. (a) Testing setup. (b) Positions of the four thermistors. (c) Measurement data.

Table 1

Design parameters of the PMB.

Parameters and Characteristics	Design Value of PMB
PMB mass, kg	2.68
Nominal air gap, mm	1.0
Outer diameter of the stator, mm	92.5
Inner diameter of the stator, mm	87.5
Outer diameter of the rotor, mm	98.5
Inner diameter of the rotor, mm	93.5
Bearing length, mm	10.0

Table 2

Design parameters of the AMB.

Parameters and Characteristics	Design Value of PMB
AMB mass, kg	1.15
Nominal air gap, mm	0.80
Maximum current in each coil, A	1.13
Number of winding turns	150
Outer diameter of the stator, mm	90.4
Inner diameter of the stator, mm	34.0
Bearing length, mm	25.6
Maximum force per direction, N	408.8
Current stiffness in center position, N/A	393.4
Position stiffness in center position, N/μm	4.18

Table 3

Copper loss calculation for the AMB (20 °C).

Coils	Measured Current	Resistance	Copper Loss
coil_1~8	0.052 A	3.4 Ω	0.0092 W
Total loss		0.074 W

Table 4

Direct fitting results of parameters for silicon steel sheets.

kh	A	kc	ke
73.0987	1.6	0.120388	1.48188 × 10−3

Table 5

Main parameters and loss calculation results of the AMB.

Iron Core	fR (Hz)	BR (T)	VR (mm3)	Iron Loss (W)
Stator core	166.7	1.4	17,780	0.488
Rotor core			10,857	0.298
Total loss				0.786

Table 6

Percentage of different iron loss types.

Iron Core	fR (Hz)	BR (T)	VR (mm3)	Iron Loss (W)
Value	0.598 W	0.187 W	1.12 × 10−3 W	0.816 W
Percentage	76.08%	23.79%	0.13%	1

Table 7

Other parameters and loss calculation results of the AMB.

Parts	Frequency	Volume	Conductivity	Max Flux Density	Iron Loss
Magnet ring	166.7 Hz	6993 mm3	2.5 × 106 S/m	0.85 T	0.159 W
PM	166.7 Hz	4995 mm3	1.1 × 106 S/m	0.85 T	0.0256 W
Return ring	166.7 Hz	18,152 mm3	2.5 × 106 S/m	0.85 T	0.258 W
	Total loss				0.443

Table 8

Main parameters of the BLDCM.

Parameters	Value
Number of stator slots	48
Number of pole pairs	8
Axial length (LPM), mm	14
Radial thickness of the PM (dPM), mm	3
Phase resistance (room temperature), Ω	0.11
Sensor	Hall

Table 9

Heat rate of the copper losses in the MLRW.

Coils	Copper Loss	Volume	Heat Rate
AMB Coil_1~8	0.0092 W	3461 mm3	2601 W/m3
BLDCM Coil_1~48	0.0606 W	794 mm3	76,826 W/m3
	Total loss		2.98 W

Table 10

Heat rate of the iron loss in the MLRW.

Parts		Iron Loss	Volume	Heat Rate
BLDCM	Outer core	0.632 W	30,657 mm3	20,615 W/m3
	Inner core	0.474 W	28,289 mm3	16,755 W/m3
	PM	3.42 W	14,602 mm3	234,214 W/m3
AMB	Stator core	0.488 W	17,780 mm3	27,447 W/m3
	Rotor core	0.298 W	10,857 mm3	27,448 W/m3
	Magnetic ring	0.159 W	6993.2 mm3	22,736 W/m3
	Return ring	0.258 W	18,152 mm3	14,213 W/m3
	PM	0.0256 W	4995.1 mm3	5125 W/m3
PMB	Stator PM	0.0070 W	28,274 mm3	247.6 W/m3
	Rotor PM	0.0085 W	30,159 mm3	281.8 W/m3

Table 11

Loss percentage of the loss types.

Losses	Copper Loss	Iron Loss	Total Loss
Value	2.98 W	5.77 W	8.75 W
percentage	34.1%	65.9%	100%

Table 12

Loss percentage of the different parts.

Losses	BLDCM Loss	AMB Loss	PMB Loss	Total Loss
Value	7.436 W	1.303 W	0.0155 W	8.75 W
percentage	84.98%	14.89%	0.17%	100%

Table 13

Error between calculated and measured temperatures.

Test Points	Point 1	Point 2	Point 3	Point 4
Calculated value	24.8	24.8	34.8	31.1
Measured Value	24.5	24.5	36.8	30.1
Error	1.2%	1.2%	5.4%	3.3%

References

1. Zhang, J.; Yuan, L.; Chen, S.L.; Liang, Y.; Yang, G. A survey on design of reaction spheres and associated speed and orientation measurement technologies. ISA Trans.; 2019; 99, pp. 417-431. [DOI: https://dx.doi.org/10.1016/j.isatra.2019.09.014] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31570163]

2. Tang, J.; Fang, J.; Wen, W. Superconducting Magnetic Bearings and Active Magnetic Bearings in Attitude Control and Energy Storage Flywheel for Spacecraft. IEEE Trans. Appl. Supercond.; 2012; 22, 5702109. [DOI: https://dx.doi.org/10.1109/TASC.2012.2218245]

3. Zhang, W.; Zhu, H.; Yang, H.; Chen, T. Experimental Analysis and Full Prediction Model of a 5-DOF Motorized Spindle. Energies; 2017; 10, 75. [DOI: https://dx.doi.org/10.3390/en10010075]

4. Mao, C.; Zhu, C. Vibration Control for Active Magnetic Bearing Rotor System of High-Speed Flywheel Energy Storage System in a Wide Range of Speed. Proceedings of the 2016 IEEE Vehicle Power and Propulsion Conference (VPPC); Hangzhou, China, 17–20 October 2016; pp. 1-6.

5. Zhao, Y.; Lu, F.; Fan, C.; Yang, J. A Novel Axial-Flux Dual-Stator Toothless Permanent Magnet Machine for Flywheel Energy Storage. Symmetry; 2022; 14, 61. [DOI: https://dx.doi.org/10.3390/sym14010061]

6. Yu, Y.; Zhang, W.; Sun, Y.; Xu, P. Basic Characteristics and Design of a Novel Hybrid Magnetic Bearing for Wind Turbines. Energies; 2016; 9, 905. [DOI: https://dx.doi.org/10.3390/en9110905]

7. Ji, W.; Ni, F.; Gao, D.; Luo, S.; Lv, Q.; Lv, D. Electromagnetic Design of High-Power and High-Speed Permanent Magnet Synchronous Motor Considering Loss Characteristics. Energies; 2021; 14, 3622. [DOI: https://dx.doi.org/10.3390/en14123622]

8. Sun, J.; Fang, J. A novel structure of permanent-magnet-biased radial hybrid magnetic bearing. J. Magn. Magn. Mater.; 2011; 323, pp. 202-208.

9. Wang, H.; Wu, J.F.; Li, Y.; Wu, Y.H. Thermal design of attitude control flywheel system for small satellites. Opt. Precis. Eng.; 2015; 23, 2265. [DOI: https://dx.doi.org/10.3788/OPE.20152308.2265]

10. Tir, M.A.; Mirimani, S.M.; Marignetti, F. A novel structure of passive magnetic bearing with axial magnetization. Proceedings of the 5th Annual International Power Electronics, Drive Systems and Technologies Conference (PEDSTC 2014); Tehran, Iran, 5–6 February 2014; pp. 311-316.

11. Tănase, N.; Morega, A.M.; Nedelcu, A.; Ilie, C. Radial-axial passive magnetic bearing system—Numerical simulation aided design solutions. Proceedings of the 2015 9th International Symposium on Advanced Topics in Electrical Engineering (ATEE); Bucharest, Romania, 7–9 May 2015; pp. 566-572.

12. Liu, Q.; Zhao, M.S.; Han, B.C.; Fan, Y.H.; Sun, J.J.; Zheng, S.Q. Energy optimization and experimental for a permanent magnet-biased redial magnetic bearing. Opt. Precis. Eng.; 2019; 27, pp. 2420-2428.

13. Achari, V.; Jagathy, R.; Madhu, G.; Bindu, G.R. Design and analysis of zero cogging Brushless DC motor for spacecraft applications. Proceedings of the 2010 ECTI International Confernce on Electrical Engineering/Electronics, Computer, Telecommunications and Information Technology; Chaing Mai, Thailand, 19–21 May 2010; pp. 254-258.

14. Chebak, A.; Viarouge, P.; Cros, J. Analytical Computation of the Full Load Magnetic Losses in the Soft Magnetic Composite Stator of High-Speed Slotless Permanent Magnet Machines. IEEE Trans. Magn.; 2009; 45, pp. 952-955. [DOI: https://dx.doi.org/10.1109/TMAG.2009.2012478]

15. Romanenko, A.; Smirnov, R.; Jastrzebski, P.; Pyrhönen, O. Losses estimation and modelling in active magnetic bearings. Proceedings of the 2014 16th European Conference on Power Electronics and Applications; Lappeenranta, Finland, 26–28 August 2014; pp. 1-8.

16. Su, W.M.; Mao, S.H.; Wang, P.J. Core losses estimation in design of high speed electric machines. Proceedings of the 2016 Progress in Electromagnetic Research Symposium (PIERS); Shanghai, China, 8 August–11 September 2016; pp. 320-323.

17. Boglietti, A.; Cavagnino, A.; Staton, D.; Shanel, M.; Mueller, M.; Mejuto, C. Evolution and Modern Approaches for Thermal Analysis of Electrical Machines. IEEE Trans. Ind. Electron.; 2009; 56, pp. 871-882. [DOI: https://dx.doi.org/10.1109/TIE.2008.2011622]

18. Huang, Z.; Fang, J.; Liu, X.; Han, B. Loss Calculation and Thermal Analysis of Rotors Supported by Active Magnetic Bearings for High-Speed Permanent-Magnet Electrical Machines. IEEE Trans. Ind. Electron.; 2016; 63, pp. 2027-2035. [DOI: https://dx.doi.org/10.1109/TIE.2015.2500188]

19. Bertotti, G. General properties of power losses in soft ferromagnetic materials. IEEE Trans. Magn.; 1988; 24, pp. 621-630. [DOI: https://dx.doi.org/10.1109/20.43994]

20. Fang, J.; Li, W.; Li, H. Self-Compensation of the Commutation Angle Based on DC-Link Current for High-Speed Brushless DC Motors with Low Inductance. IEEE Trans. Power Electron.; 2014; 29, pp. 428-439. [DOI: https://dx.doi.org/10.1109/TPEL.2013.2254499]

21. Han, B.; Zheng, S.; Wang, Z.; Le, Y. Design, Modeling, Fabrication, Test of A Large Scale Single-Gimbal Magnetically Suspended Control Moment Gyro. IEEE Trans. Ind. Electron.; 2015; 62, pp. 7424-7435. [DOI: https://dx.doi.org/10.1109/TIE.2015.2459052]

22. Binesti, D.; Ducreux, J.P. Core Losses and Efficiency of Electrical Motors Using New Magnetic Materials. IEEE Trans. Magn.; 1996; 32, pp. 4887-4889. [DOI: https://dx.doi.org/10.1109/20.539278]

23. Roshen, W. Iron Loss Model for Permanent Magnet Synchronous Motors. IEEE Trans. Magn.; 2007; 43, pp. 3428-3434. [DOI: https://dx.doi.org/10.1109/TMAG.2007.899687]

24. Kong, X.; Wang, F.; Xu, Y.; Xing, J. Analysis and calculation of iron losses of high-speed permanent magnet machines. Electr. Mach. Control.; 2010; 14, pp. 26-30.

25. Zhang, S.; Han, B.; Hong, L.; Fang, J. Thermal and structural coupling analysis of magnetically suspended flywheel rotor. Proceedings of the Sixth International Symposium on Instrumentation and Control Technology; Beijing, China, 13–15 October 2006.

26. Zhai, L.; Sun, J.; Ma, X.; Han, W.; Luo, X. Thermal–structure coupling analysis and multi-objective optimization of motor rotor in MSPMSM. Chin. J. Aeronaut.; 2019; 32, pp. 1733-1747. [DOI: https://dx.doi.org/10.1016/j.cja.2018.09.008]