Design and Optimization of Synchronous Motor Using PM Halbach Arrays for Rim-Driven Counter-Rotating Pump

: This document deals with the design of a Permanent Magnet Synchronous Motor (PMSM) to peripherally drive a counter-rotating pump inducer. The motor/pump is associated using a rim-driven principle where the motor’s active parts are located at the periphery of the inducer blades. It proposes using a Halbach array of permanent magnets for the active rotor of the motor. This solution allows the generation of a Sinusoidal Electromotive Force (EMF). Therefore, a more stable electromagnetic torque is reached. An optimum geometry suitable for the inducer speciﬁcations while respecting operational constraints is determined. The obtained geometry is then simulated using the Finite Element Method. The results are satisfactory in terms of average torque and EMF waveform. Use of the Halbach array allows a signiﬁcant improvement of the ﬂux density in the air gap compared to a designed surface-mounted machine. The experimental validation will be performed once the prototype is realized in the Laboratory of Fluid Engineering and Energy systems (LISFE).


Introduction
The unstable performance and high energy consumption of the centrifugal pumps lead to the proposal of an unconventional design to improve pump behavior.The counterrotation pump is one such innovative solution.It generally consists of an inducer and an impeller, as shown in Figure 1 [1].A study by [1] shows that driving the inductor independently and in the opposite direction of the impeller enhances the stability of the pump performance.The aim of the proposed article is to propose to design a rim-driven motor to separately drive the counter-rotating inducer.Rim-driven systems have been studied in the literature for years.They are used for marine propulsion [2,3], aeronautical applications [4], tidal [2,5,6] and wind turbines [7].The rim-driven machine is characterized by a magnetic ring located around the blades of a turbine where magnets are fixed on.A magnetic stator (which contains windings and magnetic cores) surrounds this permanent magnet ring.In comparison with more classical motors where the motor structure is connected to the turbine by a shaft, the rim-driven motor has a larger diameter and lower axial length.It can also be noted that thanks to the rim-driven structure, the turbine can be shaftless.Some comparative studies [8][9][10][11][12][13] confirm that the Permanent Magnet Synchronous Motor (PMSM) is the most suitable for designing a rim-driven system, so the earlier studies aim to improve its performance by optimizing space, maximizing torque or increasing motor efficiency.In PM synchronous machines, the use of the Halbach array in the rotor provides several advantageous features.The most important advantage is the cancellation of the flux density on one side and reinforcement on the other side of the hmagnets [14].In addition, the induced electromotive force and the magnetic flux created by the Halbach array can be quasi-sinusoidal.This last feature leads to minimizing the electromagnetic torque if the motor drive controls the sinusoidal currents (flux-oriented control).This work proposes to design a rim-driven Halbach-array-based PMSM to drive the inducer of a counter-rotating pump presented in Figure 1.The inner diameter of the driving motor is equal to the tip diameter of the inducer.The rated speed of the inducer, which has been designed for this application, is 5000 rpm, and its rated torque is 10 Nm.A first optimal design of PMSM with classical surface-mounted magnets was proposed in a previous study presented in [13].In this classical design, as the back EMF was not sinusoidal, the electromagnetic torque oscillations reached 20%, which may disturb the operation of the inducer.To overcome this issue, a design of the PMSM rim-driven motor using a PM Halbach array with progressive magnetization magnets is proposed in this paper.This solution improves the magnetic flux density waveform in the machine air gap and allows minimizing the iron used in the rotor core or even to replace it by a non-magnetic material.
The first part of the article deals with the electromagnetic study of the machine, where the set of formulas that is used for the machine design is presented.The second part presents a thermal study that allows an estimation of the temperatures inside the motor active part at its rated operating point.In a third section, a Finite Element simulation of the designed machine is presented to validate the proposed design.The results for the previously designed surface-mounted magnet structure and the new progressive magnetization structure are compared to highlight the improvement related to the Halbach Array proposed solution.The main contribution of this paper is to propose a specific design method for permanent magnet motors with a Halbach array in the rotor.The motor is used to rim-drive the turbine of a high-performance pump.The pump's turbine is directly integrated into the inner diameter of the rotor, so it requires the most stable torque possible.This method is based on an optimization process of the electrical motor taking into account in a multidisciplinary approach an electromagnetic model, Rim-Driven turbine specifications and thermal and design constraints.

