**Space Vector Modulation Strategies for Self-Sensing Three-Phase Radial Active Magnetic Bearings †**

## **Dominik Wimmer \*, Markus Hutterer, Matthias Hofer and Manfred Schrödl**

Institute of Energy Systems and Electrical Drives, Technische Universität Wien, Gußhausstr. 25, 1040 Vienna, Austria; markus.hutterer@tuwien.ac.at (M.H.); matthias.hofer@tuwien.ac.at (M.H.); manfred.schroedl@tuwien.ac.at (M.S.)


Received: 19 April 2019; Accepted: 10 May 2019; Published: 14 May 2019

**Abstract:** The focus of this study lies on the investigation of the space vector modulation of a self-sensing three-phase radial active magnetic bearing. The determination of the rotor position information is performed by a current slope-based inductance measurement of the actuator coils. Therefore, a special pulse width modulation sequence is applied to the actuator coils by a conventional three-phase inverter. The choice of the modulation type is not unique and provides degrees of freedom for different modulation patterns, which are described in this work. For a self-sensing operation of the bearing, certain constraints of the space vector modulation must be considered. The approach of a variable space vector modulation is investigated to ensure sufficient dynamic in the current control as well as the suitability for a self-sensing operation with an accurate rotor position acquisition. Therefore, different space vector modulation strategies are considered in theory as well as proven in experiments on a radial magnetic bearing prototype. Finally, the performance of the self-sensing space vector modulation method is verified by an external position measurement system.

**Keywords:** active magnetic bearing; self-sensing; radial three-phase bearing; space vector modulation

## **1. Introduction**

Magnetic bearings are of great significance for the stabilization of levitating rotors. Due to the fact that it is impossible to stabilize all degrees of freedom of a rigid body by permanent magnets, a force of a different physical origin must be applied for rotor stabilization [1,2]. Therefore, the rotor can be stabilized by electromagnets, which is stated as active magnetic bearing (AMB). AMBs require a position feedback information of the rotor to allow a stable operation. In this work, a self-sensing method is used to obtain the rotor position instead of using separate position sensors. The operation of self-sensing AMBs has been a field of research for many years [3–5] and provides advantages concerning sensor failure, construction space and production costs of the AMB. The self-sensing position determination of this study is performed by a modulation-based current switching ripple evaluation of the actuator coils. Previous studies presented different approaches for extracting the rotor position information out of the current ripple, such as current slope measurements [6,7], current ripple demodulation [8] or by the usage of artificial neural networks [9]. In this proposal, the so-called INFORM (Indirect Flux Detection by Online Reactance Measurement) method is used to determine the rotor position. This method was originally designed for the rotor angle determination of a permanent magnet synchronous motor [10]. Concerning AMBs, the INFORM method is based on a current slope measurement, detecting the inductance change depending on the rotor's eccentricity. As previously described in [11], the implementation of the self-sensing operation was based on an injection of voltage

pulses to the actuator coils. Therefore, the current controller stops in equidistant time steps and measurement pulses are applied to the coils. This approach has the drawback of a limited bandwidth of the position measurement and an interruption of the current controller. To avoid this circumstance, the required pulse pattern for position measurement is embedded in the pulse width modulation (PWM) sequence of the current controller. The use of an embedded pulse pattern, in particular the 3-Active pattern [12], is the point of origin for the space vector modulation in this work.

## **2. Self-Sensing Bearing Setup**

The bearing setup consists of two radial homopolar six pole AMBs with a common shaft as illustrated in Figure 1. The bias flux of the bearing is realized by the use of permanent magnets (PM). An axial displacement of the rotor is stabilized by a positive axial stiffness given by the bias flux and the geometry of the bearing (Figure 1b).

**Figure 1.** (**a**) Structural design of the six pole radial homopolar active magnetic bearing. (**b**) Cross section of the shaft: The bias flux (indicated by red lines) is generated by means of permanent magnets.

The stabilization of a radial rotor eccentricity is performed by the coils, which are driven in a differential configuration. The coils of the poles *U*+, *V*+, *W*+ and the opposite coils are connected in wye-configuration. For achieving low system costs, a conventional three-phase inverter is aspired for the control of the AMB [13,14]. Two opposite coils are corresponding to one phase, which enables the use of a three-phase inverter. Figure 2 shows a differential transformer at the connection point of the positive and negative phase of the bearing.

**Figure 2.** Differential current slope measurement by means of a differential open loop transformer.

By the usage of an open loop transformer in each phase, it is possible to connect the AMB to a three-phase inverter. The transformer can be realized as a separate unit or also be integrated in the printed circuit board (PCB) of the power electronics [15]. The transformer output provides the differential current slope signal, which is required for the self-sensing operation of the bearing. Figure 3a shows the current ripple caused by the PWM switching pattern. Although the stator and the rotor components were built from laminated iron sheets, it can be seen that the current slope is

distorted by eddy currents. To overcome this problem, the current slope measurement is performed as a differential evaluation of two opposite coils to suppress the influence of eddy currents. The differential current slope evaluation is not limited to the homopolar design and can be also applied to heteropolar bearings [16]. Another approach would be a model -based consideration of the eddy currents as shown in [17]. Figure 3b shows the output voltage of the open loop transformer corresponding to the differential current Δ*IU* = *IU*<sup>+</sup> − *IU*<sup>−</sup> of Figure 3a. After the settling time, the output voltage is proportional to the differential current slope *<sup>d</sup> dt*Δ*IU*.

**Figure 3.** (**a**) Current ripple of the actuator coils. The influence of eddy currents is mostly suppressed in the differential current signal Δ*IU* (*fPWM* = 20 kHz, *UDC* = 60 V). (**b**) Output signal of the open loop transformer corresponding to (**a**). The output voltage is proportional to the differential current slope.

Thus, the rotor position can be obtained by the approximation

$$
\mu(t) = L(\mathbf{x}, y)\frac{d}{dt}\Delta I(t) \to L(\mathbf{x}, y) = u(t)\left(\frac{d}{dt}\Delta I(t)\right)^{-1} \tag{1}
$$

with the coil voltage *u*(*t*) and the position-depending inductance *L*(*x*, *y*), which is described in [12].

## **3. Problem Formulation**

Measurements on a prototype of a self-sensing radial AMB (Figure 4) have shown significant power losses in the bearing in the 3-Active PWM mode. Consequently, the power losses lead to a temperature rise in the bearing, which is undesirable in many applications.

**Figure 4.** Prototype of a self-sensing radial homopolar active magnetic bearing with six poles. The bias flux of the bearing is realized by the use of permanent magnets.

Beside the control current, the current ripple of the coils causes remarkable power losses in the AMB. Hence, major power losses in the prototype occur in the flux leading paths by means of iron losses in the laminated iron sheets. One possible solution for the reduction of the eddy current losses is the use of soft magnetic composites (SMC), which is described in [18]. However, SMC materials have drawbacks like a smaller permeability and mechanical limitations [19]. This work follows an approach, which is independent of the used material. Figure 5a shows the 3-Active SVM pattern with the corresponding power losses (Figure 5b) of the prototype as a function of the DC-link voltage *UDC*.

**Figure 5.** (**a**) 3-Active space vector modulation: Each PWM period contains voltage space vectors from each phase (*U***+**, *V***+**, *W***+**). (**b**) Power losses of the bearing due to the current ripple as a function of the DC-link voltage (3-Active SVM, *fPWM* = 20 kHz).

For a fixed switching frequency, it is obvious to decrease *UDC* for achieving small power losses. On the one hand, a high value of *UDC* causes a high current ripple (Figure 5b), which provides a high magnitude of the differential current slope information. On the other hand, the iron losses are increased and therefore, undesirable power losses occur in the flux leading paths. To keep the power losses in the bearing small, a low level of *UDC* is aspired. This circumstance gives the motivation to enhance the self-sensing position measurement for dealing with small current slopes. Furthermore, it must be considered that the dynamic of the current controller depends on *UDC*. The focus of this work lies on the investigation of different space vector modulations to ensure a sufficient current dynamic as well as a high quality of the position measurement. The design of SVM contains degrees of freedom, which can be used for the development of specific pulse patterns with especially high current dynamics or high quality of the position measurement. Taking this one step further, it is possible to use different kinds of SVM during the operation of the AMB. This leads to the variable SVM and enables the combination of the properties of different modulations.

## **4. Space Vector Modulation**

Figure 6 shows the symmetrical arrangement of the magnetic poles of the bearing in the xy-plane. The spatially distributed arrangement of the poles enables the definition of six fundamental voltage space vectors, which are aligned with the magnetic poles of the bearing. The voltage space vectors can be obtained by a conventional power inverter with three half bridges by the switching states shown in Table 1 [20].


**Table 1.** Voltage space vector definition of a three-phase inverter.

HS = High-Side switch, LS = Low-Side switch; 1 = closed, 0 = open.

**Figure 6.** Cross section of the six pole homopolar radial AMB. The fundamental voltage space vectors *U***+**, *U***−**, *V***+**, *V***−**, *W***+**, *W***−** are aligned with the poles.

The fundamental space vectors shape a hexagon, which defines the possible modulation range of the voltage space vectors. Each point in this hexagon can be reached by a linear combination of six fundamental space vectors (*U***+**, *U***−**, *V***+**, *V***−**, *W***+**, *W***−**) and two zero space vectors (*Z***+**, *Z***−**). The zero space vector occurs if either all high-side or low-side switches of the three-phase inverter are closed. The calculation between the reference coordinate system (*x*, *y*) and the three phase system (*U*, *V*, *W*) is done by the Clarke-transformation [21]. The degree of freedom in the composition of the linear combination of the space vectors is restricted by the following design rules for self-sensing operation.

## *4.1. Design Rules for SVM*

The design criteria of the SVM achieve a good quality of the position measurement and a high current controller bandwidth. Therefore, the design of the SVM underlies certain restrictions to allow a self-sensing operation of the magnetic bearing:


high dynamic of the current control. For a symmetrical (angle independent) operation of the current controller, it is beneficial to limit the modulation amplitude to the in-circle of the possible modulation area. Theoretically, the symmetrical modulation amplitude can achieve a value of *Rmax* = <sup>√</sup>3/<sup>2</sup> <sup>≈</sup> 0.866 for the given AMB system.

• **Inverter:** Short pulse lengths could cause problems in semiconductor switches. Hence, the specified recovery time of the switches must be considered in the PWM pattern [22]. Furthermore, a proper operation of a potential charge pump of the gate driver must be ensured. Therefore, the pulse pattern has to provide at least one switching action in each phase.

The design procedure of the SVM is based on a representation of the desired voltage space vector *P*, which is specified by the superior current controller. Hence, the desired space vector *P* is given by means of all eight voltage space vectors.

$$\mathbf{P} = \boldsymbol{\varepsilon}\_{\mathbf{U}\_{+}} \mathbf{U}\_{\bullet} + \boldsymbol{\varepsilon}\_{\mathbf{U}\_{-}} \mathbf{U}\_{\bullet} + \boldsymbol{\varepsilon}\_{\mathbf{V}\_{+}} \mathbf{V}\_{\bullet} + \boldsymbol{\varepsilon}\_{\mathbf{V}\_{-}} \mathbf{V}\_{\bullet} + \boldsymbol{\varepsilon}\_{\mathbf{W}\_{+}} \mathbf{W}\_{\bullet} + \boldsymbol{\varepsilon}\_{\mathbf{W}\_{-}} \mathbf{W}\_{\bullet} + \boldsymbol{\varepsilon}\_{\mathbf{Z}\_{+}} \mathbf{Z}\_{\bullet} + \boldsymbol{\varepsilon}\_{\mathbf{Z}\_{-}} \mathbf{Z}\_{\bullet} \tag{2}$$

Therefore, the coefficients *cU*+, *cU*−, *cV*+, *cV*−, *cW*+, *cW*−, *cZ*+, *cZ*<sup>−</sup> define the length of the corresponding space vector. The length of the respective voltage space vector is related to the duration of the voltage pulse in the PWM pattern. The following introduced variants of SVM differ substantially in the number of the active space vectors within one PWM period. In this context, "active" refers to the fundamental space vectors with an amplitude unlike zero. In the following considerations the minimum pulse duration *tINF* is normalized to the PWM period.

## *4.2. 6-Active SVM*

The 6-Active SVM uses all fundamental space vectors for the formation of *P*. Hence, *P* is built by a linear combination of the adjoining fundamental voltage space vectors. The remaining time of the PWM period is equally distributed to all fundamental voltage space vectors to build a combined zero space vector. The maximum modulation amplitude is obtained, if the length of one active space vector drops under *tINF* and violates the timing requirements of a differential current slope measurement. Although the admissible modulation area is defined by the solid hexagon in Figure 7, the intended modulation area is limited by the minimum and maximum symmetrical modulation amplitude (*Rmin*, *Rmax*).

**Figure 7.** 6-Active SVM: The space vector *P* is formed by six fundamental space vectors (*tINF* = 0.1).

The 6-Active SVM allows six independent current slope measurements within a PWM period. For this reason, the 6-Active SVM is well suited for the self-sensing method. However, this kind of

