Design and Modeling of an Integrated Flywheel Magnetic Suspension for Kinetic Energy Storage Systems

The paper presents a novel configuration of an axial hybrid magnetic bearing (AHMB) for the suspension of steel flywheels applied in power-intensive energy storage systems. The combination of a permanent magnet (PM) with excited coil enables one to reduce the power consumption, to limit the system volume, and to apply an effective control in the presence of several types of disturbances. The electromagnetic design of the AHMB parts is carried out by parametric finite element analyses with the purpose to optimize the force performances as well as the winding inductance affecting the electrical supply rating and control capability. Such investigation considers both the temperature dependence of the PM properties and the magnetic saturation effects. The electrical parameters and the force characteristics are then implemented in a control scheme, reproducing the electromechanical behavior of the AHMB-flywheel system. The parameter tuning of the controllers is executed by a Matlab/Simulink code, examining the instantaneous profiles of both the air-gap length and the winding ampere-turns. The results of different dynamic tests are presented, evidencing the smooth air-gap changes and the optimized coil utilization, which are desirable features for a safe and efficient flywheel energy storage.


Introduction
The higher and higher penetration of renewable energy sources has significantly contributed to making more troublesome the power management in the electric networks. Indeed, the intermittent and fluctuating nature of the electrical generation calls for a more frequent adoption of energy storage (ES) devices, having the capability to balance both short-time and long-period mismatch between supplied and demanded power. To such purpose, the choice of a proper ES technology must be tuned according to typical power requests to avoid an abrupt performance and lifetime decay. The solution that seemingly overcomes such an issue consists of the combination of energy-intensive and power-intensive ES sources. A viable arrangement proposed by several contributions deals with kinetic-electrochemical hybrid systems, which can join well-assessed and scalable battery technology with the flywheel ES systems (FESS) able to provide back-up supply and effective peak-shaving function as well [1].
Indeed, the FESS prompt dynamics and operating flexibility can support several grid services, such as frequency regulation, ride-through capability and power quality [2,3]. Even for industry applications, FESS provides a favorable contribution to mitigate peak power requests and support conventional uninterruptable power supply (UPS) devices. FESS is potentially suitable for the residential sector for short-time blackouts, which can affect weak electrical grids [4]. Many studies and practical experiences recognize significant FESS features, such as the long life cycle, the moderate dependence of performance on environmental aspects, the fast switching between charge/discharge operations, the easy and reliable state of charge monitoring, and the suitability to realize modular schemes to fit a wide range of application [5,6].
The need for high-speed operation calls for advanced components ensuring enhanced reliability, compact design, and efficient operation, both at idle and during the charge and discharge cycles. The quest for higher efficiency aims towards a double target: on one side, higher energy is available for the power exchange; on the other, thermal stresses are more limited, avoiding the performance worsening of the electromagnetic converter parts. This issue is particularly relevant due to the very low pressure generally realized in the FESS enclosure.
On this matter, a completely bearingless system leads to a significant efficiency enhancement due to the friction loss removal, which at the same time results in reduced maintenance, improving reliability, lifetime, and cost savings, especially at a very high speed. Despite the low level of force densities, active magnetic bearings (AMBs) can control the flywheel position with a good dynamic response, limiting vibration and noise level [7,8]. AMBs are generally based on pairs of supplied coils wound around ferromagnetic cores to control either the radial or the axial position. Their action can be integrated with passive magnetic bearings (MBs; PMBs), based on permanent magnets (PMs) with or without ferromagnetic cores or on eddy current effects associated with the interaction between PMs and metallic conductors or superconductors [9]. Despite the absence of supplied coils, PMBs are, however, not free from power consumption (eddy current and core losses, cryogenic system supply) or from unstable conditions in one direction (e.g., PM-to-PM repulsive configuration).
The radial AMB magnetic configuration remarkably affects the cost, the driving circuit scheme, the control strategy, and the efficiency. As an example, Figure 1a shows a six-pole configuration, which is particularly convenient, as it enables the usage of a three-phase converter avoiding crosscoupling effect between the controlled axes [10]. However, because of heteropolar design, the power consumption is relevant due to the rotor losses determined by the flux variations along the stator circumference. This effect can be reduced using a homopolar design (Figure 1b), where there is no change of the flux polarity. In addition, they can control two degrees of freedom, reducing the number of power amplifiers [11]. This can benefit mostly the axial AMBs, which are demanded to generate high flux density to support the FESS rotor mass. Additional advantages derive from hybrid configurations provided by both PMs and supplied coils. The PM generates the static suspension force to keep the rated flywheel position, whereas the coils provide the supplementary force components in case of force unbalance related to a parameter's variation or other disturbances.  Different configurations of such MB types are proposed. In [12], there are two E-shaped stator cores hosting two PM rings and two control coils in suitable slots. The rotor consists of a plane disc sandwiched between the two stator iron cores. Such an arrangement enables symmetrical structure and the presence of the PM bias flux for suspension, but the configuration is rather complex, and simultaneous coil excitation is required in the presence of air-gap variations. A similar configuration is analyzed in [13], where the PMs are mounted externally with respect to the coils, slightly simplifying the magnetic circuit. However, high PM reluctance can limit the flux control capability, and the MB body seems rather bulky, affecting the flywheel rotor design.
In [14], two sets of axial PMB are used for the vertical suspension, whereas two sets of radial hybrid MB actively control the rotor at the bearing center. The major benefit is that the energy consumption for levitation is eliminated, limiting at the same time the PM amount by using Halbach magnetized arrays. Even in this case, the PMs placed on the rotor side require adequate room for their arrangement, fitting with particular rotor shapes as the ring-type one considered by the authors.
The paper proposes a novel axial hybrid AMB (AHMB) for the FESS axial suspension shown in Figure 2. Its distinctive feature with respect to other configurations relies on the integration with the structure of the steel rotor of the flywheel. Indeed, the outer flywheel rim is exploited as the secondary core of a couple of axial symmetric MBs. Both include ring control coils needed to stabilize the PM attraction force in presence of air-gap length variations. The upper side develops the force supporting the flywheel mass by a PM inserted between a pair of ferromagnetic cores. The two upper coils (UCs) are supplied to boost the PM attractive force as the air-gap increases (∆g < 0) with respect to the value g * , at which the PM by itself supports the flywheel mass (static condition). As they are series-connected, only one power amplifier is needed. The lower side consists of an electromagnet with a single coil (LC), which is supplied when ∆g > 0. Different configurations of such MB types are proposed. In [12], there are two E-shaped stator cores hosting two PM rings and two control coils in suitable slots. The rotor consists of a plane disc sandwiched between the two stator iron cores. Such an arrangement enables symmetrical structure and the presence of the PM bias flux for suspension, but the configuration is rather complex, and simultaneous coil excitation is required in the presence of air-gap variations. A similar configuration is analyzed in [13], where the PMs are mounted externally with respect to the coils, slightly simplifying the magnetic circuit. However, high PM reluctance can limit the flux control capability, and the MB body seems rather bulky, affecting the flywheel rotor design.
In [14], two sets of axial PMB are used for the vertical suspension, whereas two sets of radial hybrid MB actively control the rotor at the bearing center. The major benefit is that the energy consumption for levitation is eliminated, limiting at the same time the PM amount by using Halbach magnetized arrays. Even in this case, the PMs placed on the rotor side require adequate room for their arrangement, fitting with particular rotor shapes as the ring-type one considered by the authors.
The paper proposes a novel axial hybrid AMB (AHMB) for the FESS axial suspension shown in Figure 2. Its distinctive feature with respect to other configurations relies on the integration with the structure of the steel rotor of the flywheel. Indeed, the outer flywheel rim is exploited as the secondary core of a couple of axial symmetric MBs. Both include ring control coils needed to stabilize the PM attraction force in presence of air-gap length variations. The upper side develops the force supporting the flywheel mass by a PM inserted between a pair of ferromagnetic cores. The two upper coils (UCs) are supplied to boost the PM attractive force as the air-gap increases (Δ 0) with respect to the value * , at which the PM by itself supports the flywheel mass (static condition). As they are seriesconnected, only one power amplifier is needed. The lower side consists of an electromagnet with a single coil (LC), which is supplied when Δ 0. The main features of the proposed AHMB are as follows: • There are very limited winding losses because the steady-state suspension force is always applied by the PM flux; • There are no flux variations in static conditions, resulting in negligible rotor losses and the possibility to adopt solid magnetic cores instead of laminated stacks; • Possible eddy currents generated at dynamic condition provide a damping effect, improving the stability control; • The flywheel itself provides a high surface available for the axial force application, without introducing additional rotating cores; • There is independent control for the UCs and the LC; The main features of the proposed AHMB are as follows: • There are very limited winding losses because the steady-state suspension force is always applied by the PM flux; • There are no flux variations in static conditions, resulting in negligible rotor losses and the possibility to adopt solid magnetic cores instead of laminated stacks; • Possible eddy currents generated at dynamic condition provide a damping effect, improving the stability control; • The flywheel itself provides a high surface available for the axial force application, without introducing additional rotating cores; There is an inherent self-centering effect in the radial direction, possibly controlled by both the coil currents.
The latter property can be used to develop a two-degrees control or to supplement the action of dedicated radial AMBs, which are, however, needed for the full stability control of a five-degrees-of-freedom system. An auxiliary mechanical bearing (backup bearing) at the bottom end is, however, needed to support the rotor in the stationary condition or to enable a safe shutdown during a system fault. Both radial MBs and the auxiliary bearing are not represented in Figure 2, as they must fulfill typical arrangements adopted for different FESS applications [5,8]. It is worth pointing out that the configuration in Figure 2 is more convenient than the one analyzed in [15], because the use of a single PM simplifies its installation, and the double coil assembly allows a better magnetic utilization of the upper cores.
After the definition of the main FESS characteristics, the paper presents the procedure adopted for the AHMB electromagnetic design, mainly focused on the selection of the pole shoe geometry optimizing the axial force and the winding inductance. Such a procedure aims at proving the ability of the proposed AHMB to keep the flywheel suspension, complying with predefined air-gap and radial size constraints.
Then, a mathematical model is determined for the optimized configuration, consisting of the force and winding inductance representation as functions of the air-gap length and the coil current. The dependence on temperature is also considered, as it affects both the PM and the coil parameters. Finally, a bearing control scheme using the AHMB model is set up for the tuning of the Proportional-Integral-Derivative (PID) controller. Typical disturbances during FESS operation are investigated to obtain smooth instantaneous variations and zero ohmic losses at the stable condition. Different tests are provided and discussed in the final section, including the simulation of a seismic event to check the control operation in harsh conditions. The study does not consider eccentricity or radial displacement effects, assuming that some other MBs (generally the radial type ones) provide the needed compensations.