Electromagnetic Design
To determine the geometry of the rim-driven machine, the outer diameter of the turbine is assumed to be fixed.The main specification for the rim-driven motor is to be able to deliver the required torque to drive the pump inducer.The average value of the electromagnetic torque of a PMSM, given by Equation (1), depends on the internal diameter of the stator bore, the axial length of the machine, the maximum value of the first harmonic of the magnetic flux density created by the rotor magnets in the air gap and the type of winding used.
where k w1 is the winding factor, A L is the electric load, B 1 is the first harmonic of the magnetic flux density created in the airgap, D is the internal stator diameter, L is the active axial length and ψ is the phase shift between EMF and the phase current (in the case of maximal torque per ampere flux-oriented control of the motor, this angle is controlled to be null).The structure of the novel, progressive magnetized Halbach magnet-based rim-driven motor is given in Figure 2.Each pole is divided on n s magnets (segments).As shown in Figure 2, each of these segments is characterized by a particular magnetization direction θ m (referring to the segment axis).θ m is determined using Equation (2): where k, p are the segment numbers per pole and the pair pole number, respectively.Based on the above-mentioned characteristics of the Halbach magnet arrays, the flux in the rotor magnetic circuit can be neglected if the magnet layer thickness is not too small so that a non-magnetic material can be used instead of the iron rotor core.

Design of Halbach Array Permanent Magnets
The peak value of the first harmonic in the gap is a crucial parameter for the torque calculation, as shown in Equation (1).For a Halbach array-based structure, its expression is presented in Equation (3) [15].The use of Halbach array magnets allows us to achieve the maximum of the fundamental component of the magnetic flux density in the air gap equal to [16]: B1 = B emax .
In Equation (3), B r is the remanent flux density in the magnet, R r = D hel 2 + h ry is the rotor radius, R m = R r + h m is the external magnet radius, and γ is determined according to the material used at the back side of the Halbach array: p and l g are, respectively, the pole pair number and the gap thickness, and the expression R s = h m + R r + l g represents the internal stator core radius.These geometric parameters are shown in Figure 3.For a fixed internal diameter R r , the peak magnetic flux density B emax from Equation (1) depends on the gap thickness, the magnet height, and the pole pair number.The graph in Figure 2 represents the variation of B emax as a function of l g and h m for a pole pair number equal to three, and the dimensions and material characteristics are given in Table 1.
Table 2 shows the effect of the magnet height and airgap thickness on the peak value of the flux density in the airgap.Magnetic flux density increases with magnet thickness.Indeed, from Figure 4, the magnetic flux density tends to saturation for high values of h m .Therefore, a limit of h m is imposed to be 10 mm.   2.
Minimizing the airgap maximizes the peak value of the airgap flux density.In [17], it is shown that the airgap thickness for small rim-driven systems is a token between 1 to 2% of the internal rotor diameter.Indeed, an expression of l g is presented in [18]: where L, D are, respectively, the active length and the internal stator core diameter.

Magnetic Circuit Design
The magnetic circuit sizing consists of determining all the geometric parameters that could affect the circulation of the magnetic field.These main dimensions are shown in Figure 3 [13], where h sy , h ry are the stator and rotor yokes, h ss the slot depth, w tooth the tooth width.
The internal diameter of the stator D is given by: where D hel is the turbine external diameter, and h ry is the rotor yoke thickness

Tooth Sizing
Considering that the total flux created by magnets crosses all teeth, the tooth width could be written as: k tooth is the ratio of the tooth width to slot pitch, and h ss is the slot depth.

Yokes Thickness
The stator yoke design takes into account the maximum average flux density that crosses this part of the motor.The stator yoke thickness is given by: B sat is the saturation flux density of the used material.Assuming the rotor core is air or a non-magnetic material, the rotor yoke thickness is chosen to be smaller than the one of the stator.