modulation has the drawback of a very limited modulation amplitude *Rmax*. The impact of the small value of *Rmax* is caused by the minimal pulse width *tINF* by means of Equations (3) and (4).

$$R\_{\text{max}} = \frac{\sqrt{3}}{2} (1 - 6 \,\, t\_{\text{INF}}), \,\, t\_{\text{INF}} < \frac{1}{6} \tag{3}$$

$$R\_{\min} = \mathcal{O} \tag{4}$$

For achieving a higher value of *Rmax*, the number of space vectors is reduced in the following considerations.

## *4.3. 3-Active Low Dynamic Range SVM*

In contrast to the 6-Active SVM, the 3-Active Low Dynamic Range (LDR) SVM uses either three positive (*U***+**, *V***+**, *W***+**) or three negative fundamental space vectors (*U***−**, *V***−**, *W***−**).

**Figure 8.** 3-Active LDR SVM: The space vector *P* is formed by the positive fundamental space vectors (*U***+**, *V***+**, *W***+**). (*tINF* = 0.1).

Figure 8 shows a 3-Active LDR SVM with a linear combination of *U***+**, *V***+**, *W***+**. The corresponding coefficients from Equation (2) can be calculated by the inner product of *P* and the respective fundamental voltage space vector (Equations (5)–(7)).

$$\mathbf{c}\_{IL\_{\bullet}} = \frac{1}{3} + \frac{2}{3} \mathbf{P} \cdot \mathbf{U}\_{\bullet} \tag{5}$$

$$\mathcal{L}\_{V\_{+}} = \frac{1}{3} + \frac{2}{3} \text{ P} \cdot \text{V}\_{\text{+}} \tag{6}$$

$$\mathcal{L}\_{\mathcal{W}\_{+}} = \frac{1}{3} + \frac{2}{3} \text{ } \mathbf{P} \cdot \mathbf{W}\_{\*} \tag{7}$$

By using only three fundamental space vectors, the maximum modulation amplitude is limited to 0.5 and causes a smaller drop of *Rmax*

$$R\_{\max} = \frac{1}{2}(1 - 3 \text{ } t\_{INF}), \ t\_{INF} < \frac{1}{3} \tag{8}$$

$$\mathcal{R}\_{min} = \mathcal{O} \tag{9}$$

by an increase of *tINF* than the 6-Active SVM.

## *4.4. 4-Active SVM*

The 4-Active SVM is based on the 3-Active LDR SVM, but enhances the modulation amplitude by means of an additional fundamental voltage space vector, which is located next to the desired space vector *P* like shown in Figure 9. For an effective implementation, the desired space vector *P* is built by a linear combination of the adjoining fundamental voltage space vectors. The remaining time of the PWM period is distributed equally to the positive (*U***+**, *V***+**, *W***+**) or negative (*U***−**, *V***−**, *W***−**) fundamental voltage space vectors.

**Figure 9.** 4-Active SVM: The space vector *P* is formed by the negative (*U***−**, *V***−**, *W***−**) fundamental space vectors and *V***<sup>+</sup>** for an enhancement of the modulation amplitude (*tINF* = 0.1).

The limits of the symmetrical modulation amplitude

$$R\_{\text{max}} = \frac{\sqrt{3}}{2} (1 - 3 \, t\_{INF})\_\prime \, \, t\_{INF} < \frac{1}{5} \tag{10}$$

$$R\_{\min} = \mathbf{O} \tag{11}$$

are obtained if the length of one vector falls below *tINF*. In contrast to the 3-Active LDR SVM, the maximum modulation amplitude is enhanced by a factor of <sup>√</sup>3 (Equation (10)).

## *4.5. 3-Active High Dynamic Range SVM*

The 3-Active High Dynamic Range (HDR) SVM is designed for a maximum modulation amplitude using three fundamental space vectors and one of the zero space vectors (*Z***+**, *Z***−**). The desired space vector *P* is mainly formed by the two adjoining fundamental space vectors (*V***+**, *W***−** in Figure 10).

**Figure 10.** 3-Active HDR SVM: The desired space vector *P* is built with three fundamental space vectors (involving a space vector from each phase) and one zero space vector (*tINF* = 0.1).

In order to get a current slope information by a voltage space vector from all space axis, the third fundamental space vector (*U***<sup>+</sup>** in Figure 10) is added with the minimum length *tINF*. A zero space vector fills the remaining time of the pulse pattern and ensures the required switching action in each phase. Equations (12) and (13) show the possible modulation range

$$R\_{\text{max}} = \frac{\sqrt{3}}{2} (1 - t\_{\text{INF}}), \ t\_{\text{INF}} < \frac{1}{5} \tag{12}$$

$$R\_{\min} = 2\sqrt{3} \, t\_{INF} \tag{13}$$

for a symmetrical operation.

## *4.6. Combination of the SVMs*

The design of the SVMs shows, that each SVM has individual characteristics concerning the modulation amplitude and the usability for self-sensing operation. For an optimal operation of the bearing, it is possible to combine the properties of different SVMs. Figure 11 shows a comparison of the symmetrical modulation amplitudes as a function of *tINF*.

**Figure 11.** Comparison of the symmetrical modulation amplitude.

The 3-Active HDR has the property of a lower boundary *Rmin* of the modulation amplitude. It is obvious, that a low value of *tINF* leads to a high maximum modulation amplitude. The 3-Active LDR can operate up to *tINF* < 1/3 , which can be advantageous for applications with a high switching frequency. Concerning the performance of the self-sensing operation, each fundamental space vector gives additional information for the rotor position.

Figure 12 shows a combination of multiple SVMs for *tINF* = 0.1. The choice of the SVM is determined in a way, that the desired space vector *P* is built by the SVM, which provides the highest number of active space vectors within a PWM period. To avoid nonessential switching between different SVMs, a hysteresis can be defined by means of an overlapping modulation area.

**Figure 12.** Combination of different SVMs during operation for achieving a high modulation area and a maximum performance of the self-sensing position measurement (*tINF* = 0.1).

## *4.7. SVM Switchover*

It is possible to make a distinction between two scenarios of SVM switchover. The first scenario is a sector switchover within a SVM. In the 4-Active and 3-Active HDR SVM the fundamental voltage space vectors used depend on the sector of *P*. Thus, the fundamental space vector changes if *P* changes to a new sector of the modulation area. Therefore, a switchover of the sector causes a change in the current ripple profile. Although the steady state of the mean value of the current is not affected by this effect, the phase currents obtain a transient error by a sector switchover. The amplitude of the current error is in the range of the amplitude of the current ripple. Figure 13a shows a simulation of a sector switchover for the 4-Active SVM for a space vector with ∥*P*∥<sup>2</sup> = 0 and *ϕ* = *π*/6. It can be seen that the mean values of the phase currents differ after the sector switchover, which result in a transient deviation that decays over time.

**Figure 13.** (**a**) Transient simulation of a sector switchover within the 4-Active SVM. The voltage space vector has an amplitude of zero but changes in the angle *ϕ* = *π*/6 + Δ*ϕ* to force a sector switchover. (**b**) Transient simulation of a switchover between the 3-Active and 6-Active SVM. Both SVMs represent a voltage space vector with zero amplitude (*fPWM* = 20 kHz, *tINF* = 0.14, *UDC* = 48 V, *δ* = 0.001 rad).

The second scenario is given by a switchover between different SVMs. Figure 13b shows a switchover between the 3-Active LDR and the 6-Active SVM for the voltage space vector ∥*P*∥<sup>2</sup> = 0 and *ϕ* = 0. Although both SVM represent a zero space vector, a drift of the mean values of the currents can be observed after the switchover. This effect is caused by the different current waveforms of the SVMs. In many applications, the current ripple is significantly smaller than the control current and the drift due to sector switchover is negligible. However, if the application requires a precise current control beneath the amplitude of the current ripple, a compensation strategy is required. One conceivable solution would be the introduction of a modified PWM cycle to suppress the current drift after a SVM switchover.

## **5. Measurements**

The following measurements compare the behavior of the designed SVMs, regarding power losses, dynamic of the current controller as well as the quality of the self-sensing position measurement. The measurements were performed on the prototype presented in Figure 4 applying the test setup of Figure 14. External eddy current-based position sensors were used as a reference for position measurements.

**Figure 14.** Symbolic arrangement of the test setup with external position sensors according to Figure 4.

The test setup with two homopolar radial magnetic bearings was controlled by independent three-phase inverters. Basically, a decoupled control of the rotor gives many degrees of freedom for advanced control [23]. However, simple decentralized PIDT1 position controllers provided sufficient performance for the following considerations. Concerning position control, it was assumed that the force on the rotor is proportional to the phase currents. In general, the currents of opposite coils (e.g., *IU*+, *IU*−) are not equal due to different inductances of an eccentrically levitating rotor. Therefore, the implemented force control by means of a static current force characteristic [24] is an approximation, but it does not cause any restrictions for the subsequent measurements.

## *5.1. Power Losses of the SVM Variants*

The initial aim was to decrease the power losses in the bearing by a reduction of the DC-link voltage. Figure 15 shows a comparison of the power losses in the bearing for different SVMs at 20 kHz switching frequency. There is no significant difference between the SVMs, which allows an almost power neutral switchover between different SVMs.

**Figure 15.** Comparison of the power losses in the bearing (**a**) and the current ripple of a phase (**b**) as a function of the DC-link voltage (*fPWM* = 20 kHz, *tinf* = 0.14).

## *5.2. Dynamic of the Current Controller*

The maximum modulation amplitude *Rmax* of the particular SVM has a significant impact on the dynamic of the current controller. Figure 16a shows a dynamic comparison of the current controller by means of the step response. Due to the fact that the 3-Active HDR SVM is not able to represent a modulation amplitude smaller than *Rmin* (Figure 11), it is combined with the 4-Active SVM to obtain a steady state without oscillation.

**Figure 16.** (**a**) Step response of the current controller for different SVMs. (**b**) Corresponding magnitude of the desired voltage space vector *P* (*fPWM* = 20 kHz, *UDC* = 12 V, *tinf* = 0.14).

Figure 16b indicates the amplitude of the desired voltage space vector *P*, which corresponds to the step response. It can be seen, that the 3-Active HDR SVM has the highest modulation amplitude. Thus, the controller is only saturated for a short period of time at *Rmax*. Concerning Figure 16, the 6-Active SVM is not able to inject the desired current of 7 A to the coil. The reason is the coil resistance and the connector cable. In this case, the 6-Active modulation can only be used for small currents and shows the use case for switching to a SVM with a higher modulation amplitude.

## *5.3. Noise of the Self-Sensing Method*

The noise of the self-sensing method is of great significance for a precise control of the bearing. Figure 17a shows a noise comparison of the different SVMs with *UDC* in the range of 20 V to 60 V. The 6-Active SVM has the lowest noise level, which is caused by the high number of voltage space vectors within a PWM period. Basically, the noise level increases proportional to <sup>∝</sup> *UDC*<sup>−</sup><sup>1</sup> . However, it is possible to shift the noise level in a certain range by an adaption of the analog circuit. Therefore, a low noise level can be achieved even at low values of *UDC*. Thus, Figure 17b shows a similar noise characteristic like Figure 17a at the half DC-link voltage level (e.g., *σ<sup>x</sup>* < 0.15 μm at 19 V) by means of an adaption of the analog circuit.

**Figure 17.** (**a**) Noise comparison of the self-sensing position measurement. (**b**) Reduction of the noise at low levels of *UDC* by an adaption of the analog circuit (*fPWM* = 20 kHz, rotor fixed in the center of the bearing, standard deviation *σ<sup>x</sup>* over 3000 samples of the x-position at each setpoint of *UDC*).

## *5.4. Linearity of the Self-Sensing Method*

Beside the noise level, the self-sensing method must provide a high linearity for a proper operation of the bearing. The AMB prototype has an air gap of 800 μm and the auxiliary bearing allows a radial operational range of 400 μm. Figure 18 shows a linearity analysis of the 6-Active SVM with external position sensors for different rotor setpoints.

**Figure 18.** Linearity analysis with external position sensors for different setpoints of the rotor position (*fPWM* = 20 kHz, 6-Active SVM).

The external sensors measure the rotor position on an aluminum disc, which is mounted close to the magnetic bearing. Due to an axial and angular shift of the sensors with regard to the AMB, the position information of the sensors is transformed in the self-sensing coordinate system. Figure 19 shows the self-sensing linearity error *e* = √ *ex*<sup>2</sup> + *ey*<sup>2</sup> for different rotor positions for a SVM with 3 and 6 active voltage space vectors. It can be seen that the SVMs have a similar distribution of the linearity error with a maximum linearity error of about 16 μm.

**Figure 19.** Analysis of the self-sensing linearity error *e* = <sup>√</sup>*ex*<sup>2</sup> <sup>+</sup> *ey*<sup>2</sup> as a function of the radial rotor displacement for different SVMs: (**a**) 3-Active LDR, (**b**) 6-Active.

The linearity of the self-sensing method is important for a proper control of the bearing, especially for robustness considerations of self-sensing magnetic bearings [25,26].

## *5.5. Small Signal Behavior*