System Description
The AHMB design is developed with reference to a steel FESS based on a truncated De Laval disc with an outer rim ( Figure 2). This is a very common solution for low-and medium-speed kinetic storage, as it enables one to achieve a constant stress profile with very high shape factor K, which in turn increases the energy density [16]. The outer ring adds useful mass to the flywheel, and it can be sized to offer enough surface for the placement of an axial bearing without penalizing the stress distribution in the material. The material is an AISI 4340 steel, which presents very good mechanical properties (high ratio between ultimate stress and density σ lim /ρ). The magnetic properties are also acceptable as confirmed by the B(H) characteristic in Figure 3. storage, as it enables one to achieve a constant stress profile with very high shape factor , which in turn increases the energy density [16]. The outer ring adds useful mass to the flywheel, and it can be sized to offer enough surface for the placement of an axial bearing without penalizing the stress distribution in the material. The material is an AISI 4340 steel, which presents very good mechanical properties (high ratio between ultimate stress and density / ). The magnetic properties are also acceptable as confirmed by the characteristic in Figure 3.  Given the rated energy amount, the rotational inertia J F is given by: With ω max , ω min being maximum and minimum flywheel angular speed, respectively. The quantity s 2 − 1 represents the discharge factor, i.e., the ratio between the usable and the total stored energy. Table 1 reports the main FESS data as well as the material properties of both the flywheel rotor and the AHMB. The FESS data concern a small-scale application (i.e., UPS service or peak shaving function for a low-power renewable source), where the energy discharge must occur in a few minutes. The usage of steel flywheel in this application field has been proposed for some time now [17], and they find wide availability in the market [4]. Moreover, durability and cost issues are also factors that are evaluated in comparison with other more performant materials. The flywheel design is defined according to the procedure described in [16], assuming σ f = 1260 MPa. The additional mass is related to the rotor of the electric machine, which couples the FESS to the electric grid. The PM attraction force to suspend the flywheel mass must be F * z = 1210 N. In order to address the choice of the PM property, a preliminary rough estimation of the requested average value of the air-gap flux density is carried out using the magnetic pressure formulation, that is: According to (2), B g varies from 0.25 T to 0.3 T as k su varies from 0.7 to 0.5. Such values are hardly achievable by low-grade PMs (i.e., Ferrites) as:

•
The rated air-gap length must be large enough to deal with high-rate disturbances, the compensation of which could require high electrical and thermal stresses to the AHMB coils; Energies 2020, 13, 847 6 of 22 • Temperature increase deteriorates the PM properties with possible loss of the stable condition; • The lower the coercivity, the thicker the PM must be, therefore decreasing the useful steel surface; the consequent higher B g can be obtained only with higher PM and core volumes.
Therefore, a rare earth PM material is considered having the properties shown in Table 1. The coefficient α T describes the influence of temperature of B r and H c . The ferromagnetic cores are assumed solid with the properties of a low carbon steel AISI 1008. Its B(H) characteristic is compared to AISI 4340 steel, evidencing the higher permeability and saturation flux density. With respect to the use of laminations or soft magnetic composites, the adoption of a solid core ensures adequate mechanical stiffness, easier manufacturing, and lower production costs, without detriment to the magnetic properties.