Slot Sizing
The determination of stator thickness allows the calculation of other geometric parameters, taking into consideration the nozzle effect.The slot sizes are given as follows.For the airgap side: For the external diameter side: where Qs is the slot number, it generally depends on the chosen winding type and the pole pair number.D is the internal stator diameter shown in Figure 3.The slot depth h ss is given by [19]: J is the current density, and k r is the filing factor (copper section to slot section ratio).We can then deduce the slot area as follows: and the section occupied by the copper in the slot could be calculated by:

Machine Efficiency
It represents the output-to-input power ratio, and the output power is given by: where T m and Ω are, respectively, the torque provided by the motor and angular speed of the turbine, respectively.The machine efficiency is expressed as: P cop and P iron represent the copper and iron losses, and their expressions are detailed in the next section.

Thermal Sizing
A good evaluation of the heat sources and the different elements of the thermal circuit is necessary for the thermal modeling of the electric machine.This modeling allows an estimation of the temperatures inside the machine and especially in sensitive areas such as the stator winding.The heat in electrical machinery originates from iron losses, friction, or Joule losses.

Iron Losses
The losses in the iron depend on the variation in the flux density that crosses the magnetic circuit, the electrical frequency, the volume of the used iron, and finally, the iron characteristics.An approximative relation can be used to determine those losses using the mass power density of iron losses [20][21][22]: where f and B are, respectively, the electrical frequency and flux density amplitude in ferromagnetic steel laminations, P 0 (W/kg) represents the mass power density at a reference operating point ( f = f 0 and B = B 0 ), and the iron losses (P iron ) are then given by multiplication of the mass power density (Equation ( 17)) of iron losses by the mass of the iron in each part of the machine where the magnetic flux density is evaluated: P f teeth is the mass power density of the teeth, and P f yokes is the mass power density of yokes.M iron yokes and M iron teeth are the mass of yokes and teeth, respectively.To limit the value of the losses in the iron when the heat is generated, it is necessary to limit the frequency of the currents delivered by the inverter, which fixes a maximum limit of pole pairs.

Copper Losses
The Joule losses depend on the volume of the copper used in windings and the current density.They can be calculated by the following equation: where V copp tot is the total volume of the copper in the stator, and J is the RMS current density.

Iron Part Thermal Circuit Modeling
The modeling of solid materials in the machine is based on the T-nodal method applied to a portion of cylindrical volume [23].The principle of this method is described in several studies [23][24][25][26].A basic thermal circuit of a portion of cylindrical material is shown in Figure 5. R 1 , R 2 are thermal radial resistances of an element of a cylinder, and R 3 is a negative thermal resistance that makes the temperature in contact point T c more realistic.The calculation formulas of these resistances and temperatures are given in Appendix A.

Slot and Winding Modeling
The thermal modeling of the stator winding is more complicated because of the involvement of fluid movement and conductive heat [27].Indeed, this winding, made of copper and insulators of different conductivities, is the highest temperature source in the electrical machine.It can also be noted that the copper conductivity is strongly dependent on the temperature.The control of this temperature allows the estimation of the losses by the Joule effect in the machine to foresee an appropriate cooling technology.The determination of the nodal model of this heterogeneous medium first requires its homogenization.It consists in assimilating this medium to a solid with equivalent conductivity λ eq .Several formulas and assumptions have been used in the literature [13,23,24] to calculate the equivalent conductivity of the medium.In this article, the used model is the one proposed by [23] and whose expression is: λ c , λ is and λ bob are, respectively, the thermal conductivities of copper, of the insulator and the equivalent conductivity of the homogeneous medium.Figure 6 shows the variations of λ bob as a function of filling factor k r .For a filling factor of 50%, the conductivity of the equivalent homogenous medium is about 0.7 W/m•K.

Heat Transfer Modeling
The heat is evacuated by convection in the air gap and in the external surfaces of the machine.The convection heat transfer in these interfaces is expressed as follows: h c , S and T A are the coefficients of convection, the area of the surfaces where there is convective heat transfer and ambient temperature.The determination of the convection coefficient, as well as of the thermal resistances, is based on the following hypothesis: • The lateral surfaces are in contact with a gas (natural air or hydrogen).