Finally, a comparison of the small signal behavior of the self-sensing method was performed. For this purpose, a disturbance is applied to the self-sensing levitating rotor and the sensor signal is shown for comparison (Figure 20). Both signals show a consistent behavior, which manifests a precise self-sensing operation.

**Figure 20.** Comparison of the small signal behavior between an external sensor and the self-sensing method. A disturbance is applied at *t* = 0s(**a**) x-position (**b**) y-position (*fPWM* = 20 kHz, 6-Active SVM).

## **6. Results**

The measurements on the prototype demonstrated the characteristics of different space vector modulations concerning self-sensing and current dynamic. The combination of different kinds of modulation allows a high modulation amplitude of voltage space vectors. This feature enables the possibility of a reduction of the DC-link voltage, while still having sufficient dynamic in the current control. Measurements on a prototype showed, that AMB losses due to the current ripple have an approximate square dependence of the DC-link voltage. Therefore, already a small reduction of the DC-link voltage can reduce power losses in the AMB. Furthermore, a measurement of the power losses revealed that the power losses caused by the switching ripple are nearly independent of the used SVM, which enables almost power neutral SVM switchover. The combination of different SVMs leads to a high dynamic of the current controller, which was proven by the step response of the system. An analysis of the switchover characteristics between different SVMs indicated, that a SVM switchover has an impact to the phase currents in the range of the ripple amplitude. This effect may lead to disturbances in the current control and could be suppressed by a compensation cycle in the SVM after a switchover. Regarding the self-sensing position measurement, the experimental results of the prototype showed a low noise level, which varied with the used SVM. A linearity analysis with external sensors showed an overall error of about 16 μm for different setpoints in the admissible rotor orbit.

## **7. Conclusions and Outlook**

The focus of this study lied on an improvement of the space vector modulation of a three-phase self-sensing radial magnetic bearing. Different kinds of space vector modulations for self-sensing operation were presented in theory as well as proven in experiments on a prototype of a radial active magnet bearing. The results showed that the dynamic of the current controller and the quality of the self-sensing position control can be adjusted by the choice of the space vector modulation. The self-sensing method showed a high linearity and low noise of the rotor position for low control currents. In a next step, saturation effects of the flux leading paths will be considered regarding nonlinearities of the self-sensing operation to analyze the limitations of the self-sensing control. Further investigations will be performed at various rotational speeds to reveal potential speed-depending effects of the self-sensing control, especially with respect to a robust control of the system.

**Author Contributions:** Conceptualization, D.W., M.H. (Markus Hutterer), M.S.; methodology, D.W., M.H. (Markus Hutterer); software, D.W.; validation, D.W., M.H. (Markus Hutterer); formal analysis, D.W., M.H. (Markus Hutterer), M.H. (Matthias Hofer); investigation, D.W.; resources, M.S.; writing—original draft preparation, D.W.; writing—review and editing, M.H. (Markus Hutterer), M.H. (Matthias Hofer), M.S.; supervision, M.S.;

**Funding:** This research received no external funding.

**Acknowledgments:** This project was supported by Technische Universität Wien in the framework of the Open Access Publishing Program.

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

## **Abbreviations**

The following abbreviations are used in this manuscript:


## **References**


© 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/).

## *Article* **Thermal Behavior of a Magnetically Levitated Spindle for Fatigue Testing of Fiber Reinforced Plastic**

## **Daniel Franz \*, Maximilian Schneider, Michael Richter and Stephan Rinderknecht**

Institute for Mechatronic Systems in Mechanical Engineering, Technische Universität Darmstadt, 64287 Darmstadt, Germany; schneider@ims.tu-darmstadt.de (M.S.); richter@ims.tu-darmstadt.de (M.R.); rinderknecht@ims.tu-darmstadt.de (S.R.)


Received: 3 April 2019; Accepted: 30 April 2019; Published: 3 May 2019

**Abstract:** This article discusses the critical thermal behavior of a magnetically levitated spindle for fatigue testing of cylinders made of fiber reinforced plastic. These cylinders represent the outer-rotor of a kinetic energy storage. The system operates under vacuum conditions. Hence, even small power losses in the rotor can lead to a high rotor temperature. To find the most effective way to keep the rotor temperature under a critical limit in the existing system, first, transient electromagnetic finite element simulations are evaluated for the active magnetic bearings and the electric machine. Using these simulations, the power losses of the active components in the rotor can be derived. Second, a finite element simulation characterizes the thermal behavior of the rotor. Using the power losses calculated in the electromagnetic simulation, the thermal simulation provides the temperature of the rotor. These results are compared with measurements from an experimental spindle. One effective way to reduce rotational losses without major changes in the hardware is to reduce the bias current of the magnetic bearings. Since this also changes the characteristics of the magnetic bearings, the dynamic behavior of the rotor is also considered.

**Keywords:** active magnetic bearings; kinetic energy storage; fiber reinforced plastic; fatigue testing; thermal behavior

## **1. Introduction**

Flywheels store energy as kinetic energy of the rotor and can provide a cost efficient solution for short-term energy storage and load smoothening services in electricity grids (e.g., [1,2]). Their advantages lie in the high possible number of cycles and the low initial cost per unit of power. The main drawbacks are their high stand-by-losses and the relatively low energy density compared to other storage technologies. In order to utilize the advantages, the energy density of one system should be increased while decreasing the stand-by losses at the same time.

## *1.1. Outer-Rotor Fywheel Design*

One possible flywheel design is an outer-rotor setup (e.g., [3–5]). It promises high energy densities due to the large radii of the rotor and high rotational speeds. A realized full-scale system is described in [6]. A model of the system is shown in Figure 1a. The rotor is a magnetically levitated hollow cylinder. The outside of the rotor is made out of fiber reinforced plastic (FRP) with a circumferential fiber orientation, using carbon fibers. The inside of the rotor consists of different rotor parts of the active and passive components. Active magnetic bearings (AMBs) are used for radial levitation and

a passive magnetic bearing is used for axial levitation. A permanent magnet synchronous machine (PMSM) accelerates and decelerates the rotor. All rotating components of the magnetic bearings and the PMSM are integrated on the inner circumference of the FRP rotor. Under rotation, these components press against the FRP, resulting in radial compressive stress transversal to the fiber orientation, which is superimposed by circumferential stress in the fiber direction (see Figure 1b). The energy density of the system increases with the radii of the rotor and its rotational speed, both factors lead to an increased stress in the FRP [7]. Consequently, high stress in the material is necessary to reach high energy densities. Charging and discharging the flywheel leads to cyclically varying mechanical stresses. To investigate the cyclic stability and lifetime of the FRP rotor, cyclic material tests, such as cyclic transverse tensile tests (derived from [8]), cyclic transverse compressive tests (derived from [9]), and cyclic four-point bending tests are performed on material samples. These samples are thin walled and relatively easy to fabricate with a high quality, whereas the rotor of the flywheel is much thicker and harder to produce in an industrial fiber winding process. Furthermore, the thermal expansion of carbon fiber and the epoxy plastic matrix used differ, leading to inner stress when a thick walled FRP structure cools down from the curing temperature. Hence, there is a much higher probability of imperfections and defects in the thick walled FRP structure of an outer-rotor flywheel than in the thin material samples. Therefore, the applicability of test results obtained with thin material samples to the rotor of the flywheel has to be investigated by tests with thick walled rotors. A specialized test rig was designed and set up to perform cyclic tests on thick-walled FRP specimens, which represent the outer-rotor of kinetic energy storage. These specimens are smaller, cheaper and have a lower energy content, reducing the danger during destructive testing.

**Figure 1.** (**a**) Model the outer-rotor Flywheel described in [6] with a halved rotor; (**b**) Top view on the rotor showing the segments pressing against the fiber reinforced plastic (FRP) under rotation.

## *1.2. Testing Procedure*

The ratio between the circumferential and the radial stress in a test specimen should be as close as possible to the one present in a full-scale flywheel. In the system described in [6], the radial transverse compressive stress in the FRP is 60 MPa and the circumferential longitudinal stress is 389 MPa at its maximum speed of 15,000 rpm. In the design of the flywheel, the maximum stress was limited to half of the expected strength of the FRP, to account for the high data uncertainty. To test the limits of the FRP, the stress in the specimen is doubled compared to the flywheel. To create this state of stress in the test specimen, a circumferentially segmented steel ring is placed inside the FRP cylinder of the specimen, and rotated at 30,000 rpm. This results in a transverse compressive stress of 100 MPa and a longitudinal tension of 775 MPa in the FRP. With an outer diameter of 190 mm, the surface speed of the specimen reaches approximately 298 m/s. The fatigue test will be conducted by alternating between a maximum speed of 30,000 rpm and a minimum speed of 15,000 rpm. At minimum speed the transversal compressive stress is 26 MPa and the longitudinal tension is 200 MPa. Because of the well-known high tensile strength of carbon fiber, the chance of fiber fracture is low at this state of stress. The probability of matrix fracture is expected to be much higher [10]. It is planned to perform up to 200,000 cycles per specimen. With a cycle time of 30 s this will take about 70 days. Additionally, overload tests are planned, where the specimen is accelerated to 40,000 rpm, resulting in a surface speed of nearly 398 m/s, a circumferential longitudinal tension of 1,400 MPa and a radial compressive stress of 180 MPa. The latter should lead to failure of the matrix. In the longitudinal direction, the FRP can withstand static tensions over 2,000 MPa. To reduce air drag and the consequential heating of the specimen, all tests are performed under vacuum conditions, with a pressure below 0.05 Pa.

## *1.3. Test Rig Description*

Figure 2a shows a section view of the designed test rig, which was introduced in [11]. A hub and a shaft coupling connected the specimen to the driving spindle. This configuration, rather than an outer-rotor setup, was chosen to protect the active parts of the spindle in case of a failure of the specimen at high speed. Each of the eight segments of the steel ring inside the FRP weighed 1.027 kg, and the distance between their center of gravity (CoG) and the rotation axis is about 55 mm. Hence, at 40,000 rpm one segment contained around 30 kJ of kinetic energy. The surrounding containment was designed to absorb this energy in case of specimen failure. The design criteria for the containment were derived from [12]. The containment also serves as a vacuum chamber. A detailed view of the spindle is shown in Figure 2b. The motor in the middle is a water cooled PMSM with four poles, a maximum torque of 9.6 Nm, and a maximum power of 30 kW. In order to avoid excessive wear of the high-speed drive, the rotor is supported by AMBs. The radial AMBs were designed in a heteropolar configuration with laminated cores on rotor and stator, both made out of 0.2 mm thick sheets of non-oriented silicon steel (NO20). The radial position of the rotor is measured with four eddy current sensors at each bearing, two for each radial direction.

**Figure 2.** (**a**) Cross section of the test rig; (**b**) Cross section of the spindle.

To lift the combined mass of the spindle rotor and the specimen of roughly 20 kg, the axial AMB needs a large pole area, which leads to a big outer diameter of the thrust disk. Because of the high stress at high speed, a monolithic design was chosen (see [13]). The stator of the axial AMB was fabricated of the soft magnetic composite Siron S400b of the PMG Füssen GmbH [14]. An inductive sensor at the upper end of the rotor detects the axial position. Two-rowed hybrid spindle bearings with a static load rating of 20 kN are used as backup bearings. The rotor temperature is measured with two infrared sensors—one between the upper radial AMB and the thrust disk and one next to the lower radial position sensors, below the lower AMB. The moment of inertia of the specimen is about ten times higher than that of the rotor, and the surface speed of the specimen is more than two times higher. To protect the rotor and the spindle in case of a specimen failure, a predetermined breaking point was included in the hub (see Figure 3). This tapering also uncouples the dynamics of the rotor from the specimen to some extent.

**Figure 3.** (**a**) Cross section of the specimen; (**b**) Assembled specimens.

One goal of the test rig operation is to minimize the total test duration for one specimen, maximizing the number of specimens that can be tested in a given amount of time. The test rig should thus be operated with as much power as possible. High power normally leads to high losses, as well as subsequent heating. Losses on the stator can be effectively cooled with water, but without a medium for convection or thermal conduction through physical contact, the rotor can only transfer heat to the stator via radiation. Hence, even small power losses in the rotor can lead to a high rotor temperature. The temperature of the rotor is crucial for the feasibility of the cyclic fatigue testing, since the temperature of the specimen affects the strength of the FRP and can change the test result. Furthermore, the magnets of the PMSM should not be operated over a certain temperature. Thus, the overall goal of this paper is to minimize the test duration under the boundary condition of a maximum rotor and specimen temperature. For this purpose, the losses in the rotor during operation and their influence on the rotor temperature have to be identified. Furthermore, feasible methods of minimizing the most crucial losses in an existing test rig are derived.

## **2. Method**

With the goal of finding the shortest total test duration, this paper proposes a structured approach to evaluate the heating of a magnetically levitated rotor in a vacuum, and deduces measures to avoid critical temperatures. An overview on the approach is shown in Figure 4. First, losses in the rotor of both radial AMBs, the axial AMB, and the PMSM were calculated for different rotational speeds, using transient electromagnetic 2D finite element (FE) models implemented in ANSYS Maxwell 2015. The highest and lowest operating speeds of the spindle followed from the demanded stress in the FRP, as described in Section 1.2. For each component, a polynomial was fitted to the related losses, to obtain