AHMB Electromagnetic Design
The design procedure is carried out by magneto-static finite element analyses (FEAs) using a commercial code (Maxwell ® 2D, Ansys, Inc., Canonsburg, PA, USA). The purpose is to define the geometric configuration by parametric analyses, focused on the improvement of the force performances and the winding inductance values. The potential high number of geometric variables can be reduced by adopting some fair design assumptions. With reference to Figure 4, the following conditions are defined: (d) equal annular cross section of the core legs S ul,h = S ul = k u S up (h = 1, . . . , 4) for the upper side and S ll,j = S ll = k l S lp (j = 1,2) for the lower side; k u and k l are decreasing factors to enable enough room for the coil placement; their choice must consider local magnetic saturation, even with the coil mmf contribution; (e) equal cross section S uc = h uc w uc = h uc w uc of the upper coils considering that they should provide the same mmf with the same current density; (f) PM operation near the maximum energy point (little higher than B r /2) at the rated air-gap length with unexcited coils; (g) fixed variation range for the PM height h m,min , h m,max and radius r mu,min , r mu,max for manufacturing and cost issues. (a) radial bounds defined by and ; (b) equal annular cross section for all the pole shoes, that are = , (h = 1, …, 4) for the upper side and = , (j = 1, 2) for the lower side; such a condition should ensure uniform flux density distribution, resulting in evenly distributed axial force among the pole shoes on the same side; are decreasing factors to enable enough room for the coil placement; their choice must consider local magnetic saturation, even with the coil mmf contribution; (e) equal cross section = ℎ = ℎ of the upper coils considering that they should provide the same mmf with the same current density; (f) PM operation near the maximum energy point (little higher than /2) at the rated air-gap length with unexcited coils; (g) fixed variation range for the PM height ℎ , , ℎ , and radius , , , for manufacturing and cost issues.

Upper Side
The above conditions (a)-(g) lead to the definition of five independent variables related to the pole shoes' coordinates and to the PM sizes. This condition would require too wide an exploration,