•
The internal surfaces are in contact with the liquid pushed by the blades of the rotor.
The convection coefficient can be generally calculated by the equation: where l is the length of the surface of the heat exchange and N u is the Nusselt number, which can be determined as follows: • For the external surface, the Nusselt number is given by: R e is the Reynolds number, and P r represents the Prendt factor.• For the internal side (gap), the determination of N uex is more complex.In the literature, there are several formulas that are used in the case of a thin space.In this paper, the expression used is the one described in [19].h c is decomposed into an axial component h ax and tangential component h tan .The Nusselt number for the axial convection, N uax , is given by Equation ( 23), and for the tangential convection, the Nusselt number N u tan is expressed in Equation ( 24) T a represents the Taylor number, which is given as a function of rotor speed, the geometry of the machine and the characteristics of the fluid that crosses the airgap.
The calculated coefficients of convection are 17.35 W/m 2 • • C and 5.26 × 10 4 W/m 2 • • C for external and internal (airgap side) surfaces, respectively.

Strategy
An optimization process was used to determine the geometry with the minimum amount of material volume possible, which leads to an average electromagnetic torque equal to the one imposed by the load (inducer).The optimal design must also respect thermal and structural constraints.The methodology used here is similar to that one described in [28].This process, described in Figure 7, starts from parameters initialized by the user (Table 3).By testing different configurations of the sizing variables, the optimization process selects the most accurate configuration to meet the specifications while respecting the constraints imposed.The fmincon function, which is a gradient-based nonlinear constraint, is used.Fmincon is a locally constrained nonlinear programming Solver from the Matlab optimization Toolbox; it uses a mathematical solution package that computes a quasi-Newton approximation of the Hessian.A detailed description of this method is given in [29,30].This method, which allows taking into account nonlinear constraints, is suitable for electrical machine design optimization, as shown in [30].

Design Variables
The design variables are the current density J, the electric load A L , the pole pair number, the magnet height h m and the axial length.The range of variation of these parameters is given in Table 4.The optimal values of design variables and Equations ( 1)-( 24) allow us to determine the flux density, the electromagnetic torque, the whole geometry of the motor and then the temperature in its active parts.

Objective
The objective function of the design optimization is to minimize the global volume of the motor's active parts, which may reduce its global cost and facilitate integration.The objective function, also called the fitness function, is then given by Equation ( 25): where V copper , V iron and V magnets are the copper, the iron, and the magnet volumes.

Constraints
Design constraints are integrated into the optimization process to ensure the feasibility of the system.For this study, six constraints are considered.The first, presented in Equation (26), requires that the average electromagnetic torque, T average , should be greater than the resisting torque T m corresponding to the turbine-rated torque.
The second, shown in Equation ( 27), concerns motor efficiency.It is fixed as greater than or equal to η min = 96%.
The third constraint, as expressed in Equation ( 28), avoids the oversizing of the converter.Therefore, the power factor must be higher than a minimal value fixed to 0.9: FP(Input var ) ≥ FP min = 0.9 (28) The fourth one, highlighted in Equation ( 29), emphasizes that the temperature in the stator windings must not exceed the maximum temperature value of the insulation T lim .In the case of the machine under investigation, a Class F winding is used at a maximum temperature of 155 ºC.
T max (Input var ) ≤ T lim (29) The fifth constraint, shown in Equation (30), concerns the teeth and slot widths, which must be considered to obtain rigid teeth and avoid too thin slots.The ratio of tooth width to the slot pitch is included between a maximum value R max and a minimum value R min .
The sixth constraint, presented in Equation (31), is introduced to avoid demagnetization of the magnets.For this purpose, it is essential to ensure that the reverse magnetic flux density applied to the magnet remains above the minimum coercive field.
The last constraint, given in Equation (32), limits the input voltage, which allows the determination of the number of conductors per slot N cd .

Optimization Results
The optimal values are obtained for a minimal volume of the active parts of the machine.Table 5 shows the optimal parameters and the volumes of the active parts of the optimized motors.These parameters were used to calculate the other geometric parameters and the average electromagnetic torque, as presented in Table 6.