a mean loss value, which was used for the calculation of the steady state temperature. Next, air friction losses were calculated analytically, based on the kinetic theory of gases. Finally, to calculate the rotor temperature, a 3D-FE-model of a quarter of the test rig was set up and evaluated using ANSYS Workbench 19.1. Because of the water cooling, losses on the stator were not considered; instead the stator temperature in the simulation was predefined with measured data. A transient evaluation of the model was compared to measured temperatures on the test rig, to adjust thermal and loss calculation. To compute the rotor temperature after a long time of cycling, a steady state evaluation of the adjusted model was used. From this model, the most influential components on the rotor temperature can be derived. The focus was on components that allow for an improvement of the thermal behavior without major changes in the hardware. One option, which will be discussed as an example, is the reduction of the bias current in the radial AMBs [15]. Since this also changed the behavior of the rotor in the AMBs, a brief analysis of the rotor dynamics was performed using measurements and a mechanical 3D-FE-model of the rotor implemented in ANSYS Workbench 19.1. Another loss reduction can be achieved by reducing the control activity of the AMB [15]. Most of these strategies, like unbalance compensation (e.g., [16]) or a linear-quadratic-Gaussian-control (e.g., [17]), were not feasible at this point due to limitations of the AMB controller hardware. If no further loss reduction can be achieved in the components, the next step is to adjust the cycle time. Since rotational losses increase with speed and losses in the PMSM also increase with acceleration, the cycle time has a major influence on the mean losses and therefore on the rotor temperature. If the rotor temperature is below the critical values, the cycle time can be reduced. If the critical values are exceeded, the cycle time has to be increased. These options will be discussed further on.

**Figure 4.** Applied approach to minimize the test duration by calculating and reducing rotor losses for the existing test rig.

## **3. Modeling of the Thermal Stability of the Rotor**

The following sections discuss the loss calculation of the active magnetic bearings for the radial and axial direction and the electric machine, as well as losses caused by air friction. Subsequently, the thermal model focusing on the rotor temperature is described. For the model validation and potential adjustment, the calculation results are compared with measurements.

## *3.1. Loss Calculation of the Radial Active Magnetic Bearings*

Both radial AMBs have eight poles, which are evenly distributed on the inner circumference of the stator and arranged in the sequence N–S–S–N–N–S–S–N. Flux leakage between two poles with equal flux direction—that is, N–N or S–S—is small, and thus two neighboring poles with differing flux direction, called a pole pair, function as one independent electromagnet. One pole pair of the upper radial AMB with an idealized average path of the magnetic flux ϕ is shown in Figure 5a. The magnetic flux of the AMB passes through the NO20-metal sheets that surrounded the solid rotor, which is shown in the bottom left corner. The nominal air gap between the rotor and the surrounding stator in the centered position is 0.4 mm. To linearize the actuator characteristic, a differential winding design is used, where a common bias coil and a control coil, with associated currents, are used for two counteracting pole pairs [18]. In the test rig, two coils are mounted on each pole of the radial AMBs, the inner for bias current *IB* and the outer for the control current *IY*. Both coils have the same number of turns. All bias coils are connected in series for one AMB and the control coils for each direction in one AMB. The magnetic flux in a specific point in the rotor changes direction each time it passed a pole with a different flux direction. The flux also drops between two poles with equal flux direction. This remagnetization leads to hysteresis losses *Ph* and excess losses *Pe*, as well as eddy currents and subsequent losses *Pc*. To calculate the losses, a transient electromagnetic 2D-FE-model of the rotor was set up and evaluated. Figure 5b shows the simulated magnetic flux density B distribution in a quarter of the upper radial AMB at a rotational speed of 30,000 rpm. Only the bias current of *IB* = 5.67 A magnetized the rotor; the control current was set to zero. The flux density lay at 0.7 T on average, but locally exceeded the saturation flux density of the used NO20 sheets of 1.1 T.

**Figure 5.** (**a**) Quarter of the cross section of the upper radial active magnetic bearings (AMB); (**b**) Magnetic flux density in the upper radial AMB with only the bias current as excitation. The rotor turns clockwise.

For the loss calculation, the Bertotti formula [19] in the form of Equation (1) was applied, in which the three loss mechanisms *Pc*, *Ph*, and *Pe* depend on the amplitude of the flux density *Bm* and the frequency *f* with which the magnetic flux changes:

$$\begin{aligned} P\_B &= P\_\varepsilon + P\_h + P\_{\varepsilon\prime} \text{ with} \\ P\_\varepsilon &= \int k\_\varepsilon \left( f B\_m \right)^2 dV, \ P\_h = \int k\_h f B\_m^2 \, dV, \text{ and } P\_\varepsilon = \int k\_\varepsilon \left( f B\_m \right)^{1.5} \, dV. \end{aligned} \tag{1}$$

The corresponding loss constants *kc*, *kh*, and *ke* were derived by fitting data from the manufacturer of the electric steel obtained by standardized measurements [20] to Equation (1), which yields:

$$k\_{c,AMB} = 0.23 \, \frac{\text{Ws}^2}{\text{T}^2 \text{m}^3} \text{ and } k\_{h,AMB} = 193.6 \, \frac{\text{Ws}}{\text{T}^2 \text{m}^3}.$$

For the NO20 sheets used in the radial AMBs, *ke*,*AMB* could not be significantly determined and therefore was set to zero. The calculated total iron losses for the upper and lower radial AMB are shown in Table 1. Additionally, the switching of the amplifiers leads to further high frequency changes in the magnetic flux. The resulting rotor losses can be calculated with the same model and yielded 0.11 W for the upper AMB and 0.09 W for the lower AMB.

**Table 1.** Rotational losses in the radial AMBs.


## *3.2. Loss Calculation the of Axial Active Magnetic Bearing*

The magnetic flux in the axial AMB is mostly symmetrical to the rotational axis, so no significant remagnetization losses should occur due to rotation. Nevertheless, the switching of the amplifier leads to changes in the magnetic flux and associated rotor losses. Equation (1) is only applicable for thin sheets, but in the solid thrust disk of the axial AMB, without the usage of a laminated core or soft-magnetic composites, eddy current losses dominate. Therefore, the other loss mechanisms were neglected. To calculate eddy current losses a transient electromagnetic 2D-FE-model was analyzed. Figure 6a shows the magnetic flux density B in the axial AMB, with an air gap of 0.4 mm and a coil current of 3.15 A. The resulting force is 200 N, which is required to lift the rotor with the specimen, with a total weight of 20.3 kg. Changes in the magnetic flux in the rotor result in eddy currents and therefore in a non-zero current density *J* in the rotor. Eddy current losses *Pc*,*<sup>J</sup>* were calculated as the integral of the square of *J* over the rotor volume *V*, divided by the conductivity of the material σ*el* (see Equation (2)). The conductivity of X14CrMoS17 is <sup>σ</sup>*el*,*rot* = 1.43 <sup>×</sup> <sup>10</sup><sup>6</sup> <sup>1</sup> <sup>Ω</sup>*<sup>m</sup>* .

$$P\_{\mathfrak{c},l} = \frac{1}{\sigma\_{cl}} \int f^2 dV. \tag{2}$$

**Figure 6.** (**a**) Magnetic flux density in the axial active magnetic bearing (AMB); (**b**) Current density in the surface area of the thrust disk during switching of the axial AMB.

Since eddy currents oppose their provoking field, high frequency changes of the magnetic field only affect the outer layers of the conducting material. Hence, eddy currents are also limited to a thin region on the surface of the rotor (see Figure 6b). The losses of the axial AMB in the rotor due to the amplifier switching, derived by Equation (2), did not exceed 0.1 W.

## *3.3. Loss Calculation of the Permanent Magnet Synchronous Machine*

The last active component that induces losses in the rotor is the PMSM. Here, the magnetic field rotates synchronously with the rotor. Remagnetization of the rotor occurs because of flux drops at slots in the stator and non-harmonic changes of the motor current due to switching of the inverter. Losses were calculated via a transient electromagnetic 2D-FE model. In the solid permanent magnets, losses were calculated using Equation (2), with an electric conductivity of the magnets of <sup>σ</sup>*el*,*mag* = 1.1 <sup>×</sup> 106 <sup>1</sup> <sup>Ω</sup>*<sup>m</sup>* . The rotor is laminated underneath the permanent magnets to further reduce losses. Remagnetization losses in these sheets were calculated according to Equation (1). The corresponding coefficients are:

$$k\_{\rm c,PWM} = 0.12 \,\frac{\text{Ws}^2}{\text{T}^2 \text{m}^3}, \, k\_{\rm h,PWM} = 166.7 \,\frac{\text{Ws}}{\text{T}^2 \text{m}^3}, \, \text{and} \, k\_{\rm c,PWM} = 3.24 \,\frac{\text{Ws}^{1.5}}{\text{T}^{1.5} \text{m}^3}.$$

Results for different phase currents and rotational speeds are shown in Figure 7. The losses increased nearly linearly with the speed above 10,000 rpm. For root mean square values of the phase current (rms) below 10 A, losses were nearly independent of the phase current.

**Figure 7.** Rotor losses of the permanent magnet synchronous machine.

## *3.4. Air Friction Losses*

For the calculation of air friction losses, it has to be evaluated if the gas in the system can be modeled with continuum dynamics. This was done by calculating the mean free path *l* of the gas. *l* is the avarage distance a molecule moves before it collides with another molecule of the gas. For a gas with pressure *p*, temperature *T*, and a mean atom diameter *dm*, the mean free path *l* is calculated by Equation (3) (see [21]), where *<sup>k</sup>* <sup>=</sup> 1.381 <sup>×</sup> <sup>10</sup>−<sup>23</sup> <sup>J</sup> <sup>K</sup> is the Boltzmann constant,

$$
\vec{l} = \frac{kT}{\sqrt{2}\pi p d\_m^2}.\tag{3}
$$

The mean molecule diameter of air is *dm* = 3.559 <sup>×</sup> 10−<sup>10</sup> <sup>m</sup> [21]. During rotation, the maximum pressure in the test rig is *p* = 0.05 Pa and the minimum temperature *T* = 290 K, resulting in a minimum mean free path of *l* = 0.144 *m*. Since *l* is much bigger than most of the air gaps in the test rig, an air molecule is more likely to collide with the walls of the test rig than with other air molecules. Hence, continuum dynamics are not applicable for this system. Instead, the kinetic theory of gases was utilized, where molecules are modeled as randomly moving elastic spheres. In [22], air friction losses *Pa* were derived by means of the momentum excange beween a spinning rotor and air molecules. For an axial disc with an outer radius *ro* and an inner radius *ri* this yields [22]:

$$P\_{a, \text{ax}} = \pi p \sqrt{\frac{m\_W \pi}{2kT}} (r\_o^4 - r\_i^4) \alpha^2 \,, \tag{4}$$

for a radial cylindrical surface with a radius *rr* and a height *hr* [22]:

*Pa*,*rad* = 4π*p mW*π <sup>2</sup>*kT hrr* 3 *<sup>r</sup>*ω2, (5)

and for a truncated cone with a height *hc*, a lower radius *rl* and an upper radius *ru*:

$$P\_{a, \text{rad}} = \pi p \sqrt{\frac{m\_W \pi}{2kT}} \frac{\sqrt{\left(r\_l - r\_u\right)^2 + h\_c^2}}{r\_l - r\_u} \Big(r\_l^4 - r\_u^4\Big) \omega^2,\tag{6}$$

where ω is the rotational speed of the rotor and *mW* is the mean molecule mass of the gas. Hence, air friction losses increase linearly with *p*, quadratically with ω, and decrease with *T*. For dry air, the mean atom mass is approximately *mW* <sup>≈</sup> 4.81 <sup>×</sup> 10−<sup>26</sup> kg [21]. The losses were calculated for each surface of the rotor and the specimen and were subsequently superimposed. For the worst case, with a maximum pressure of *p* = 0.05 Pa and a minimum temperature of *T* = 290 K, the sums of the calculated losses for the rotor and the specimen at different speeds are listed in Table 2.

**Table 2.** Air friction losses in the test rig with *p* = 0.05 Pa and *T* = 290 K.


## *3.5. Simulation of Thermal Rotor Behavior*