Upper Side
The above conditions (a)-(g) lead to the definition of five independent variables related to the pole shoes' coordinates and to the PM sizes. This condition would require too wide an exploration, especially if more performance parameters are to be evaluated with and without excited coils. Therefore, the design of the AHMB upper side has been developed in two steps. The former determines the PM geometry to achieve a predefined air-gap flux density, which fulfills the force target and avoids excessive magnetic saturation in the core yokes. The latter performs a parametric analysis to determine the more convenient widths of the core poles to optimize some performances related to the axial forces, as well as to the winding inductances.
The preliminary FEAs vary the PM sizes in the ranges h m,min , h m,max = {2 mm, 10 mm} and r mu,min , r mu,max = {150 mm, 190 mm} to achieve the suspension force target F * z at the rated air-gap g * and operating temperature θ * . The coil section is fixed considering limits related to the maximum current density uM , winding inductance increase (more relevant with higher height-to-width ratio), and overall geometric dimensions. The other variables are simply adjusted to fulfill the condition (f). Table 2 gives the main data of the most favorable configuration for the flux density distribution. Figure 5 shows the flux density in the cores and in the air-gap, respectively.   Figure 5a shows the nearly uniform distribution in the cores with values low enough to avoid high saturation when UCs are excited. Even the air-gap flux density in front of the pole shoes meets the target value ( Figure 5b). It is worth noting that a remarkable contribution to the magnetic pressure is also present in front of the slots, which can hardly be evaluated by analytical formulations.
The selection of the current polarity in the two coils is carried out by calculating the force-current characteristic at an augmented air-gap length , = * Δ . This approach aims to evaluate the current density when the winding exerts the maximum upward force. In this case, with as the coil net section, the resultant current density must fulfill the maximum value . Figure 6 shows the four possible cases. The flux direction highlighted in the figure denotes that concordant current polarities tend to strengthen the PM flux on the outer poles (A) or on the inner ones (B); differently, discordant current polarities (C, D) generate a mixed condition of the preceding ones.  Figure 5a shows the nearly uniform distribution in the cores with values low enough to avoid high saturation when UCs are excited. Even the air-gap flux density in front of the pole shoes meets the target value ( Figure 5b). It is worth noting that a remarkable contribution to the magnetic pressure is also present in front of the slots, which can hardly be evaluated by analytical formulations.
The selection of the current polarity in the two coils is carried out by calculating the force-current characteristic at an augmented air-gap length g u,max = g * + ∆g M . This approach aims to evaluate the current density u when the winding exerts the maximum upward force. In this case, with k f S uc as the coil net section, the resultant current density must fulfill the maximum value uM . Figure 6 shows the four possible cases. The flux direction highlighted in the figure denotes that concordant current polarities tend to strengthen the PM flux on the outer poles (A) or on the inner ones (B); differently, discordant current polarities (C, D) generate a mixed condition of the preceding ones.  Figure 7 shows the resulting axial force as a function of the current density. The supply strategy A is clearly the most efficient, as it allows for the target force achievement at the lowest current density ( ≅ , corresponding to a total ampere-turns = 460 A for each coil).   Figure 7 shows the resulting axial force F z as a function of the current density. The supply strategy A is clearly the most efficient, as it allows for the target force achievement at the lowest current density ( u uM , corresponding to a total ampere-turns A u = 460 A for each coil).
Energies 2020, 13, 847 9 of 22 polarities; the orange and black arrows indicate the permanent magnet (PM) and coil flux directions, respectively. Figure 7 shows the resulting axial force as a function of the current density. The supply strategy A is clearly the most efficient, as it allows for the target force achievement at the lowest current density ( ≅ , corresponding to a total ampere-turns = 460 A for each coil). Furthermore, the force profile is almost linear, which can be useful for the design of the control system. The strategies C and D behave alike, with unacceptable current densities (≅15 A/mm 2 ). It is worth noting that strategy B provides a demagnetizing effect with a minimum for ≅ 8 A/mm 2 .
Though B is clearly not convenient for the augmented air-gap, the switching from A to B by the coil current inversion could be profitable to weaken attractive force if the flywheel should approach the AHMB upper side Δ 0 . In order to verify such a possibility, the force values are recalculated using the strategy B at the minimum upper air-gap value , = * − Δ (conf. B1). The results (strategy B1 in Figure 6) evidence that the AHMB would develop values too high to balance the suspended mass whatever the current density, therefore making inapplicable such a supply strategy.
As for the design of the magnetic core, the pole shoe geometries remarkably affect both the airgap flux density profile and the value of the coil inductance that in turn determines the voltage rating and the promptness of the supply circuit. Larger poles provide a lower magnet volume and reduced magnetic saturation, but they generally increase the leakage coil flux. Therefore, a more extensive parametric analysis is needed to explore different combinations of the pole geometric variables, keeping unvaried PM sizes and coil sections.
With reference to the geometric parametrization shown in Figure 4a, the radial coordinates of the pole shoes are expressed by the per unit (p.u.) variables , , . First, the radial coordinates of the endpoints are defined, by parametrizing Δ as a function of 0 1 : Figure 7. Force comparison between the coil supply strategies; A, B, C, D: excitation polarity represented in Figure 6; B1: excitation polarity as in B at a reduced air-gap length.
Furthermore, the force profile is almost linear, which can be useful for the design of the control system. The strategies C and D behave alike, with unacceptable current densities ( 15 A/mm 2 ). It is worth noting that strategy B provides a demagnetizing effect with a minimum for u 8 A/mm 2 .
Though B is clearly not convenient for the augmented air-gap, the switching from A to B by the coil current inversion could be profitable to weaken attractive force if the flywheel should approach the AHMB upper side (∆g > 0). In order to verify such a possibility, the force values are recalculated using the strategy B at the minimum upper air-gap value g u,min = g * − ∆g M (conf. B1). The results (strategy B1 in Figure 6) evidence that the AHMB would develop values too high to balance the suspended mass whatever the current density, therefore making inapplicable such a supply strategy.
As for the design of the magnetic core, the pole shoe geometries remarkably affect both the air-gap flux density profile and the value of the coil inductance that in turn determines the voltage rating and the promptness of the supply circuit. Larger poles provide a lower magnet volume and reduced magnetic saturation, but they generally increase the leakage coil flux. Therefore, a more extensive parametric analysis is needed to explore different combinations of the pole geometric variables, keeping unvaried PM sizes and coil sections.
With reference to the geometric parametrization shown in Figure 4a, the radial coordinates of the pole shoes are expressed by the per unit (p.u.) variables a 1u , a 2u , a 3u . First, the radial coordinates of the endpoints are defined, by parametrizing ∆r 1u as a function of a 1u (0 < a 1u < 1): Then, the following conditions are derived by imposing that the annular areas are equal and that no overlapping occurs between the middle pole shoes: with b 1u = a 2u a 3u ∆r and b 3u = (1 − a 2u )a 3u ∆r. By the elaboration of (4), the quantity ∆r is therefore obtained as: The remaining parameters to define the pole shoe geometry are finally calculated as: The quantities considered for the performance evaluation are: • the suspension force F * z0 at the air-gap g * due to the only PM; • the suspension forces F z0 and F z with augmented air-gap g u,max , due to the PM only and in the presence of current excitation (current density uM ), respectively; • the coil per-turn inductance L uc in the same condition as F z is calculated; because of the series connection and the chosen current polarity, it is L uc = L uc + L uc + 2M uc , with L uc , L uc being the inner and outer coil inductances, respectively, and M uc being the mutual inductance coefficient.
For the purpose of obtaining a homogenous comparison, the following p.u. quantities are defined for each k-th configuration (k = 1, . . . , n con f ): The chosen formulation aims at finding the optimal configuration as the one minimizing all the quantities in (7). For instance, the parameter ∆ f z is related to the force F z with coil excitation, which must be as high as possible for a given current. Figure 8a compares the quantities (7) for n con f = 126 configurations defined in the ranges a 1u = 0.65 ÷ 0.75, a 2u = 0.3 ÷ 0.9, and a 3u = 0.05 ÷ 0.3, considering as reference configuration the one obtained from the preliminary analysis. It is worth remarking the very limited improvement in the p.u. forces, while the p.u. inductance has a large variation. As no configuration approaches simultaneously the minimum value for all the quantities, a performance index P is defined to achieve a global ranking by its minimization. This index is expressed as a weighted function, which integrates the performances described by (7): Energies 2020, 13, 847 10 of 21 The coefficients < 1 ( = 1, 2, 3) give the desired importance to each performance, fulfilling the condition ∑ = 1. The subscript denotes the values for the reference configuration, in order to have = 1 as the reference one.
(a) (b) The function (8) also includes a penalty for the configurations which do not fulfill the static force reference; in this case, the unit value is increased by a contribution depending on the force deficit weighted by the coefficient .
The index is evaluated for the same = 126 configurations assuming the coefficient set , , , = 0.3, 0.4, 0.3, 100 , which is considered a good trade-off among the different The coefficients α i < 1 (i = 1, 2, 3) give the desired importance to each performance, fulfilling the condition i α i = 1. The subscript r denotes the values for the reference configuration, in order to have P = 1 as the reference one.
The function (8) also includes a penalty for the configurations which do not fulfill the static force reference; in this case, the unit value is increased by a contribution depending on the force deficit weighted by the coefficient β 1 .
The index P is evaluated for the same n con f = 126 configurations assuming the coefficient set α 1 , α 2 , α 3 , β = {0.3, 0.4, 0.3, 100}, which is considered a good trade-off among the different performances. Figure 8b highlights a limited reduction, as expected, of the index P because of the very small variations of the p.u. forces. The optimal configuration (#72) improves P only by about 2% with respect to the reference one (#60), however, confirming the effectiveness of the preliminary design. As shown in Table 3, the most evident improvement concerns the reduction of L uc by ≈ 8%, obtained mainly by increasing the slot opening of the inner core. It is worth remarking that for g u = g * , L uc 1.2 L uc and M uc 0.1 L uc ; the difference between the inner and outer coil parameters would make the force control troublesome if a parallel connection were used.  Figure 9 shows the resulting PM force with no excitation as a function of the air-gap length g u . The nonlinear dependence on g u with a steeper increase as ∆g > 0 can be noticed. This behavior will result in higher loading for the LC since it must generate a balancing force F * zl > F * zu with the air-gap variation ∆g = ∆g M being the same. However, this is not an issue because the AHMB lower side allows for an easier coil placement as no PM is present.

Lower side
The electromagnetic configuration of the AHMB lower side must be sized to provide a suitable downward attraction force on the flywheel when a disturbance tends to lower the upper air-gap below the rated value * . As observed from Figure 9, the LC loading condition is higher than the UC one for the same air-gap variation. At the PM temperature = 80 °C and Δ = Δ , the rated force for the AHMB lower part is * = * − = −580 N, with as the PM force at , ; for the

Lower side
The electromagnetic configuration of the AHMB lower side must be sized to provide a suitable downward attraction force on the flywheel when a disturbance tends to lower the upper air-gap below the rated value g * . As observed from Figure 9, the LC loading condition is higher than the UC one for the same air-gap variation. At the PM temperature θ = 80 • C and ∆g = ∆g M , the rated force for the AHMB lower part is F * zl = F * z − F z0 = −580 N, with F z0 as the PM force at g u,min ; for the AHM upper part, at g u,max , it is F * zu = 330 N. The design procedure has been developed investigating the sensitivity of some relevant geometric parameters on the axial force and the coil inductance, as well as on the flux density distribution in the core legs. Indeed, it is worth remarking that heavy saturation needs to drain more current from the electrical supply for a given winding arrangement, therefore affecting the supply ratings and the winding thermal stress. Therefore, both the coil shape and the core widths are examined as well. In order to restrict the number of variables, the design assumptions a-d are considered, including two further conditions: • constant coil section S lc with variable aspect ratio k lc = h lc /w lc ; it follows that w lc = S lc /k_lc and h lc = k lc w lc ; • variable pole shoe height h lp = h lp + h lp with fixed value for h lp .
The variable definition develops like the upper side procedure. With reference to the quantities in Figure 4c, first, the radial sizes are calculated: where a 1l is a p.u. variable describing the variation of the inner pole shoe width. By elaborating the expression of the leg cross sections as: the core leg widths can be obtained by imposing the equality of the annular areas: Therefore, by combining (9)-(11), it can be derived that the parameter set a 1l , k lc , h pl , r ml enablesone to define uniquely the geometric configuration. A preliminary analysis proves that the target force is achieved with the set a 1l , k lc , h lp , r ml = {0.95, 4.25, 6 mm, 166 mm} and current density l 6 A/mm 2 (total ampere-turns A l = 720 A). Moreover, the force characteristic is approximately linear, and the coil inductance is almost constant for l up to 8 A/mm 2 because of the limited magnetic saturation. Then, the net axial force ∆F zl = F zl − F * zl and the per-turn coil inductance L lc as functions of k lc and a 1l are examined, imposing that A l = 750 A and keeping unvaried the remaining parameters. Figure 10a shows the appreciable decrease of L lc as k lc increases, becoming nearly constant for k lc ≥ 4. On the contrary, the force varies very little, even if a tendency to decrease appears for k lc ≥ 2. Figure 10b highlights the importance of the parameter a 1l . Its value must be high enough to balance the PM attraction force (a 1l ≥ 0.8), however, leading to a noticeable L lc increase.
for up to 8 A/mm because of the limited magnetic saturation. Then, the net axial force Δ = − * and the per-turn coil inductance as functions of and are examined, imposing that = 750 A and keeping unvaried the remaining parameters. Figure 10a shows the appreciable decrease of as increases, becoming nearly constant for ≥ 4. On the contrary, the force varies very little, even if a tendency to decrease appears for ≥ 2. Figure 10b highlights the importance of the parameter . Its value must be high enough to balance the PM attraction force ( ≥ 0.8), however, leading to a noticeable increase.
(a) (b) Figure 10. Net force and inductance of the AHMB lower side; (a) dependence on the coil aspect ratio ; (b) dependence on the p.u. variable .
As a final investigation, the parameters and ℎ are jointly varied to examine possible trade-offs between force and inductance increase, as well to guarantee uniform flux density distribution inside the core legs. For such a purpose, the standard deviations Δ and Δ are calculated from the axial flux density sampled on the lines ℓ and ℓ′′, drawn in Figure 4c. The parameters vary within the ranges of 140 mm 180 mm and 5 mm ℎ 12 mm, mainly established to deal with geometrical limitations. Table 4 reports the quantities kept constant during the parametric sweep. For the sake of comparison, besides Δ , the p.u. quantities (defined as in (7), replacing with ), Δ , and Δ are evaluated. The latter two quantities are expressed as ratios of Δ and Δ with respect to their minimum value among all the examined configurations.  Figure 11 shows the progress of the selected quantities. As increases up to a certain threshold (≈156 mm, configuration #24), there is a general improvement, achieving Δ < 0 as well as low and similar flux density distributions for the two core legs. The value Δ < 0 means that the balancing condition is achieved with a lower current than in the initial configuration. From a general overview, configuration #28 represents an acceptable performance trade-off. Its design parameters are summarized in Table 4. As a final investigation, the parameters r ml and h lp are jointly varied to examine possible trade-offs between force and inductance increase, as well to guarantee uniform flux density distribution inside the core legs. For such a purpose, the standard deviations ∆B z and ∆B z are calculated from the axial flux density sampled on the lines and , drawn in Figure 4c. The parameters vary within the ranges of 140 mm ≤ r ml ≤ 180 mm and 5 mm ≤ h lp ≤ 12 mm, mainly established to deal with geometrical limitations. Table 4 reports the quantities kept constant during the parametric sweep. For the sake of comparison, besides ∆F zl , the p.u. quantities l coil (defined as in (7), replacing L uc with L lc ), ∆b z , and ∆b z are evaluated. The latter two quantities are expressed as ratios of ∆B z and ∆B z with respect to their minimum value among all the examined configurations.  Figure 11 shows the progress of the selected quantities. As r ml increases up to a certain threshold (≈156 mm, configuration #24), there is a general improvement, achieving ∆F zl < 0 as well as low and similar flux density distributions for the two core legs. The value ∆F zl < 0 means that the balancing condition is achieved with a lower current than in the initial configuration. From a general overview, configuration #28 represents an acceptable performance trade-off. Its design parameters are summarized in Table 4.
Energies 2020, 13, 847 13 of 21 Figure 11. Results of the parametric analysis on the pole shoe shape of the AHMB lower side.

Control System
The assessment of the control system and the analysis of the dynamic operation of the FESS suspension were carried out by a Matlab ® Simulink (The Mathworks, Inc., Natick, Massachusetts, USA) code. This code implements both the force balance equation governing the flywheel vertical motion and the electrical equation reproducing the coil supply condition. The solution of such equations requires, in turn, the definition of the electromagnetic models describing the axial forces and the coil inductances as functions of the air-gap length and of the coil ampere turns.

Model of the Suspension Dynamics
The resultant force in the vertical direction governs the dynamic variation of the air gap according to the following equation: The time variation of the axial coordinate is associated with the corresponding air-gap variation Δ with respect to the width * at the stable condition. It is worth noting that Δ is affected by the vertical movements of both the flywheel and the AHMB fixed parts, the latter is related to FESS enclosure oscillations (i.e., seismic activity or soil-transmitted vibrations).
The vertical force due to the PM and due to the UC are dependent on the PM temperature as well. This dependence is quite important, mostly for installations in a low-pressure environment where the heat removal is mainly based on radiation effect, making critical the thermal operation, particularly for small rated systems [18]. To deal with the temperature variation, the total force , = + of the AHMB upper side has been modeled by the combination of two terms [15]. The former 〈 , 〉 = , ( ) consists of the value related to the mean temperature = ( + )/2 with , minimum and maximum operating temperature during a FESS cycle; the latter Δ , = , ( ) − 〈 , 〉 is a differential force considering the contribution of the temperature excursion with respect to for a given current and air-gap length. The force function is formulated as: with ( ) weighting coefficient having a linear variation with temperature. Equation (13) represents an interpolating function providing the upper force at any temperature, knowing the values at and . Such values are obtained by means of a sequence of magnetostatic FEAs, and then the values 〈 , 〉 and Δ , as functions of Δ are calculated and Figure 11. Results of the parametric analysis on the pole shoe shape of the AHMB lower side.

Control System
The assessment of the control system and the analysis of the dynamic operation of the FESS suspension were carried out by a Matlab ® Simulink (The Mathworks, Inc., Natick, MA, USA) code. This code implements both the force balance equation governing the flywheel vertical motion and the electrical equation reproducing the coil supply condition. The solution of such equations requires, in turn, the definition of the electromagnetic models describing the axial forces and the coil inductances as functions of the air-gap length and of the coil ampere turns.

Model of the Suspension Dynamics
The resultant force F z in the vertical direction governs the dynamic variation of the air gap according to the following equation: The time variation of the axial coordinate z is associated with the corresponding air-gap variation ∆g with respect to the width g * at the stable condition. It is worth noting that ∆g is affected by the vertical movements of both the flywheel and the AHMB fixed parts, the latter is related to FESS enclosure oscillations (i.e., seismic activity or soil-transmitted vibrations).
The vertical force F z0 due to the PM and F zu due to the UC are dependent on the PM temperature θ as well. This dependence is quite important, mostly for installations in a low-pressure environment where the heat removal is mainly based on radiation effect, making critical the thermal operation, particularly for small rated systems [18]. To deal with the temperature variation, the total force F zu,t = F z0 + F zu of the AHMB upper side has been modeled by the combination of two terms [15]. The former F zu,t = F zu,t (θ mean ) consists of the value related to the mean temperature θ mean = (θ max + θ min )/2 with θ min , θ max minimum and maximum operating temperature during a FESS cycle; the latter ∆F zu,t = F zu,t (θ min ) − F zu,t is a differential force considering the contribution of the temperature excursion with respect to θ mean for a given current and air-gap length. The force function is formulated as: with k 0 (θ) weighting coefficient having a linear variation with temperature. Equation (13) represents an interpolating function providing the upper force at any temperature, knowing the values at θ mean and θ min . Such values are obtained by means of a sequence of magnetostatic FEAs, and then the values F zu,t and ∆F zu,t as functions of ∆g are calculated and stored in 2D lookup tables. The corresponding plots are presented in Figure 12. The curve of F zu,t at A u = 0 represents the PM force contribution, which clearly increases with A u .
Energies 2020, 13,847 14 of 21 stored in 2D lookup tables. The corresponding plots are presented in Figure 12. The curve of 〈 , 〉 at = 0 represents the PM force contribution, which clearly increases with .
(a) (b) However, at the lower air-gap values, such increase is to some extent limited by the magnetic saturation. The same effect occurs for Δ , which decreases by over 10%; however, its variation does not remarkably affect the total force value. As for the AHMB lower side, the LC force is independent of the temperature because no PM is used.

Electrical Model
The supply system of the UCs and LC controls the flywheel suspension by exciting alternatively the winding to compensate for the changes in the air-gap length. An electrical model is therefore needed to address the choice of the voltage rating and of the current regulator parameters governing the dynamic performances of the whole suspension system. The voltage equation is basically the same for both the windings, including the ohmic voltage drop as well the inductive and motional emf related to the current and air-gap variations, respectively. For the sake of an easier implementation, it can be written considering the ampere-turns as the state variable, that is: with as the supply voltage, as coil turns, and , , and as the per-turn resistance, inductance, and PM flux linkage, respectively. For the LC, it is / = 0 as no PM is present. The differential Equation (14) is clearly nonlinear with variable coefficients, also taking into account that both and are dependent on temperature. This property requires the implementation in a Simulink subsystem, the results of which are subsequently passed to the air-gap control scheme. Figure 13 shows the complete block diagram. The input signals consist of the air-gap and the vertical velocity / representing the flywheel dynamic state, the reference ampere-turns * from the air-gap control regulator, and the AHMB temperature . * is saturated to the value providing the maximum allowable current density. The current control loop governed by a saturated PI regulator adjusts the ampere-turns according to the reference value. The PI saturation limits the voltage to a maximum value even during transient conditions, complying with the supply converter rating. The scheme does not consider the converter control strategy (e.g., PWM or hysteresis control), as its time constant is generally negligible in comparison to the winding electrical one.
The activation function elaborates the current signal to activate alternatively the coil excitation. Two different expressions are used, which are ( ) = (| | + )/2 assigned to the UC control scheme and ( ) = (| | − )/2 to the LC one. In this way, according to the reference * However, at the lower air-gap values, such increase is to some extent limited by the magnetic saturation. The same effect occurs for ∆F zu,t which decreases by over 10%; however, its variation does not remarkably affect the total force value. As for the AHMB lower side, the LC force F zl is independent of the temperature because no PM is used.

Electrical Model
The supply system of the UCs and LC controls the flywheel suspension by exciting alternatively the winding to compensate for the changes in the air-gap length. An electrical model is therefore needed to address the choice of the voltage rating and of the current regulator parameters governing the dynamic performances of the whole suspension system. The voltage equation is basically the same for both the windings, including the ohmic voltage drop as well the inductive and motional emf related to the current and air-gap variations, respectively. For the sake of an easier implementation, it can be written considering the ampere-turns A as the state variable, that is: with v as the supply voltage, N as coil turns, and R, L, and λ pm as the per-turn resistance, inductance, and PM flux linkage, respectively. For the LC, it is dλ pm /dz = 0 as no PM is present. The differential Equation (14) is clearly nonlinear with variable coefficients, also taking into account that both R and λ pm are dependent on temperature. This property requires the implementation in a Simulink subsystem, the results of which are subsequently passed to the air-gap control scheme. Figure 13 shows the complete block diagram. The input signals consist of the air-gap g u and the vertical velocity dz/dt representing the flywheel dynamic state, the reference ampere-turns A * from the air-gap control regulator, and the AHMB temperature θ. A * is saturated to the value providing the maximum allowable current density.
Energies 2020, 13, 847 15 of 21 sign, UC ( * 0) or LC ( * < 0) is excited. However, a reset of the PI integrator is also needed when * = 0 to avoid simultaneous conduction of both the windings. As for the electrical parameters, the resistance is calculated using the classical formulation dependent on the winding temperature : with , as the mean radius of the upper coils (inner and outer ones, respectively), , as copper resistivity at , and thermal coefficient as Δ = − and Δ = − . The UC and LC inductances and and the PM flux are determined for different airgap and ampere-turn values by the same sequence of FEAs used for the force components. As for the dependence on temperature, is assumed to decrease linearly as temperature increases according to the material temperature coefficient; differently, is assumed independent of . The inductance curves are reported in Figure 14.
is appreciably affected by the current, confirming the influence of magnetic nonlinearity (≈−15% with = , ). The dependence on the air-gap estimated by interpolation is ~ / , deviating from the conventional relation ~ [19], likely due to magnetic saturation and leakage effects. denotes values comparable with , however, being almost insensitive to the coil current. The current control loop governed by a saturated PI regulator adjusts the ampere-turns according to the reference value. The PI saturation limits the voltage to a maximum value v max even during transient conditions, complying with the supply converter rating. The scheme does not consider the converter control strategy (e.g., PWM or hysteresis control), as its time constant is generally negligible in comparison to the winding electrical one.
The activation function G elaborates the current signal to activate alternatively the coil excitation. Two different expressions are used, which are G(u) = (|u| + u)/2 assigned to the UC control scheme and G(u) = (|u| − u)/2 to the LC one. In this way, according to the reference A * sign, UC (A * > 0) or LC (A * < 0) is excited. However, a reset of the PI integrator is also needed when A * = 0 to avoid simultaneous conduction of both the windings.
As for the electrical parameters, the resistance is calculated using the classical formulation dependent on the winding temperature θ: with r uc , r uc as the mean radius of the upper coils (inner and outer ones, respectively), ρ min , α as copper resistivity at θ min , and thermal coefficient as ∆θ uc = θ uc − θ min and ∆θ lc = θ lc − θ min . The UC and LC inductances L uc and L lc and the PM flux λ pm are determined for different air-gap and ampere-turn values by the same sequence of FEAs used for the force components. As for the dependence on temperature, λ pm is assumed to decrease linearly as temperature increases according to the material temperature coefficient; differently, L uc is assumed independent of θ. The inductance curves are reported in Figure 14. L uc is appreciably affected by the current, confirming the influence of magnetic nonlinearity (≈−15% with g u = g u,min ). The dependence on the air-gap estimated by interpolation is ∼ g −1/2 u , deviating from the conventional relation ∼ g −1 u [19], likely due to magnetic saturation and leakage effects. L lc denotes values comparable with L uc , however, being almost insensitive to the coil current.  Figure 15 shows the overall scheme of the AHMB control system, incorporating both the outer loop regulating the air-gap dynamics and the inner loops regulating the electrical supply of the UCs and the LC. The reference for the ampere-turns is generated by a PID regulator according to the airgap deviation Δ obtained from the integration of (12). The transfer function ( ) of the PID regulator has a classical formulation, including a derivative filter coefficient to limit the high spikes typical of the pure derivative action. Therefore, it is ( ) = + / + /(1 + / ).

Overall Model
By an anti-windup integration, the air-gap should be progressively adjusted to set the coil current to zero. Two separated integrators are used to deal with the different UC and LC electrical parameters. The integrator coefficients must be tuned to limit transient energy consumption, however, avoiding steady-state oscillations as well. The actual value of the upper air-gap obtained as = − Δ ( : initial value) is passed to the force and electrical models described in the preceding sections. The temperature should consider the thermal operation of both the flywheel affected by the electrical machine and the AHMB itself. The Simulink code enables the analysis of several types of external disturbances to verify the flywheel dynamic operation in the vertical direction, considering the action of the proposed AHMB control. Among them, three cases are implemented in the diagram: • an external force added to the force balance; in practical cases, such contribution can derive from vertical attraction forces between the rotor and the stator of the electrical machine;   Figure 15 shows the overall scheme of the AHMB control system, incorporating both the outer loop regulating the air-gap dynamics and the inner loops regulating the electrical supply of the UCs and the LC. The reference for the ampere-turns is generated by a PID regulator according to the air-gap deviation ∆g obtained from the integration of (12). The transfer function F PID (s) of the PID regulator has a classical formulation, including a derivative filter coefficient K N to limit the high spikes typical of the pure derivative action. Therefore, it is F PID (s) = K P + K I /s + K D K N /(1 + K N /s).  Figure 15 shows the overall scheme of the AHMB control system, incorporating both the outer loop regulating the air-gap dynamics and the inner loops regulating the electrical supply of the UCs and the LC. The reference for the ampere-turns is generated by a PID regulator according to the airgap deviation Δ obtained from the integration of (12). The transfer function ( ) of the PID regulator has a classical formulation, including a derivative filter coefficient to limit the high spikes typical of the pure derivative action. Therefore, it is ( ) = + / + /(1 + / ).

Overall Model
By an anti-windup integration, the air-gap should be progressively adjusted to set the coil current to zero. Two separated integrators are used to deal with the different UC and LC electrical parameters. The integrator coefficients must be tuned to limit transient energy consumption, however, avoiding steady-state oscillations as well. The actual value of the upper air-gap obtained as = − Δ ( : initial value) is passed to the force and electrical models described in the preceding sections. The temperature should consider the thermal operation of both the flywheel affected by the electrical machine and the AHMB itself. The Simulink code enables the analysis of several types of external disturbances to verify the flywheel dynamic operation in the vertical direction, considering the action of the proposed AHMB control. Among them, three cases are implemented in the diagram: • an external force added to the force balance; in practical cases, such contribution can derive from vertical attraction forces between the rotor and the stator of the electrical machine;  By an anti-windup integration, the air-gap should be progressively adjusted to set the coil current to zero. Two separated integrators are used to deal with the different UC and LC electrical parameters. The integrator coefficients must be tuned to limit transient energy consumption, however, avoiding steady-state oscillations as well. The actual value of the upper air-gap obtained as g u = g 0 − ∆g (g 0 : initial value) is passed to the force and electrical models described in the preceding sections. The temperature θ should consider the thermal operation of both the flywheel affected by the electrical machine and the AHMB itself.
The Simulink code enables the analysis of several types of external disturbances to verify the flywheel dynamic operation in the vertical direction, considering the action of the proposed AHMB control. Among them, three cases are implemented in the diagram: • an external force F zd added to the force balance; in practical cases, such contribution can derive from vertical attraction forces between the rotor and the stator of the electrical machine; • an air-gap variation ∆g d caused, for instance, by a sudden motion of the AHMB fixed parts; • a superimposed vertical velocity dz/dt, simulating, for instance, the effect of an earthquake (instantaneous data of a seismic velocity).
The first two cases mainly aim to tune the controller parameters and examine the effectiveness of both the regulators and the anti-windup action. As for the earthquake simulation, the purpose is not to reproduce exactly the flywheel vertical trajectory, but to check the system capability to react promptly in the presence of a high variability disturbance, keeping the air-gap excursion within the predefined limits.

Simulation Results
Some examples of application are proposed to prove the control effectiveness for various types of disturbances. The winding electrical parameters are calculated considering the coil's geometrical characteristics as well as the ratings of a typical DC supply system ( Table 5). The design current density is higher for the lower coil because of the PM absence. In addition, its resistance is lower because of the lower coil turns with higher wire diameter. The parameters of the current and air-gap regulators are set after some preliminary tests starting from the use of the built-in procedures for the regulator self-tuning implemented in the Simulink block. Differently, the anti-windup integrators are specifically tuned by considering air-gap step variation. The target was to achieve the current cancellation by 50 s with superimposed air-gap excursion ∆g = ±0.5 mm. Basically, a trade-off between two opposite targets has been investigated. On one hand, a prompt response is required to constrain the flywheel motion; on the other hand, a slower recovery of the steady state is desirable to avoid excessive electrical stresses.
The first test analyzes the system dynamics during a temperature variation caused by a sequence of charge and discharge cycles of the FESS and a subsequent idle condition. The temperature profile θ(t) is approximated by a trapezoidal profile starting from θ min and achieving θ max after 120 s.
Then, the temperature is held constant (120 s) and finally decreased to θ min again (180 s). Moreover, it is assumed that θ uc = θ lc = θ. Figure 16a presents the dynamic variation of the air-gap and of the AHMB winding ampere-turns. At first, there is a settling-down period since the initial position (2 mm) does not fulfill the force balance at θ min . Therefore, the LC is supplied to enable a recovery of the equilibrium point with no excitation (≈2.37 mm). The initial peak of the ampere-turns aims at the limitation of the air-gap derivative. During the temperature variation, the air-gap follows a smooth transition towards the steady state values, thanks to a proper excitation first of the upper and then of the lower coils. After their action, the currents are always driven to zero in a few tenths of a second, confirming the effectiveness of the anti-windup control block. The final test investigates the effects of the vibration transmitted to the flywheel enclosure by a seismic event. The superimposed vertical speed consists of the data measured by the Italian Geophysics and Volcanology Institute during the earthquake in L'Aquila (Italy) in 2009, which was about magnitude 6.0 in the Richter scale [20]. The main purpose is to assess that the control system does not amplify the oscillations and finally drives the air-gap to its steady state value in spite of such a critical phenomenon. The disturbance speed and the air-gap instantaneous profile are presented in Figure 17a. Though the air-gap denotes a high-frequency noise, the variation amplitude is moderate, and the drift tendency has been successfully compensated by the AHMB excitation control. In this case, the AHMB windings are alternatively excited (Figure 17b), with a predominant contribution of the UC mainly in the final interval to restore the steady-state air-gap length. The zoomed box in the same figure confirms that only one coil is supplied at a time, with a consequent decrease of the power consumption. Moreover, the ampere-turn values are always within the allowable limits.  The second test analyzes the response to a step force disturbance F zd = −150 N at θ = θ max as shown in Figure 16b. Even in this case, the coil current reacts promptly to compensate for the downward overload, allowing an air-gap decrease to the value (≈1.68 mm), at which the force balance is restored with no coil excitation.
The final test investigates the effects of the vibration transmitted to the flywheel enclosure by a seismic event. The superimposed vertical speed consists of the data measured by the Italian Geophysics and Volcanology Institute during the earthquake in L'Aquila (Italy) in 2009, which was about magnitude 6.0 in the Richter scale [20]. The main purpose is to assess that the control system does not amplify the oscillations and finally drives the air-gap to its steady state value in spite of such a critical phenomenon. The disturbance speed and the air-gap instantaneous profile are presented in Figure 17a. Though the air-gap denotes a high-frequency noise, the variation amplitude is moderate, and the drift tendency has been successfully compensated by the AHMB excitation control. In this case, the AHMB windings are alternatively excited (Figure 17b), with a predominant contribution of the UC mainly in the final interval to restore the steady-state air-gap length. The zoomed box in the same figure confirms that only one coil is supplied at a time, with a consequent decrease of the power consumption. Moreover, the ampere-turn values are always within the allowable limits. The final test investigates the effects of the vibration transmitted to the flywheel enclosure by a seismic event. The superimposed vertical speed consists of the data measured by the Italian Geophysics and Volcanology Institute during the earthquake in L'Aquila (Italy) in 2009, which was about magnitude 6.0 in the Richter scale [20]. The main purpose is to assess that the control system does not amplify the oscillations and finally drives the air-gap to its steady state value in spite of such a critical phenomenon. The disturbance speed and the air-gap instantaneous profile are presented in Figure 17a. Though the air-gap denotes a high-frequency noise, the variation amplitude is moderate, and the drift tendency has been successfully compensated by the AHMB excitation control. In this case, the AHMB windings are alternatively excited (Figure 17b), with a predominant contribution of the UC mainly in the final interval to restore the steady-state air-gap length. The zoomed box in the same figure confirms that only one coil is supplied at a time, with a consequent decrease of the power consumption. Moreover, the ampere-turn values are always within the allowable limits.

Conclusions
In this paper, a novel configuration of an axial hybrid magnetic bearing (AHMB) is analyzed for the suspension of steel flywheels. It enables several benefits in terms of integration with the flywheel design, control flexibility, and low electromagnetic losses that are essential to reduce power consumption in energy storage applications. By means of parametric analyses using finite element simulations, an optimized design of the magnetic bearing parts is obtained, proving that the required force is achieved without any change to the rotor shape and that the winding inductance values match voltage rating and control requirements. Moreover, the resulting configuration presents limited PM volume, compact sizes, and moderate ampere-turn values for the flywheel suspension even with large air-gap variations.
Afterward, the implementation of the AHMB nonlinear electromagnetic model in a numerical code for the dynamic analysis addressed the tuning of the regulator parameters related to the axial force control. The simulation results concerning different types of disturbances evidence the control system promptness without excessive oscillations for both the air-gap and the coil currents. In particular, the validity is demonstrated of the current control technique which manages efficiently the upper and lower bearing supply to enable separate excitation between the two circuits during transients as well as no excitation at steady-state. This behavior is confirmed also in the presence of a high-stress condition like a seismic event, indicating the suitability of the proposed AHMB for practical applications.
Author Contributions: A.T. wrote and reviewed the paper, developed the magnetic bearing design and modeling and investigated the dynamic behaviour. M.A. conceived the magnetic bearing configuration, developed the control strategy, reviewed the paper and supervised the development activities. R.B. conceived the project, acquired the founding resources, reviewed and approved the final paper. All authors have read and agreed to the published version of the manuscript.
Funding: This work was financially supported by the Department of Industrial Engineering-University of Padova under the project entitled "Flywheels and Li-ion Batteries for Stationary Hybrid Energy Storage Systems".

Conflicts of Interest:
The authors declare no conflict of interest.
List of symbols α 1 , α 2 , α 3 , β 1 weighting coefficients of the performance index P a 1u , a 2u , a 3u per-unit variables used for the parametrization of the upper pole shoes geometry a 1l per-unit variable used for the parametrization of the lower pole shoe geometry A u , A l upper and lower coil ampere-turns b 1u , b 3u slot openings of the upper cores B r , H c , α T PM remanence, coercivity, and temperature coefficient B g mean air-gap flux density (axial component) B s , µ r saturation flux density and magnetic permeability ∆B z , ∆B z standard deviations of the axial flux density in the lower core legs ∆F zl net axial force with lower coil excitation ∆g, ∆g M air-gap deviation with respect to g * and its maximum value ∆g d air-gap disturbance ∆r 1u , ∆r 3u , ∆r 5u , ∆r 7u pole shoe widths of the inner and the outer cores (upper side) ∆r 1l width of the inner pole shoe of the lower side core F z , F z0 , F zu , F zl resultant, PM, upper side and lower side force contribution F * z requested force for the total mass suspension F * z0 rated value of the PM suspension force F z0 , F z0 PM suspension forces at g * + ∆g M and g * − ∆g M ; F z upper side suspension forces at g * + ∆g M with current excitation; F * zu , F * zl PM force deviation at the maximum and minimum air-gap bounds; F zu,t , ∆F zu,t contributions of upper side force model dealing with temperature variation F zd force and air-gap disturbances g * air-gap length at balanced condition with no coil supply g u , g l upper and lower air-gap length h up , h lp height of the upper and lower pole shoes h uc , w uc height and width of the inner coil (upper side) h uc , w uc height and width of the outer coil (upper side) h lc , w lc height and width of the lower coil h m , l m , V m PM height, length and volume K, W F , J F flywheel shape factor, rated energy amount, and rotational inertia k 0 weighting coefficient to model the force with temperature variations k f , S uc , S lc coil filling factor, upper and lower coil cross sections L, R, λ pm per-turn inductance, resistance, and PM flux L uc , L lc upper and lower coil inductance M F , M a flywheel and additional masses n max , s maximum flywheel speed and speed ratio P performance index θ, ∆θ uc , ∆θ uc operating temperature and temperature coil rise ρ, σ lim , σ f mass density, ultimate and low-cycle fatigue stresses ρ min , α copper resistivity at θ min and thermal coefficient u , l upper and lower coil current density r 1u , . . . , r 6u radial coordinates of the upper pole shoes r 1l , r 2l radial coordinates of the lower pole shoes r mu mean PM radius r ml mean coil radius of the lower side coil r o , r i outer and inner radius of the flywheel rotor rim r uc , r uc mean coil radius of the upper side coils R uc , R lc upper and lower coil ohmic resistance S r , k su rim surface and surface utilization factor S up,k , S lp,k annular cross section of the k-th pole shoe (up: upper side, lp: lower side) S ul,k , S ll,k annular cross section of the k-th core leg (ul: upper side, ll: lower side) v, N supply voltage and coil turns ω max , ω min maximum and minimum flywheel angular speed w 1u , w 2u , w 1u , w 2u width of the inner and of the outer core legs (upper side) w l , w l width of the lower side core legs w 3u , w 3u thickness of the inner and outer yokes (upper side)