Finite Elements Method
To validate the quality of the design model presented in the previous section, the Finite Element Method (FEM) is used.The electromagnetic behavior of the RD machine in 2D is evaluated using FEMM4.2software.The main objective of this part is to validate the results obtained by the proposed optimization algorithm.Therefore, the obtained results are for the rim-driven machine with Surface-Mounted Magnets (SMM) and the one using Halbach arrays with Progressive Magnetized Segments (HAPMS).Table 5 resumes the input values used for the progressive magnetization and the surface-mounted machines.These values are the output of the optimization process for both structures.
The design parameters obtained from the optimization algorithm are used as input parameters for simulating the motor with a 2D Finite Element Method ( FEMM4.2 software), as shown in Figure 8.In addition to the parameters given in Table 6, Table 7 gives the characteristics of the material used for the simulation.

Flux Density and Back Electromotive Force
The distribution of the flux in the active parts of the machine is shown in Figure 9.We suppose that the median magnet is facing toward the tooth.For this situation, the maximum values of the flux are 1.1 T in teeth and 1.68 T in the stator yoke, respectively.These values are approximatively matching with the ones obtained from the optimization algorithm.
Figure 10 compares the back EMF waveforms for the HAPM and MMS structures.The HAPM is quasi-sinusoidal, which leads to a stable electromagnetic torque.
Table 8 represents the comparison between the back EMF RMS values obtained from the optimization algorithm and by the Finite Elements Method (FEM) for a HAPM structure.Then, the two RMS values are about the same, which validates the results of the algorithm.

Electromagnetic Torque and Cogging Torque
The variation of electromagnetic torque, deduced as an interaction of the generated back EMF with sinusoidal currents, is given in Figure 11.Its average value is 11.1 N.m.The ripples are limited to less than 1%.It leads to a significant improvement in the results found in [13].
The cogging torque is dependent on the geometric parameters of the stator, the shape and thickness of the magnets, and their arrangement in the machine.Indeed, it is proportional to the flux density in the airgap.Halbach arrays with progressive magnetization in the rotor of the motor were designed to minimize the amplitude of this torque.Figure 11 represents the motor torques obtained by FEM for MMS and HAPMS structures.The comparison between the subfigures in Figure 11a,b shows that the cogging torque is reduced thanks to the sinusoidal waveform of the flux distribution.

Conclusions
This paper deals with the design and optimization of the active parts of a PMSM motor.The motor is used to rim-drive the inducer of a contra-rotating pump.The use of Halbach array permanent magnets in the rotor led to a quasi-sinusoidal waveform of back EMF.In addition, the HAPMS allowed for a reduction in the cogging torque.The simulation of the motor on the Femm 4.2 software confirms the advantage of the HAPMS structure compared to the MMS structure.The EMF generated by the machine using progressive magnetization is quasi-sinusoidal.Hence, the electromagnetic torque given in Figure 11a is relatively constant, with a ripple of less than 1%.This result is very satisfactory compared to more than 20% obtained for a surface-mounted magnet using a radial magnetization structure.
To validate the operation of the rim-driven inducer, a dynamic model of the counterrotation pump will be analyzed and simulated.As the model pump considers the thermal, electromagnetic, and hydraulic behavior of the pump, it is proposed to use the bond graph modeling as in [31][32][33][34][35][36].

Figure 3 .
Figure 3. Geometric parameters of the rim driven motor.

Figure 4 .
Figure 4. Peak magnetic flux density in the airgap as a function of l g and h m .at a characterstic points A, B and C. Magnetic flux density at points A, B and C is given in Table2.

Figure 5 .
Figure 5. Basic model of a cylinder portion with the T nodal method.

Figure 6 .
Figure 6.Variation of the equivalent conductivity of the homogenate medium as a function of the filling factor k r .

Figure 8 .
Figure 8. Structure with Halbach array magnet with progressive magnetization (a) and its 2D mesh (b).

Figure 9 .Figure 10 .
Figure 9. Distribution of flux density in the active parts of the motor (a), and Density Plot (b).

Figure 11 .
Figure 11.Torque, electromagnetic toque and cogging torque for HAPMS (a) and MMS (b) as a function of electrical angle.

Table 2 .
Variation of Bemax as a function of l g and h m .

Table 4 .
Range of variation of optimization variables.

Table 5 .
Optimal parameters and optimized volume.

Table 6 .
Inputs to the Femm code for MMS and HAPMS.

Table 8 .
Comparison of analytical values by the FEM of back EMF RMS.