The thermal simulation was performed using a 3D-FE model of a quarter of the test rig, both with and without the specimen. The rotor transmits thermal energy to the stator via radiation, and vice versa. It was assumed that the emissivity and absorptivity of the rotor and stator surfaces do not depend on the wavelength or direction of the radiation. Under this assumption, the emissivity and absorptivity of the surface are equal [23]. For real materials, only a part of the radiation that hits a surface is absorbed while the rest is reflected, hence, the coefficient of emission is smaller than one. When considering the heat exchange between two surfaces through radiation, both surfaces have to be taken into account. Consequently, for the simulation of the rotor heating, the inner surface of the stator has to be included in the model. The stator and rotor were painted with lacquer on the surface of the AMBs and the PMSM. The coefficient of emission was assumed to be 0.9 for all painted parts and the FRP. For the unpainted, blank parts of the steel rotor, a coefficient of emission of 0.3 was assumed, and for the aluminum parts of the stator and the specimen a value of 0.05 was assumed [23]. The temporal stator temperature profile of the active components was roughly approximated with three measured temperatures, which were linearly interpolated. The temperature values can be found in Appendix A. This transient temperature change has only a small impact on the simulation results; using a constant mean temperature value for all stator components increases the calculated rotor temperature only about 1 ◦C. For the rest of the stator, a constant temperature of 25 ◦C was predefined. To obtain the temperature distribution in the rotor, heat conduction in the rotor was also taken into account in the model, with a conductivity of 25 W/Km. Calculating the steady state rotor temperature during the test cycles requires the mean loss value per cycle for each actuator. As a baseline, a cycling time of 30 s was defined—i.e., during one cycle, the rotor constantly accelerates for 15 s from 15,000 rpm to 30,000 rpm, and then constantly decelerates with the same slope back to 15,000 rpm. The rotational speed and the speed dependent losses of the active components are illustrated in Figure 8a. The mean values of all calculated losses are shown in Figure 8b.

**Figure 8.** (**a**) Rotational speed and calculated rotational losses of the active components during one cycle; (**b**) Resulting mean losses of all calculated loss mechanisms in the rotor during cycling.

## *3.6. Temperature Measurements*

To validate the described models, the rotor temperature was measured during cycling without a specimen. Two infrared sensors measured the rotor temperature: one between the upper radial AMB and the thrust disk and one next to the lower radial position sensors (see Figure 2b). The utilized sensors can measure the temperature of metallic surfaces above 50 ◦C, with an uncertainty of ±2 ◦C. First, the cooling of the levitated but not spinning rotor from a known temperature was analyzed. This excluded the rotational losses and the switching losses in the PMSM. Hence, only the small switching losses of the AMBs were present. The results of a transient simulation and the measurement are shown in Figure 9a. With a maximum deviation of 1.85 ◦C, the simulation results lie in the range of the measurement uncertainty. Consequently, the switching losses of the AMBs and the modeling of the heat transfer are sufficiently accurate. Next, the thermal behavior of the rotor during rotation was analyzed. Before the measurement, the system was turned off for 12 h, to ensure the rotor was cooled down to room temperature. During the measurement, the PMSM continuously accelerated and decelerated the rotor without a specimen from 15,000 rpm to 30,000 rpm. In the transient simulation, the mean losses from Figure 8b were used. The initial speed up to 15,000 rpm, which takes about 15 s, was neglected. Figure 9b compares the measured temperatures at both sensor positions with the simulation. While the rotor temperature at both positions was nearly the same in the simulation, it differed considerably in the measurement. The model underestimated the rotor heating at both measurement positions. At the lower position, the deviation was 5.2 ◦C, while it was 22.4 ◦C at the upper position. Due to the good agreement of the rotor cooling, the modeling error is expected to lie in the loss calculation rather than the thermal model.

**Figure 9.** (**a**) Simulated and measured rotor cooling without the specimen; (**b**) Simulated and measured rotor temperatures without the specimen during cycling.

One reason why the losses could be underestimated can be found in the loss coefficients of the NO20 sheets in the AMBs and PMSM. The loss measurements done by the manufacturer that were used to identify the loss coefficients are normally performed on unmachined sheets after a precise annealing process. However, the material behavior can change during milling and heating. For example, a damaged sheet isolation can increase eddy current losses by up to 30% [24] and hysteresis losses can vary more than 50% depending on their heat treatment [25]. Because these changes can hardly be quantified, the calculated losses are adjusted to fit the temperature measurement. At the lower AMB, the calculated losses were increased by 10% to fit the temperature at the lower temperature sensor. Since the sheets at the upper and lower AMB were milled and heated similarly, the losses in the upper AMB increased by 10% as well. To further match the simulation with the measured temperature at the upper sensor, 8 W of additional losses in the axial AMB were required in the simulation. For the loss calculation, it was assumed that the axial AMB is perfectly symmetrical with respect to the rotational axis. In reality, small circumferential variations in the material of stator and rotor lead to flux changes during rotation and hence to higher rotor losses. Furthermore, control activities in the AMBs were neglected in the simulation, but during the measurement, considerable movement in the axial direction was observed. The adjusted losses are summarized in Table 3. A transient thermal simulation with the adjusted loss values showed good accordance with the measured data (see Figure 10a).

**Table 3.** Adjusted mean losses in the rotor during cycling.

**Figure 10.** (**a**) Simulated and measured rotor temperatures with adjusted AMB losses without the specimen during cycling; (**b**) Simulated steady state temperature of the rotor with and without the specimen.

With the adjusted loss values, a steady state thermal simulation of the test rig with and without the specimen was performed. Without the specimen, the motor needs an rms phase current of 3.8 A to accelerate the rotor during the described cycle, but with the specimen an rms phase current of about 37 A is required. The losses of the PMSM were considered accordingly in the simulation. The other losses in the rotor were assumed to be equal for both cases. The simulated temperature distribution in the center of the rotor over the axial length is shown in Figure 10b. In the specimen, the temperature between the FRP and the steel segments is displayed. The maximum temperature without a specimen was 176 ◦C, located at the upper AMB. Even though higher rotor losses occur with the specimen, the rotor temperature was 10 ◦C lower than without it, as the large surface area of the FRP benefits heat transfer to the surrounding containment. The thin tapering at the hub led to a large temperature drop between rotor and specimen, preventing the temperature of the FRP from exceeding 40 ◦C. The operation temperature of the FRP in the flywheel must not exceed 80 ◦C. The magnets of the PMSM should not be operated over 120 ◦C. In order to have a safety margin, the magnet temperature should not exceed 110 ◦C. The simulation results, however, show a temperature of the magnets of over 140 ◦C. Thus, to avoid reducing the cycle time, rotor losses have to be reduced.

## **4. Loss Reduction**

As seen in the previous section, the highest power losses in the rotor were due to rotational losses in the radial AMBs and PMSM, as well as switching losses in the PMSM, which all a have similar order of magnitude. Losses due to control activities in the axial AMB were slightly lower. Switching in the AMBs and air friction can be neglected. Hence, the main focus should be on the reduction of rotational losses in the radial AMBs and the PMSM as well as switching losses in the PMSM. In this paper only the radial AMBs will be further discussed.

## *4.1. Reducing Losses in the Radial Bearings*

One effective way to reduce the rotational losses of an AMB is to reduce its bias current *IB*. For example, reducing *IB* from 5.67 A to 4 A reduces the rotational losses in the rotor by more than 50%, to 8.2 W in the upper AMB and 6.8 W in the lower AMB. A thermal simulation and measurement with *IB* = 4 A was performed, where the rotor was again cycled according to Figure 8a. The results are shown in Figure 11. The transient simulation and the measurement of the test rig without a specimen again show sufficient correspondence (see Figure 11a). The steady state simulation predicted a maximum temperature of 129 ◦C at the axial AMB and 123 ◦C at the PMSM without a specimen, as well as 127 ◦C and 117 ◦C with a specimen (see Figure 11b). Changing *IB* also changes the rotor dynamics in the AMBs. To evaluate the possibility of a bias reduction, the dynamics of the levitated rotor will be investigated in the next section.

**Figure 11.** (**a**) Simulated and measured rotor temperatures with a reduced bias current of 4 A in both AMBs without the specimen during cycling; (**b**) Steady state rotor temperature with a reduced bias current of 4 A in both AMBs with and without the specimen.

## *4.2. Rotor Dynamics*

A mechanical 3D-FE model of the rotor was created and evaluated. For a modal analysis of the model, the AMBs were modeled as springs with a constant stiffness of 2.5 <sup>×</sup> 105 N/m. The first bending mode of the rotor without a specimen lay at 1,030 Hz, thus 353 Hz above the highest rotational frequency. The calculated second bending mode was at 2,170 Hz. Measured Campbell diagrams of the rotor without the specimen are shown in Figure 12. They were derived from the signals of the radial position sensors of the AMBs during an acceleration of the rotor from 0 rpm to 40,000 rpm in 110 s. The signals were filtered with a 5 kHz low-pass filter and sampled with 10 kHz. To create the Campbell diagram, the measurement was divided into 4096 sample long windows, and for each window a Fourier transformation was performed. Figure 12a shows the Campbell diagram of the position data

of the upper AMB with *IB* = 5.67 A. The lines starting at the lower left corner are harmonics of the rotational speed, the thickest one being the first harmonic. Odd harmonics of the rotational speed are clearly visible, whereas even harmonics are barely detectable. The approximately horizontal lines are eigenfrequencies (EFs) of the system. At about 90 Hz lies the rigid body tilting mode of the rotor in the AMBs. The translational rigid body mode is not visible in the diagram. The first bending mode of the rotor can be seen at 1057 Hz, the second at 2063 Hz. The second splits into a forward and backward mode. While the calculated first bending mode only differed by 27 Hz, the second bending mode was about 100 Hz lower than the simulation predicted. Underneath each Campbell diagram, the mean deflection of the rotor for each window is shown. The deflections were less than 20% of the backup bearing clearance of 200 μm. Peaks can be seen when the rotational speed or its harmonics hit an EF. The first two peaks were due to the rigid body modes. The next peak at 25,000 rpm was caused by the fifth harmonic reaching the forward mode of the second bending EF, and the last peak at 39,000 rpm was caused by the third harmonic reaching the backward mode of the second bending EF. Smaller peaks can also be seen when higher harmonics reach the second bending EF, where every second odd harmonic excited its forward mode and every other odd harmonic excited the backward mode. A notch filter was placed at 1 kHz in the AMB control, so there was no additional excitation of this EF due to harmonics of the rotational speed. Figure 12b shows the same evaluations for *IB* = 4 A. The harmonics were not influenced by *IB* and the EFs only changed slightly. However, the amplitude of the deflection differed noticeable. The peak where the first harmonic crossed the second rigid body mode was about 30% higher with a reduced *IB*, whereas the peaks at higher speeds were now below 20 μm. Because the normal speed range starts at 15,000 rpm, the rigid body mode has to be crossed only once per test, while the peak at 25,000 rpm is reached twice every cycle. Consequently, without a specimen, reducing *IB* does not only reduce the rotational losses, but also reduces the control activity of the AMB during cycling.

**Figure 12.** (**a**) Measured Campbell diagram and deflection of the rotor with a bias current of 5.67 A without specimen; (**b**) Measured Campbell diagram and deflection of the rotor with a reduced bias current of 4 A without specimen.

Measurements with a specimen were not performed yet, but Figure 13a shows the calculated EFs and Figure 13b the simulated Campbell diagram of the rotor with the specimen. The first elastic mode was at 9 Hz, where only the specimen tilts at the tapering. At speeds higher than 9 Hz, the CoG of the specimen did not move translationally. In the second, third, and fourth modes, the specimen was tilting around its CoG and in the fifth mode it does not move at all. The rotor showed elastic deformations beginning with the fourth mode. Due to the high inertia ratio of the specimen, the higher EFs showed a strong dependency on the rotational speed, which can be seen in the Campbell diagram in Figure 13b. The first three EFs were below 250 Hz, where the fatigue testing begins, and have to be passed while starting a test. The fourth EF increased in such a way that it is not reached by speed-synchronous excitations within the operation range. However, due to the close proximity of the fourth mode, the bandwidth of the AMBs has to include this frequency. The fifth EF had a minimum distance of 290 Hz to the rotational speed and therefore it will have to be filtered out from the position signal.

**Figure 13.** (**a**) Calculated bending eigenmodes below 1 kHz of the rotor with sample; (**b**) Calculated Campbell diagram of the rotor with the specimen. Each eigenmode shown in (**a**) splits in a forward and a backward mode.

For the differential winding design, which is used for the radial AMBs, the control current has to be smaller than *IB*. To allow control currents bigger than *IB* the AMBs would have to be rebuilt for a differential control design and a nonlinear control might be needed (see [15,26–28]). Both would require major changes in the hardware. Hence, with the existing hardware, reducing *IB* also reduces the maximum force the actuators can generate. In the centered position, one AMB can generate a force of 278 N with *IB* = 5.67 A, but only 144 N with *IB* = 4 A. For the rotor without a specimen this is unproblematic, because there are nearly no external loads in the test rig and the rotor can be precisely balanced to minimize unbalance forces. However, the initial unbalance of an assembled specimen is big and balancing it precisely is difficult. The force generated by the AMB has to be big enough to counteract the unbalance force while passing the first three EFs. This has to be further investigated.

## *4.3. Increasing the Cycle Time*

When no further loss reduction can be achieved, the cycle time has to be changed. In the upper example, even with a bias current of 4 A, the rotor temperature were above 110 ◦C but below 120 ◦C, as shown in Section 4.1. If further bias reduction or other loss reduction strategies are not applicable, a less favorable possibility to reduce the rotor temperature is to increase the cycle time. Here, two variations were considered, which are illustrated in Figure 14a. In the first variation, the acceleration of the rotor is reduced to lower load dependent losses in the PMSM, in the second variation, after a cycle, the rotational speed is held constant at 15,000 rpm for a while to allow the rotor to cool down. How long the cycle has to be increased depends on the associated loss reduction and the maximum allowable rotor temperature. For example, the cycle time could be increased from 30 s to 37 s to reduce the rms motor phase current from 37 A to 30 A, and thereby the rotor losses in the PMSM by 1.7 W to 17.4 W. Losses in the AMBs stay the same. Alternatively, after a 30 s cycle, the speed is held constant at 15,000 rpm for 7 s. Here, the mean losses can be calculated by averaging the mean losses for the first 30 s, as shown in Sections 3.5 and 3.6, and the losses at a constant speed. To keep the rotor without a specimen at a constant speed of 15,000 rpm, an rms phase current of 1.8 A is needed in the PMSM. Hence, it was assumed that an rms phase current of less than 10 A is needed to keep the rotor with a specimen at constant speed. Below 10 A, losses were nearly independent of the phase current (see Section 3.3). The resulting mean losses are summarized in Table 4. For the axial AMB constant, rotor losses of 8.1 W were assumed for both variations. The total cycle time in both variations was 37 s and the whole test would take the same amount of time—86 days to perform 200,000 cycles.

**Figure 14.** (**a**) Variations of increasing the cycle time to reduce the rotor temperature; (**b**) Steady state rotor temperature of both cycle variations with a reduced bias current of 4 A in both AMBs with a specimen.


**Table 4.** Mean rotor losses during cycling variation 1 and 2 with a reduced bias current of 4 A.

Since losses in the PMSM showed only a small dependency on the load, the total mean losses in variation 2 were smaller than in variation 1. The calculated rotor temperature is shown in Figure 14b. For variation 1, the maximum calculated temperature in the permanent magnets of the PMSM was 115.6 ◦C, and for variation 2 it was 111.0 ◦C. Thus, with a slightly increased cooling time at 15.000 rpm, variation 2 should not reach the critical magnet temperature, whereas variation 1 has no big benefit for the rotor temperature. However, in the model the thermal conduction between rotor and magnets of the PMSM is assumed as ideal. In reality, small gaps could decrease the conductivity and therefore increase the dependency of the magnet temperature on the internal losses of the PMSM. This would make the loss reduction in the PMSM critical compared to the other components, giving variation 1 higher relevance.

## **5. Discussion**

This paper proposed a structured approach to evaluate the heating of a magnetically levitated rotor in a vacuum, and deduced measures to avoid critical temperatures. First losses and the subsequent rotor heating were calculated. The rotational losses in the radial AMBs and the PMSM were identified as the main rotor loss mechanisms in the test rig. Losses in the axial AMB during rotation are smaller but also in the same order of magnitude. Switching losses in the AMBs and air friction losses only have a very small contribution to the total losses. One way to achieve loss reduction in the radial AMBs, without major changes in the spindle hardware, is to reduce the bias current. While doing so, the controllability of the rotor has to be ensured. The applicability of this approach was shown for the rotor without a specimen, but was not yet shown with a specimen. Stability problems could arise at low rotational speeds where many EFs can be excited. An adjustment of the bias with the rotational speed might be necessary. Further loss reduction can be accomplished by reducing the control activity. A centralized control of the radial AMBs should reduce the activities in the upperradial AMB. To realize a centralized control, however, the amplifier needs to be replaced. An overview of loss reduction strategies in AMBs can be found, for example, in [15]. Losses in the PMSM can be reduced by increasing the switching frequency of its inverter or by employing an output filter. Both options also require changes in the hardware. If no further loss reduction can be achieved, the cycle time has to be adjusted, so that the critical rotor temperature is exactly reached. In the present case, the cycle time had to be increased. With a cycle time of 37 s and a bias current of 4 A in both radial AMBs, a calculated steady state temperature of 111 ◦C can be reached, however, this has the drawback of an increased test duration. To account for modeling errors, the cycle time can further be adjusted during operation, depending on the rotor temperature, in a closed control loop. For the calculations, a constant thermal conductivity in the rotor was assumed. This assumption is especially critical between the magnets and the rotor. In reality, the thermal coupling will be poorer, resulting in a higher magnet temperature.

**Author Contributions:** Conceptualization, D.F., M.S. and M.R.; methodology, D.F., M.S.; validation, D.F.; formal analysis, D.F.; investigation, D.F.; resources, D.F.; data curation, D.F.; writing—original draft preparation, D.F.; writing—review and editing, D.F., M.S., M.R. and S.R.; visualization, D.F..; supervision, S.R.; project administration, D.F. and S.R.

**Funding:** This research was funded by the German Federal Ministry for Economic Affairs and Energy, grant number 03ET6064A.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

## **Appendix A**

The following Table A1 shows the measured stator temperature of the test rig during cool down. Table A2 shows the same temperatures during cycling with a bias current of 5.67 A and Table A3 shows them during cycling with a bias current of 4 A. These temperatures were used in their corresponding simulations.


**Table A1.** Measured stator temperature during cool down.

**Table A2.** Measured stator temperature during cycling with a bias current of 5.67 A.



**Table A3.** Measured stator temperature during cycling with a bias current of 4 A.

## **References**


© 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/).

## *Article* **Design and Analysis of a 1D Actively Stabilized System with Viscoelastic Damping Support †**

**Josef Passenbrunner 1,\*, Gerald Jungmayr <sup>1</sup> and Wolfgang Amrhein <sup>2</sup>**


Received: 5 March 2019; Accepted: 15 April 2019; Published: 17 April 2019

**Abstract:** Passively magnetically stabilized degrees of freedom yield the benefit of reduced complexity and therefore costs. However, the application of passive magnetic bearings (PMBs) also features some drawbacks. The poor damping capability leads to exaggerated deflection amplitudes when passing the resonance speeds of the applied system. This results in the necessity of external damping. Complying with the goal of costs and complexity, viscoelastic materials offer a suitable solution. However, these materials show high frequency and temperature dependent properties which induce the necessity of a proper model. Thus, the design of systems, as presented in this paper, requires accurate modeling of the dynamic behavior including the nonlinear characteristic of damping elements to predict the system displacements. In the investigated setup only two degrees of freedom remain to be controlled actively. These are the axial rotation and the axial position of the rotor which are controlled by the motor and an active magnetic axial bearing (AMB). This article focuses on the rotor dynamic modeling of a radial passively magnetically stabilized system especially considering the nonlinear behavior of viscoelastic damping elements. Finally, the results from the analytic model are verified by measurements on a manufactures test system.

**Keywords:** passive magnetic bearing; viscoelastic material; damping; modeling; rotor dynamics; active magnetic bearing; cost reduction

## **1. Introduction**

Nowadays, the number of applications for rotating machines is huge. In general, therefore, the rotor is supported with traditional mechanical ball or slide bearings, which feature the drawback of applied lubricants, mechanical friction and resulting wear. Magnetic bearing technology is a contemporary research area that can be described as relatively young. The principle is based on the generation of suspension forces by magnetic fields which allows a contact-free operation of the motor. Magnetic bearings have their advantages especially in areas that require high-purity operation or very high speeds. These applications can be found in machining technology as milling spindles, in vacuum technology as turbomolecular pumps, in medical applications as blood pumps as well as in flywheel accumulators [1,2] and gas compressors. Additionally, the aspect of service life is superior, due to the fact that in magnetically suspended systems it is only determined by the electronic components.

In addition, as is usual in the field of magnetically stabilized system design, the costs play a very important role. The high complexity compared to standard ball bearings has thereby a major impact on the price. With reference to fully actively stabilized systems, which are still applicable in commercial

areas, investigations regarding passively stabilized systems are rare [3,4]. The first works on passive permanent magnetic (PM) ring bearings can be found in 1976 [5]. In 1981, Yonnet [6] characterized the different possible bearing principles for the first time. Marinescu and Yonnet presented the basics for dimensioning of PM ring bearings in 1979 and 1980 [7,8]. A great advance in analytical computability was presented by Lang [9] for attractive ring bearings. An extension to repulsive ring bearings was provided by Jungmayr [10]. Today the filed of applications for passive PM ring bearings include spinning centrifuges, turbomolecular pumps, flywheels and fully magnetically stabilized fans.

The use of passive magnetic bearings (PMB) offers a very interesting approach to reduce the complexity of the magnetic bearing to a minimum and thereby makes the application attractive for cost-sensitive areas. Hence, in such systems the demand for power electronics and position sensors is minimized. The passive bearing arrangement by means of permanent magnets has to be emphasized, as it enables a nearly loss-free bearing arrangement with very little effort. Furthermore it features the advantage of easy determination of stiffness and static load carrying capacity. The determination of the stiffness of permanent magnet bearings can be performed analytically in arrangements without ferromagnetic material. For applications with ferromagnetic material finite element calculations are necessary. However, a major disadvantage of permanent magnetic bearings is the extraordinary low damping of the stabilized degrees of freedom, which lead to exaggerating deflections when passing rigid body resonances during a run-up process. So additional damping has to be induced into the system [11–14]. Especially systems with large unbalances need a very precise consideration of the dynamic behavior.

A possibility to provide damping is offered by viscoelastic materials which meet the targets of low costs and simplicity. Viscoelastic materials are currently being used in various applications for damping and suppression of vibrations. The selection of the damping material itself is often empirical. One reason for this lies in the complex dynamic behavior of this type of materials. In most cases the stiffness and damping of elastomers feature a high dependency on the excitation frequency and temperature. Using viscoelastic materials it is no longer possible to actively influence the system dynamics subsequently. Thus, a precise consideration of the material characteristics is necessary. This implies an exact knowledge of the overall system, especially the rotordynamic behavior, to ensure proper operation. Occurring deflections, due to the rotor dynamics, have to meet certain restrictions, to avoid contact between the rotor and the stator. This requires the derivation of the systems equations of motion and a proper model of the frequency-dependent damping elements behavior.

The concept examined in this paper is based on an industrial application, which is located in the area of mass production. With a large number of manufactured systems the price plays a decisive role. Thus, only system concepts come into consideration, which can drastically reduce costs. The avoidance of actively stabilized degrees of freedom and the resulting hardware savings (electronics, sensor technology, bearing coils) are crucial. Concepts with passively magnetically stabilized degrees of freedom are favorable due to their simplicity and low costs. However, a solution can only be achieved if the above-mentioned problems of damping can be solved in a cost-effective way. Although the modelling of viscoelastic damping elements is very complex, it offers an approaches to set up a proper system description.

The aim of this work is the design and optimization of a low-cost magnetic bearing drive system taking into account the special boundary conditions of the application. The high reliability of the concept is to be demonstrated, helping to implement contactless bearings in the industrial field more often.

In Section 2 the system setup and components are described. Section 3 determines the viscoelastic material model followed by Section 4, where the derivation of the equations of motion is presented. With the derived model an optimization is performed in Section 5. Finally, the optimized system is verified by measurements in Section 6.

## **2. System Setup**

Based on the application considered, the rotor features a vertical alignment. The speed range is 10.000–30.000 rpm, which, however, is constantly interrupted. Hence, a frequent start and stop of the system is demanded. A very critical point are unbalances, which will lead to high deflections in the passive bearings. These unbalances occur mainly due to constantly changing additional process masses, which are attached to the rotor during operation. Thereby the unbalance maximum is defined by a value of 20 gmm and can be applied to the rotor at any position. The drive system itself can be balanced in order to reduce the unbalance to a minimum. However, the unbalance of the process mass must be taken as the occurring unbalance. Especially in connection with a PMB concepts, this represents an immense challenge and requires an exact consideration of the damping elements.

Geometric restrictions exist on the outer diameter of the rotor which is given by a maximum value of 42 mm. This requires a very compact design of the system. The limited space conditions lead to a stacked vertical construction. Figure 1 shows the setup of the considered magnetically levitated system.

**Figure 1.** Setup of the investigated system.

In the shown configuration, the radial deflections and the tilting are stabilized by two radial permanent magnet bearings. These ringbearings are of repulsive type and use axial magnetized magnets, which are easy to produce and lower the costs for the bearings. Furthermore, this configuration offers the opportunity of magnet stacking [15] to achieve higher bearing stiffness with the same cross section and thereby simplifies the rotordynamic design procedure.

Without ferromagnetic material near the PM the stabilizing radial stiffness of the bearings can be analytically calculated [9]. This drastically simplifies the bearing design for the prototype, because no finite element simulations are necessary. Due to Earnshaw's theorem [16], which is adapted for passive ring bearings in [17], the destabilizing axial stiffness *sz* of the bearings is given by

$$\mathbf{s}\_z = -\mathbf{2} \cdot \mathbf{s}\_{\mathbf{r}} \tag{1}$$

where *sr* describes the stabilizing radial stiffness.

Hence, at least one direction needs to be stabilized actively. In between the active axial bearing, the motor and the position sensor [18,19] are placed. As damping elements viscoelastic ring elements are used. These elements are located between the stator and the system housing. The viscoelastic ring elements are placed and glued in two aluminum rings to allow easy mounting. With this concept and the right dimensioning, it is possible to achieve improved dynamic system properties [20]. The rotor is designed as an exterior rotor. In principle, interior rotor concepts are also possible, but it should be avoided that the critical bending Eigenfrequencies of the rotor occur within the speed range. Moreover, the supporting rod of the motor, the active bearing and the sensor system become thin. Hence, the flexible behavior might as well affects the dynamics of the system. However, the considered motor is constructed with a slotted stator and ferrite magnets for cost reduction. As AMB a reluctance force bearing is used, whose force density is higher compared to a Lorenz force bearing.

## **3. Model of the Viscoelastic Behavior**

Viscoelastic materials feature a frequency dependent behavior of the stiffness and damping values. To describe the characteristics of such materials often a generalized Maxwell model is used [21–23]. As shown in Figure 2, such a model consists of a single spring with the equilibrium modulus *E*<sup>0</sup> and several Maxwell units in parallel. This spring *E*<sup>0</sup> describes the material response after infinite time. Each Maxwell consists of a single spring and a single damper in series, reproducing the frequency dependency by adding different time constants *τn*. So a quasi non-linear behavior can be represented by superposing linear elements. For harmonic excitations in the frequency domain a representation of the generalized Maxwell model with a complex modulus

$$E(\omega) = E'(\omega) + jE''(\omega) \tag{2}$$

is useful. Thereby, *E* stands for the storage module and *E* for the loss module. The quotient gives the loss factor

$$
\eta = \tan \delta = \frac{E''}{E'}.\tag{3}
$$

Converted to the Prony parameters [24] of the generalized Maxwell model the components of the complex modulus result in

$$E'(\omega) = E\_0 + \sum\_{n=1}^{N} E\_n \frac{\omega^2 \tau\_n^2}{1 + \omega^2 \tau\_n^2} \tag{4}$$

and

$$E^{\prime\prime}(\omega) = \sum\_{n=1}^{N} E\_n \frac{\omega \tau\_n}{1 + \omega^2 \tau\_n^2} \tag{5}$$

with the time constants *τ<sup>n</sup>* = *dn En* .

In addition, the Maxwell model uses *N* inner states *yn*. A main advantage of this model is its easy integration as it can be described by a system of linear differential equations.

**Figure 2.** MAXWELL-Modell of the viscoelastic materials dynamic behavior representing one viscoelastic ring support.

## *Determination of Material Parameters*

To describe the thermo-viscoelastic behavior usually master curves are used [24]. Thereby, the theory of temperature-time-correspondence is applied to combine measurements at different temperatures and draw conclusions for other frequencies. The basic approach is depicted in Figure 3.

**Figure 3.** Application of the temperature-time-correspondence.

In theory [25] viscoelastic material shows the same storage modulus for a certain temperature *T*<sup>1</sup> and frequency *f*<sup>1</sup> as for a different temperature *T*<sup>2</sup> and a referring scaled frequency *f*<sup>2</sup> according

$$E(f\_1, T\_1) = E(f\_1 \cdot a\_{T\_{2-1}, \prime} T\_2) = E(f\_2, T\_2) \tag{6}$$

with the shift factor *aT*<sup>2</sup>−<sup>1</sup> . With lower temperatures the storage modulus is increasing and falling with increasing temperatures. As only a limited frequency range can be measured, one measurement is not sufficient to represent the required frequency range. Hence, it is necessary to combine the measurements. The shift factors *aT* thereby perform a vertical shift of the measured data, so that a smooth overall mastercurve is resulting. If unfilled elastomers are used, the kinetic-theory-factor [25] has to be applied, which is obtained from the theory of entropy elasticity. The storage modulus is thereby horizontally shifted to the reference temperature *Tr* using

$$E(T\_r) = \begin{cases} E(T\_m) \* \underbrace{\frac{T\_r}{T\_m}}\_{k\_T} \* \underbrace{\frac{\rho\_r}{\rho}}\_{\approx 1} & T\_m \ge T\_{\%} \\ E(T\_m) & T\_m < T\_{\%} \end{cases} \tag{7}$$

with *Tm* as measured temperature and *Tg* as the glass transition temperature. The factor *ρr*/*ρ* normalizes the specific volume at a temperature *Tm* to the reference temperature *Tr*. Therefore, the horizontal shift is only applied to temperatures lower than the glass transition temperature, which lies for technical elastomers normally far below zero degrees Celsius. The loss factor values are also shifted in that process as it is the quotient of *E* and *E*.

Reliable master curve data is hardly provided by most manufacturers and even if available it is essential to know the exact measurement conditions. For example a measurement under precompression shows a very different behavior compared to the same measurement without precompression. Hence, it was decided to measure the data on our own. Even though the measurement principle looks simple, it is a very tricky task. In Figure 4 the employed measuring method is shown. It is based on the so called Dynamic-Mechanic-Temperature-Analysis (DMTA).

**Figure 4.** Scheme of the employed measuring method for the viscoelastic materials.

A specimen is excited at different temperatures by harmonic deformation *x*(*t*) at various frequencies. In the process the reaction force *F*(*t*) is measured by a load cell. Thereby, the ratio between force and displacement reflects the stiffness of the material. Due to dissipation effects a phase shift *δ* occurs between the force and the displacement, which represents the damping capability of the material. The measurement leads to the so-called isotherms. Thereafter, the measured stiffness values [26] have to be converted to the material significant storage modulus *E* and loss modulus *E* using the geometric parameters of the specimen. For a circular specimen with the height *hs* and diameter *ds*the relation between the measured axial stiffness *kax* and the material modulus *E* is given by

$$k\_{ax} = E \cdot \frac{A}{h\_{\mathfrak{s}}} \cdot (1 + 2S^2) \quad \quad S = \frac{d\_{\mathfrak{s}}}{4h\_{\mathfrak{s}}} \quad A = \frac{d\_{\mathfrak{s}}^2 \* \pi}{4}.\tag{8}$$

In Figure 5a,b the measured isotherms of a selected butyl rubber with a hardness of Shore A40 are shown. It can be seen that the storage modulus is increasing with lower temperature. However, the loss factor shows, at the beginning, an increase with lower temperature till it reaches a maximum value. Afterwards, the loss factor falls again.

**Figure 5.** Storage modulus (**a**) and loss factor (**b**) of the measured isotherms.

For this material also measurement data from the manufacturer is available. It was observed that the measured data shows a 50% higher stiffness and a 30% lower loss factor compared to the given data provided by the manufacturer. That proves the difficulty to get proper material specifications. The isotherms are afterwards shifted using Equation (6) processing the vertical shift and (7) for the horizontal shift of the data below the glass transition temperature to obtain the master curve considering a smooth gradient of the stiffness and loss factor. However, there are other possibilities to shift the data. Basically the procedure to find the optimal shift parameters and consequently the identification of the Prony-parameters requires the solution of a nonlinear minimization problem. In this work a genetic algorithm is used varying the shift parameters and determineing the Prony-parameters by minimizing the mean square deviation. The results are plotted in Figure 6. The fitted model shows a very good compliance with the measured data.

**Figure 6.** Fitted curves of the storage modulus (red) and the loss factor (blue) in comparison to the measured data (×) at a reference temperature of 25 ◦C, whereby the different colors mark the various measurement temperatures.

## **4. Rotordynamic Model**

In the considered system, high rotor unbalances are induced by variable process masses acting on the rotor. Obviously, the relative deflection between stator and rotor in the PMB planes are important. To determine the movements of the system bodies the equations of motion have to be derived. Due to the system setup, the AMB, the motor and the upper PMB, are located on a thin shaft. Therefore,

also the bending of this part have to be considered. The equations of motion are conducted by using the projection equation [27], which projects the (generalized) forces into the unconstrained space, where the motion takes place. These forces are given by

$$\sum\_{i=1}^{N} \left[ \left( \frac{\partial\_R \boldsymbol{\sigma}\_{si}}{\partial \dot{\boldsymbol{q}}} \right)^{T} \left( \frac{\partial\_R \boldsymbol{\omega}\_{si}}{\partial \dot{\boldsymbol{q}}} \right) \right] \begin{bmatrix} (\_{R}\dot{\boldsymbol{p}} + \_{R}\boldsymbol{\tilde{\omega}}\_{IR\ R} \boldsymbol{p} - \_{R}\boldsymbol{f}^{s})\_{i} \\ (\_{R}\boldsymbol{\mathcal{L}} + \_{R}\boldsymbol{\tilde{\omega}}\_{IR\ R} \boldsymbol{\mathcal{L}} - \_{R}\boldsymbol{\mathcal{M}}^{s})\_{i} \end{bmatrix} + \left( \frac{\partial V}{\partial \boldsymbol{q}} \right)^{T} + \left( \frac{\partial R}{\partial \dot{\boldsymbol{q}}} \right)^{T} = 0,\tag{9}$$

with *<sup>R</sup>vsi* and *<sup>R</sup>ωsi* describing the bodies velocities and angular velocities, while *<sup>R</sup> p* and *<sup>R</sup>L* represent the impulse and angular momentum in the reference system *R*. Thereby,

$$\left(\frac{\partial\_R \boldsymbol{\sigma}\_{si}}{\partial \dot{\boldsymbol{q}}}\right)^T \qquad \text{and} \qquad \left(\frac{\partial\_R \boldsymbol{\omega}\_{si}}{\partial \dot{\boldsymbol{q}}}\right)^T \tag{10}$$

represent the Jacobian matrices and the therms

$$\left(\frac{\partial V}{\partial q}\right)^T \quad \text{and} \quad \left(\frac{\partial R}{\partial \dot{q}}\right)^T \tag{11}$$

are used to consider potential forces (e.g., springs, dampers, elastic potential and gravitation). *<sup>R</sup> f <sup>s</sup>* and *<sup>R</sup>M<sup>s</sup>* describe forces and torque acting on the center of mass.

The rotor is modeled as rigid body, whereas the stator is split into a rigid lower part and a flexible upper part, because the moment of resistance of the lower part is much higher due to the increasing diameter. The mass of the AMB and the motor are integrated into the rigid part of the stator to reduce the complexity of the model. The upper PMB is modeled as point mass. As flexible beam model a Ritz approach based on a cubic function *u*(*z*) = *a*<sup>0</sup> + *a*1*z* + *a*2*z*<sup>2</sup> + *a*3*z*<sup>3</sup> is used to describe the occuring displacements of the flexible stator part. Figure 7 shows the comparison of the exact and approximated first and second Eigenmodes of a semibeam. As can be seen, the cubic approach gives a very good representation of the first bending mode which is the most important for the investigated system.

**Figure 7.** Comparison of the exact and approximated Eigenmodes of a semibeam bending.

For modelling a linear spring elements the potential energy is used. In the system such linear springs can be used as representation for


The potential energy of a linear spring is determined by

$$\mathcal{V}\_{\mathcal{F}} = \frac{1}{2} \Delta \mathfrak{x}(q) \ K \,\Delta \mathfrak{x}(q), \tag{12}$$

where Δ*x*(*q*) describes the relative extension vector of the spring from the force-free position as a function of the generalized coordinates *q* and *K* represents the tensor of spring constants. The same approach can also be used to describe the influence of elastic elements to their attachments. The speed-proportional potential energy of a damper is given by

$$\mathcal{R}\_D = \frac{1}{2} \Delta \dot{\mathbf{x}}(\mathbf{q}, \dot{\mathbf{q}}) \text{ D } \Delta \dot{\mathbf{x}}(\mathbf{q}, \dot{\mathbf{q}}). \tag{13}$$

*D* describes the tensor of damping constants and Δ*x*˙(*q*, *q*˙) the relative velocity vector of the damper dependent on *q* and ˙*q*.

The equation of motion obtained from the projection equation generally shows the structure

$$M(q)\ddot{q} + \mathbf{g}(q, \dot{q}) - \mathbf{Q}(q, \dot{q}) = 0. \tag{14}$$

and is of non-linear character. *M*(*q*) is the skew symmetric mass matrix. *g*(*q*, *q*˙) contains the remaining terms of the acceleration like centrifugal- and coriolis forces. In *Q*(*q*, *q*˙) the applied spring-, damperand weightforces as well as the actuating forces and moments take place. If only small deflections from the rest position *q* = *q*<sup>0</sup> + Δ*q* are assumed, Equation (14) can be developed by a TAYLOR series. A development up to terms of first order is sufficient and only the linear differential equation for small deflections remains

$$\mathbf{M}\_0 \Delta \ddot{q} + \mathbf{P}\_0 \Delta \dot{q} + \mathbf{Q}\_0 \Delta q = 0 \tag{15}$$

with

$$\mathcal{M}\_0 = \mathcal{M}(q\_0), \quad \mathcal{P}\_0 = \left. \frac{\partial(\mathbf{g} - \mathbf{Q})}{\partial \dot{\mathbf{q}}} \right|\_{q\_0}, \quad \mathcal{Q}\_0 = \left. \frac{\partial \mathcal{M} \ddot{\mathbf{q}}}{\partial q} \right|\_{q\_0} + \left. \frac{\partial(\mathbf{g} - \mathbf{Q})}{\partial q} \right|\_{q\_0}. \tag{16}$$

The matrices *P*<sup>0</sup> and *Q*<sup>0</sup> can be split in their symmetric and antisymmetric parts whereby the linearized equation of motion is as follows:

$$
\lambda \mathfrak{M} \mathfrak{o} \Delta \ddot{q} + (\mathcal{G} + \mathcal{D}) \Delta \dot{q} + (\mathcal{N} + \mathcal{K}) \Delta q = 0 \tag{17}
$$


Due to the symmetry of the bearing arrangement, it is possible to halve the system order by introducing complex coordinates. The translations are defined by

$$\mathbf{x}\_{\!\!x} = \mathbf{x} + \mathbf{j}\mathbf{y} \tag{18}$$

and the tilting by

$$
\underline{\mathcal{P}}\_{\mathcal{L}} = \mathfrak{a} - j\beta. \tag{19}
$$

The transformation *zr* = *Tc***Δ***q* leads to the system valid in the complex coordinates *z*:

$$\mathbf{M}\_{\varepsilon}\ddot{z}\_{r} + (\mathbf{G}\_{\varepsilon} + \mathbf{D}\_{\varepsilon})\dot{z}\_{r} + (\mathbf{N}\_{\varepsilon} + \mathbf{K}\_{\varepsilon})z\_{r} = 0. \tag{20}$$

The system matrices are transformed according

$$A\_{\mathcal{C}} = T\_{\mathcal{C}} \cdot A \cdot \operatorname{Re} \{ T\_{\mathcal{C}}^{\top} \},\tag{21}$$

whereby *A* stands for a general matrix.

To include the viscoelastic behavior the dynamic system model with *zr* states has to be extended with the inner states *z<sup>i</sup>* of the material model

$$\mathbf{z} = \left[ \begin{array}{c} \mathbf{z}\_r \end{array} \Big| \begin{array}{c} \mathbf{y}\_i \ \mathbf{z}\_i \end{array} \right]^T = \left[ \begin{array}{c} \mathbf{z}\_r \end{array} \Big| \begin{array}{c} \mathbf{z}\_i \end{array} \right]^T. \tag{22}$$

There *y<sup>i</sup>* are the translational and *υ<sup>i</sup>* the rotational inner states. The first order differenzial system is then given by

$$
\begin{bmatrix} E & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & M\_{\mathcal{L}} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & D\_{i} \end{bmatrix} \begin{bmatrix} \dot{z}\_{r} \\ \dot{z}\_{r}^{p} \\ \dot{z}\_{i} \end{bmatrix} + \begin{bmatrix} \mathbf{0} & \mathbf{-}E & \mathbf{0} \\ K\_{\mathcal{L}} & D\_{\mathcal{L}} + \mathbf{G}\_{\mathcal{L}} & K\_{ic} \\ K\_{ci} & \mathbf{0} & K\_{i} \end{bmatrix} \begin{bmatrix} z\_{r} \\ z\_{r}^{p} \\ z\_{i} \end{bmatrix} = \mathbf{0} \tag{23}
$$

with *<sup>z</sup>*˙*<sup>r</sup>* = *<sup>z</sup><sup>p</sup> <sup>r</sup>* . This transformation in a first order system is not trivial, because the inner states are massless and thus the mass matrix *M* is not invertible. As the inner states are only present in the first derivative only the general states *z<sup>r</sup>* have to be transformed. The system of the inner states can be included in the last row and column of (23).

With the derived system the deflections of the system as well as the eigenfrequencies can be calculated. For the eigenmodes the linear differential equation

$$
\dot{\mathbf{x}} = \mathbf{A}\mathbf{x} \tag{24}
$$

can be solved using an exponential approach

$$\mathfrak{x} = \mathfrak{k}\mathfrak{e}^{\lambda t} \tag{25}$$

leading to the characteristic system

$$\left[\lambda E - A\right]\hat{\mathbf{x}} = \mathbf{0}.\tag{26}$$

This system has a non-trivial solution when

$$\det\left(\lambda E - A\right) = 0.\tag{27}$$

The eigenvalues resulting from the solution of the polynomials have thereby the shape

$$
\lambda\_i = \delta\_i + j\omega\_i. \tag{28}
$$

*δ<sup>i</sup>* describes the damping constant and *ω<sup>i</sup>* the frequency of the corresponding eigenmode. The resulting motions of the system in its eigenmodes are depicted in Figure 8.

The green points represent the PMB's and the viscoelastic damping connections. The red cylinder is thereby showing the rotor and the blue shaft indicates the flexible stator beam connected to the grey rigid part of the stator. As it is shown in the four modes the stator is tilting as a entity, where as the rotor shows a tilting around the lower PMB in the first mode and around the upper PMB in the second mode. In the third rigid mode the rotor points a translational movement. The fourth mode indicates the semibeam bending of the flexible stator beam.

Figure 9 shows the effect of the flexible stator shaft of a designated system for different shaft materials. In the top PMB the deflections in the first resonance at about 40 Hz are increasing with softer material. Furthermore, the bending frequency, which is present in the range of 400–500 Hz is situated in the operation range. The stator shaft material has a huge influence on the relative rotor deflections in the first bending mode. Due to the required space limitation for the AMB, the motor and the position sensor system, the minimum stator shaft length is given with 116 mm as used in Figure 9. As the deflection, due to the bending, cannot be influenced by the stiffness of the PMBs and the stiffness and positioning of the damping elements, the only suitable material for the shaft, that meets the requirements of a maximum deflection in the PMB planes of 500 μm, is steel.

**Figure 8.** Schematic representation of the system motions in the three rigid body modes (**a**–**c**) and the first flexible mode of the stator beam (**d**).

**Figure 9.** Relative rotor deflection at the upper PMB of not optimized systems with different shaft materials.

It is our goal to design a system that does not exceed a maximum deflection of 500 μm in both PMBs to prevent contact with the fault bearings in case of a process-unbalance of 20 gmm. To achieve this the stiffness and the position of the PMBs and the damping elements are varied for optimization.

## **5. Optimization**

In the first optimization a given system setup was used and only the position and dimensions of the viscoelastic damping elements were optimized by minimizing the relative rotor deflections to the stator in the PMB planes. For optimization the software tool "SymSpace", developed at the Linz Center of Mechatronics (LCM), was used. It uses a genetic algorithm described in detail in [28]. In Table 1 the main parameters of the system setup are given. The mass of the rotor includes the additional mass which is mounted on the rotor. The variable parameters of the optimization with their referring range limits are shown in Table 2.


**Table 1.** Main parameters of the investigated system.



For excitation, two different unbalances with an unbalance value of 20 gmm are assumed. They are located at the top and the bottom of the rotors linked part. Finally, it turned out that the employment of two different damping materials leads to the best results. For the bottom viscoelastic element the material with a hardness of Shore A40, which was measured, is preferable. For the upper damping viscoelastic element a softer material with a hardness of Shore A20 shows the best results. However, for this material only the manufacturer data was available leading to an uncertainty since the material parameters play a major role for the rotor dynamic behavior. In Figure 10 the Pareto front of the optimization is shown.

It can be seen that the relative rotor deflections can be minimized to a value of 160 μm for both unbalance positions, which is much smaller than the PMB air gap of 500 μm. This proves the high potential of elastomer materials. The optimized values for the damping elements are summarized in Table 3.

Thereby, the lower harder element has the positive effect that it will carry the axial pre-load generated from the system.

**Figure 10.** Pareto-front of the optimization referring the maximum deflections arising due to the specified unbalances of 20 gmm.


**Table 3.** Optimized Variables.

## **6. Verification of the System Model**

To verify the system model a prototype was built with the optimized damping elements. To ensure a comparable initial situation to the simulation the rotor was initially balanced without any linked part. As the unbalance of the process mass is not predictable, the verification is done with a system without any process mass. Afterwards, two defined unbalances given in Table 4 were placed at two specified positions also used for the balancing of the system. As it is not easily possible to measure the relative deflections directly, only the absolute rotor movement is compared. The measurements are conducted using laser triangulation sensors with a sample time up to 100 kHz and a resolution of 25 nm. Hence, also high rotational speeds can be quantified with high resolution.


To eliminate the effects of the remaining rotor unbalance, the deflections of the balanced rotor were subtracted from the deflections with the defined unbalances. It should be noted that this subtraction has to be performed properly considering the phases of the measured data and is only permitted assuming a linear behavior of the system, which is already assumed in the modeling.

In Figure 11 the measured absolute rotor deflections in the PMB planes (see Figure 1) over to the rotational rotor frequency is shown. It can be noticed that the highest displacement is occuring in the lower PMB due to the unbalance in the upper balance plane. The difference to the low optimized deflections shown in Figure 10 is caused by the lower rotor mass of the verification system as it is measured without the additional mass, which has an negative influence on the system adjustment as shown in [20].

**Figure 11.** Measured deflections of the system according to applied unbalances as given in Table 4.

Thus, the measured rotor displacements need to be compared with an adapted simulation. For this the rotor has to be modified to characterize the one in the measured system without the process mass. Without the process mass the rotor features the data given in Table 5.

**Table 5.** Parameters of the varified rotor.


Comparing the adapted simulation with the measurement, as shown in Figure 12, revealed that the rigid body modes of the real system are situated at lower frequencies. That indicates a system featuring a lower overall stiffness. To identify the source of the difference first the stiffness of the PMBs are investigated. In a separated measurement of the PMBs it turned out, that the real stiffness is 20% lower than calculated. This is due to a lower PM remanent flux density, which was also confirmed with a dipole measurement of the PM rings.

**Figure 12.** Comparison of the measured data with the simulated deflections with "Unbalance1" (**a**) and "Unbalance2" (**b**) specified in Table 4.

Furthermore, the difference in the Eigenfrequencies indicate that the material of the upper damping element offers a lower stiffness, either. To verify this assumption measurements of the used material with a hardness of Shore A20 at a temperature of 20 ◦C and 25 ◦C were performed and compared to the related stiffness and damping factor characteristics provided by the manufacturer.

Figures 13 and 14 show that the real stiffness and damping behavior is smaller than expected. The measured stiffness values is only be reflected when setting the manufacturer data to a 5 ◦C higher temperature. Regarding the loss factor of the material even about 12 ◦C temperature rise has to be performed. That result shows again the necessity for proper material data, since the damping directly affects the occurring deflections of the rotor. The large deviation is also due to the measuring method used by the manufacturer, which uses a prestressing of the material in comparison to the self-used measurement process without prestressing.

**Figure 13.** Comparison of measured probe stiffness with the manufacturer data.

**Figure 14.** Comparison of measured loss factor with the manufacturer data.

In Figure 15 the comparison of the adjusted simulation is presented. The adopted temperature was set to a 12 ◦C higher temperature as the loss factor is the most important parameter considering the deflection amplitudes. The frequency range is limited to 200 Hz since the bending of the stator shaft is hardly apparent. The deflections show a very good agreement with the measurements.

**Figure 15.** Comparison of the measured data with the simulated deflections with "Unbalance1" (**a**) and "Unbalance2" (**b**) specified in Table 4 considering the lower PMB stiffnesses and a adopted temperature of 37 ◦C for the upper damping element.

## **7. Discussion**

Bringing magnetic bearing technology to a broader field of application requires the reduction of complexity and costs. The idea to use passive magnetic bearings, to achieve this goal, seems promising. So in the best case only one degree of freedom has to be controlled actively what saves costs for required actors, sensors and power electronics. Moreover, passive permanent magnetic bearings feature high reliability in case of power loss and failures. Disadvantageous proves the fact, that these type of bearings require a precise design regarding their field of application. Based on the nearly undamped system behavior, unbalances play a critical role. Different to active magnetic bearings no damping can be induced and adjusted by control, which is one of the main reasons for the limited use in commercial products. However, there are many ways to induce damping into such systems. To keep the cost advantage of passive magnetic bearings, a solution for damping must be found which only causes low additional costs. Using viscoelastic damping elements often failed in the past due to the very complex modeling and the hardly available and reliable manufacturers data. Even if available, the data is strongly influenced by the measurement process. As shown in this work, useful results can only be obtained when the exact characteristics of the damping material is well known. Unfortunately its determination is very time-consuming and an expensive task.

Moreover, a proper model of the rotordynamic behavior has to be derived to evaluate the occurring deflections of the system parts. Thereby, the most important deflections, which has to be considered, are the relative deflections in the PMB's. Since the air gap between the outer rings and the inner rings is limited, it has to be ensured, that any contact is prohibited. Thus, the generation of a proper all-over system model is complex. As there are numerous parameters influencing the system performance, the demand for an optimization is obvious. As is shown, the right configuration leads to very good results according the stabilizability of unbalances. However, only the knowledge of the exact system parameters guarantee a proper prediction of the dependent variables.

**Author Contributions:** Conceptualization, J.P.; validation, J.P.; formal analysis, J.P.; investigation, J.P.; writing—original draft preparation, J.P.; measurements, J.P.; writing—review and editing, G.J.; supervision, G.J.; funding acquisition, W.A.

**Funding:** This research received no external funding.

**Acknowledgments:** The presented research work has been supported by the Linz Center of Mechatronics GmbH (LCM), which is part of the COMET/K2 program of the Federal Ministry of Transport, Innovation and Technology and the Federal Ministry of Economics and Labor of Austria. The authors would like to thank the Austrian and Upper Austrian Government for their support.

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

## **References**


© 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/).

*Article*
