Magnetic FEA Direct Optimization of High-Power Density, Halbach Array Permanent Magnet Electric Motors

Abstract

Hybrid electric aero-propulsion requires high power-density electric motors. The use of a constrained optimization method with the finite element analysis (FEA) is the best way to design these motors and to find the best solutions which maximize the power density. This makes it possible to take into account all the details of the geometry as well as the non-linear characteristics of magnetic materials, the conductive material and the current control strategy. Simulations were performed with a time stepping magnetodynamic solver while taking account the rotor movement and the stator winding was connected by an external electrical circuit. This study describes the magnetic FEA direct optimization approach for the design of Halbach array permanent magnet synchronous motors (PMSMs) and its advantages. An acceptable compromise between precision and computation time to estimate the electromagnetic torque, iron losses and eddy current losses was found. The finite element simulation was paired with analytical models to compute stress on the retaining sleeve, aerodynamic losses, and copper losses. This type of design procedure can be used to find the best machine configurations and establish design rules based on the specifications and materials selected. As an example, optimization results of PM motors minimizing total losses for a 150-kW application are presented for given speeds in the 2000 rpm to 50,000 rpm range. We compare different numbers of poles and power density between 5 kW/kg and 30 kW/kg. The choice of the number of poles is discussed in the function of the motor nominal speed and targeted power density as well as the compromise between iron losses and copper losses. In addition, the interest of having the current-control strategy as an optimization variable to generate a small amount of flux weakening is clearly shown.

1. Introduction

Today, reducing fuel emissions is a global concern that poses many challenges. In the civil aviation industry, manufacturers have pledged to halve aircraft carbon emissions by 2050 compared with 2005. Electric motors are becoming an essential alternative to the propulsion system because of their potential to reduce aircraft emissions. However, their use requires batteries, converters and high power density electrical machines [1,2]. Permanent magnet synchronous motors (PMSMs) are widely used in these applications due to their high efficiency, high power density, and small size, making it a key energy conversion element and an efficient alternative to conventional motors [3]. In addition, when these motors have an array of permanent magnets with continuous changes in the direction of magnetization, further advantages are added, such as magnetic flux concentration effect and a more sinusoidal radial magnetic field in the airgap. It also provides higher torque when a large amount of magnet is used. This array of magnets is known as the Halbach array and the machines with this design are called surface mounted permanent magnet motors with a Halbach array [4]. They are widely used in aircraft applications, as shown in the high power density motor presented by [5] and the motor developed by the National Aeronautics and Space Administration (NASA) [6].

Using an iterative optimization process is the best way to design a machine that meets requirements and specifications, such as power density, volume, torque and weight. The design method consists in finding the best feasible solution according to the application and determines its topology, dimensions, characteristics of magnetic materials, conductive materials, and current-control strategy [7]. During this iterative step, the optimization algorithm uses the information on the performance estimation of each motor configuration to propose a new set of motor dimensions that improves the performance [1]. For that, it is necessary to use accurate numerical models able to estimate the motor’s performances according to their size and their control method [8].

Coupled multiphysics modelling methods must be used to properly assess motor performance [9,10]. Indeed, there are many interactions between the physical phenomena in a machine and they must be considered as explained in [11]. Figure 1 shows the multiphysical couplings and the connection between them. However, the trade-off between the accuracy of these models and computational time is a critical factor when using an iterative optimization procedure that may require a significant number of iterations before converging to a solution [12].

In [9], the authors present an electric motor optimization process for aircraft propulsion. They develop an optimization tool for applications of PMSMs with high power density, considering multiphysics couplings such as electromagnetic, thermal, and mechanical models, with the ability to explore wide search spaces and compare different types of electrical motors. Four radial PMSMs were optimized for a speed range of 8 to 20 krpm considering different pole numbers corresponding to an upper electrical frequency limit of 1.5 kHz. The power of the machine was 1 MW and the optimization target raised the power density of the motor. In the same order of ideas for aircraft propulsion, ref. [10] proposes an optimization procedure for the design of permanent surface magnet motors using a cooling system based on the airflow produced by the speed of the aircraft. To solve this problem, the authors integrate magnetic, electric, thermal and mechanical models in an optimization procedure using a differential evolution algorithm. Two optimization problems were studied: single objective problem minimizing losses on a 60-kW motor with a nominal speed of 3000 rpm and multi-objective problems, maximizing the power density with the restrictions imposed. This shows the feasibility of a direct optimization of an electric machine by including a finite element model, but the resolution of the magnetic field calculation was performed in magnetostatic. Therefore, some additional losses like eddy currents in solid conductive parts were neglected such as in the magnets or the retaining sleeve. An FEA with an external circuit coupling and a magnetodynamic resolution makes it possible to evaluate more precisely the effects of a current control with flux weakening and to take them into account in the optimization process.

Knowing the target specifications, the first step to develop an optimization process is to identify the key parameters that needs to be computed. Examples of specifications for a SMPM motor are the nominal speed, the nominal torque, the maximal heating, and the efficiency [13]. Key parameters that must be evaluated are the electromagnetic torque, iron losses, copper losses, eddy current losses in the magnets, magnetic flux densities in the teeth and in the stator yoke, the torque ripple, mechanical constraints on the retaining sleeve, windage losses, peak temperature of the winding and peak temperature of the magnets [14]. Indeed, depending on the class of their insulation, winding must not operate over a certain temperature [15]. Likewise, magnets have a Curie temperature, and their performances decrease with temperature [16]. Magnetic, electrical, mechanical, and thermal models are then required [11].

The aim of this study is to present a direct magnetodynamic time stepping FEA optimization process using the MATLAB “sqp” algorithm through the constrained nonlinear optimization function “fmincon” and discuss its computational cost and implementation. This method has the advantages that iron losses, magnet losses, magnetic flux densities and electromagnetic torque estimations are more accurate. Moreover, magnetic non-linearity is considered with FEA and an original approach is to set the sinusoidal current control strategy (module and angle) as an optimization variable. The implementation of this method makes it possible to find general design rules such as the choice of the number of poles according to the nominal speed to achieve a target power density.

2. Implantation and Modelling Overview

Optimization processes are composed of an optimization method and modelling techniques so the algorithm can compute the key parameters (or constraints) at every iteration. As explained in [1], there are two main families of constrained nonlinear optimization methods: methods with gradient-based algorithms and methods with intelligent optimization algorithms. The first family includes the conjugate gradient algorithm and the sequential quadratic programming algorithm while the second family includes evolutionary algorithms as genetic algorithms and multi-objective optimization algorithms. These methods and their advantages are widely discussed in [17]. Figure 2 shows the families of optimization methods and their modelling techniques.

From Figure 2 can be observed that there are several iterative optimization methods that use FEA and/or analytical models, among these are analytical modelling, direct FEA modelling, space mapping [18] and metamodeling [19] for constraints and objective function evaluations.

Analytical modelling consists, as its name suggests, of an analytical model. It has a low numerical cost, and it is generally easy to implement but less accurate than other modelling techniques. The author of [20] shows a direct analytical model using an equivalent magnetic network where magnetic non linearities are considered paired with a mathematical model analysis (MMA) optimization algorithm. The authors concluded that there is consistency between both methods with a maximum error of 5% for torque results and up to 14.5% for results in the winding temperature model.

Direct FEA modelling consists of automated parametric FEA simulations that are driven by an iterative optimization algorithm. These can be in 2D or in 3D and can be dynamic or static. These models are computationally expensive but accurate. For instance, [21] shows a two variable direct 3D FEA optimization process for a PMSM with a novel rotor structure and [22] proposes a modified robust design optimization (RDO) paired with a parametric FEA model for brushless DC motors. The model direct FEA employs design analysis models to evaluate the objectives and constraints in the model, and by using an algorithm (for example GA), optimizes all parameters at the same time [1]. It is characterized by its simplicity in implementation and can be integrated within different software, for example, MATLAB and EXCEL. The finite element field calculation software such as ANSYS MAXWELL or CEDRAT FLUX can be driven by MATLAB or EXCEL to carry out multiphysics optimization. Finite elements analysis (FEA) is widely used in the design analysis and calculation in electrical machines due to its good accuracy but requires a large computation time and is harder to implement. For example, consider a PMSM with 10 optimization parameters by using genetic algorithms, with an approximate number of 200 iterations to converge and an overall population size equal to 50 (10 × 5), this implies that around 10,000 (200 × 50) optimization model evaluations are required, which translates into high computing times [1].

Space mapping combines a computationally cheaper model with a correction based on a more expensive model that helps in the optimization using a parameter space transformation [11,18]. The interest of a space mapping approach is also clearly shown in [23]. It is slower than direct analytical models because of FEA simulations but is more accurate. Metamodeling uses a population of FEA simulation results to generate a model using, for example, neural networks [24] or the Kriging method [25]. It is also possible to use multiple metamodeling methods and compare their performances [8]. This method is numerically expensive when it comes to the FEA simulation results generation, but it is otherwise accurate and fast.

The design of PMSMs by the FEA and piloted by an optimization method helps to realize the optimal conception of the motor by preselecting different topologies, evaluating for each of them the key parameters, the restrictions imposed, and the objective function posed. This approach helps to find the best solution and helps to establish general rules for the design of PMSMs.

As it can be seen in Figure 3, the “fmincon” function of the MATLAB software with the “sqp” algorithm option is the main optimization algorithm. Analytical models for the stress on the retaining sleeve, the copper losses and the windage losses were also implemented in MATLAB. Moreover, a direct FEA coupled electromagnetic model was implemented on FLUX2D software, a parametric FEA software. FEA simulations were called from MATLAB and ran by the FEA software via python files.

3. Parametric Motor Simulations with FEA Software

FEA electromagnetic simulations are performed with a parametric FEA software (Flux2D version 10). Python scripts were written so machines with different numbers of poles, number of slots and other geometric parameters can be constructed, meshed, and simulated by the FEA software but launched from the MATLAB software.

3.1. Parametric Motor Construction and Meshing

First, a python script constructs the machine in the FEA software by placing points, connecting them with lines, filling in the created faces with materials and meshing the domain. The main input parameters of the machine construction, and meshing are presented in

Table 1. Input parameters of the parametric construction and meshing of an electric motor.

Description	Symbol	Units	Value in Figure 4
Internal Dirichlet boundary condition diameter	Ddiri	m	0.02
Rotor internal diameter	Dint	m	0.044
Rotor back iron thickness	Tb-iron	m	0.007
Magnets thickness	Tmag	m	0.01
Thickness of the retaining sleeve	Tsleeve	m	0.002
Thickness of the air gap	Ta	m	0.002
Thickness of the stator teeth wedge	Tw	m	0.003
Stator slots depth	ds	m	0.012
Stator yoke thickness	Ts	m	0.01
Axial length of the active part	L	m	0.0935
Inner slot opening angle	θe	rad	0.15
Outer slot opening angle	θe2	rad	0.165
Wedge-opening angle	θw	rad	0.05
Number of winding layers in a slot	Nlayers	-	2
Number of turns in a coil	Nturns	-	1
Number of magnets in a rotor pole	NsecMag	-	6
Total number of slots of the machine	Nslots	-	24
Number of slots in the study domain	Nslots2	-	6
Total number of rotor poles pairs	P	-	2
Factor between the external diameter of the machine and the diameter of the external boundary condition	fdiri	-	1.1

and Figure 4. The chosen input parameters are general so that it is possible to construct machines with different topologies.

The python script also has other parameters including the choice between rounded or rectangular bottom slots, magnetization orientation of the magnets and the distribution of phase coils in every layer of every slot. Figure 4 and Figure 5 respectively show an electric motor with rounded bottom slots and with rectangular bottom slots. The magnetization orientations of the magnets for the Halbach rotor are calculated exactly as explained in [26]. The outer slot opening angle, θ_e2, can be calculated so the width of the stator teeth is constant. With rectangular bottom slots, θ_e2 can be calculated in the function of the motor geometry as follows:

$${\theta }_{e2}=\frac{2\pi }{N_{slots}}-\frac{D_w\left(\frac{2\pi }{N_{slots}}-{\theta }_e\right)}{D_{es}}$$

(1)

where

$$D_w=D_{int}+2T_{b-iron}+2T_{mag}+2T_{sleeve}+2T_a+2T_w$$

(2)

$$D_{es}=D_w+2d_s$$

(3)

For rounded bottom slots, θ_e2 respects Equation (4).

$${\theta }_{e2}=\frac{2\pi }{N_{slots}}-\frac{D_w}{D_{es}}·\left(\frac{2\pi }{N_{slots}}-{\theta }_e\right)·\left(1+\sin \left(\frac{{\theta }_{e2}}{2}\right)\right)$$

(4)

As it can be seen in Figure 4, boundary conditions must be specified. The internal Dirichlet condition is optional, but the external boundary condition is mandatory. To speed up the simulation, the study of the 2D geometrical domain was minimized to a portion the transversal plane. A cyclic or anti-cyclic boundary condition was added in the function of the number of poles in the studied domain. If the number of poles drawn is even, the periodicity must be cyclic. On the contrary, if the number of poles drawn is odd, the periodicity must be anti-cyclic.

The mesh was performed with an automatic mesh generator and the number of nodes (nn) was used to adjust the density of the mesh as shown in Figure 5.

The choice of the number of nodes makes it possible to find a compromise between the computation time and the precision of the magnetodynamic simulations. It was found that the magnetic losses and electromagnetic torque computed with a mesh having 4124 nodes (Figure 5b) are very close (less than 0.05%) to those obtained with a mesh with 8595 nodes (Figure 5a). The main difference is on the computation time which is divided by two. Therefore, for the iterative optimization process, it is very advantageous to reduce the number of nodes and to use the mesh of Figure 5b.

A finer mesh can also be imposed in the magnet to compute eddy current losses in the magnets. Magnets skin depths were meshed with three elements along the skin depth d_magSkin. An approximate skin depth was calculated with Equation (5). This computation assumes that the magnet material is linear, that magnets are a semi-infinite block and that there is no wave reflection at the interface between the magnets and the rotor back iron [27].

$$d_{magSkin}=\sqrt{\frac{{\rho }_{mag}}{\pi {\mu }_{mag} f_{mag} }}$$

(5)

where ρ_mag is the resistivity of the magnets in [Ω × m], μ_mag is their permeability in [H/m] and f_mag is the highest field harmonic frequency with a significant magnitude in [Hz]. As discussed in [28] and [29], this frequency often depends on the inverter PWM frequency. However, since the phase currents imposed here were perfectly sinusoidal, this frequency was supposed equal to the slot frequency of the stator, as shown in (6). N_rpm is the rotor speed in [rpm].

$$f_{mag}=\frac{N_{rpm}}{60}{N}_{slots}$$

(6)

3.2. Electric Circuit and Stator Winding Supply

Once the machine geometry is designed with the help of FEA software, a python script must specify the external circuit of the machine with sinusoidal current sources and impedances. Input parameters for the stator current waveforms are shown in

Table 2. Input parameters of equivalent electric circuit.

Parameter Description	Symbol	Units
Electrical angle between phase A axis and the pole axis for the initial rotor position	θac	rad
RMS value of the phase current	Is	A
Electrical angle from current “Is” to no-load emf “E”	Ψ	rad

. The motor circuit is studied in receptor convention.

The current waveforms are sinusoidal with an RMS value of I_s and are forming a three-phase balanced circuit. As can be seen in Figure 6, θ_ac is the electric angle between the rotor magnetic pole axis in its initial position and the magnetic axis of phase A. For example, the value of θ_ac for the machine drawn on Figure 6 is 0.2618 rad or 15°.

A phase shift must be added to the phase currents in function of θ_ac and Ψ. The three current waveforms can be expressed as:

$$i_a=\sqrt{2}{I}_s\sin \left({\omega }_{s}t-{\theta }_{ac}-\mathsf{\Psi }\right)$$

(7)

$$i_b=\sqrt{2}{I}_s\sin \left({\omega }_{s}t-{\theta }_{ac}-\mathsf{\Psi }+\frac{2\mathsf{\pi }}{3}\right)$$

(8)

$$i_c=\sqrt{2}{I}_s\sin \left({\omega }_{s}t-{\theta }_{ac}-\mathsf{\Psi }-\frac{2\mathsf{\pi }}{3}\right)$$

(9)

where ω_s is the electric frequency in rad/s. RMS current values in the d axis and q axis can be obtained from I_s and Ψ with the following relations.

$$I_D=I_s\sin \left(\mathsf{\Psi }\right)$$

(10)

$$I_Q=I_s\cos \left(\mathsf{\Psi }\right)$$

(11)

3.3. Full Load Magnetodynamic Simulations

Electromagnetic torque, iron losses and magnetic inductions in the teeth and the stator yoke were calculated by a 2D magnetodynamic simulation launched by a python script where the rotor was rotating at its nominal speed. The initial position of the rotor was set to zero and its final position was set to three polar steps for the computation the iron losses with the Bertotti model.

4. Analytical Models

This section presents the analytical models used in parallel with the FEA simulations to compute copper losses, aerodynamic losses, the mechanical stress on the retaining sleeve.

4.1. Retaining Sleeve Stress Analytical Model

Input parameters required to calculate the mean stress on the retaining sleeve are indicated in the

Table 3. Input parameters for retaining sleeve stress analytical model.

Parameter Description	Symbol	Units
Magnet material density	ρm	kg/m3
Angular width of a magnet	θm	rad
Angular width of the back iron below a magnet	θs	rad
Retaining sleeve density	ρs	kg/m3
Radius measured below the magnets	ra	m
Radius measured below the retaining sleeve	rla	m
Radius measured outside of the retaining sleeve	rs	m
Thickness of the sleeve	Tsleeve	m

.

Stress on retaining sleeves of high-speed SMPM machines has been discussed in many articles such as [30,31]. The mean stress on the retaining sleeve is [32]:

$${\sigma }_{\theta mean}={\mathsf{\Omega }}^2\left[\frac{\left({\rho }_m{\theta }_m\left({r_{la}}^3-{r_a}^3\right)+{\rho }_s{\theta }_s\left({r_s}^3-{r_{la}}^3\right)\right)}{\left(6T_{sleeve}sin\left(\frac{{\theta }_s}{2}\right)\right)}\right]$$

(12)

where Ω is the rotor’s angular speed in rad/s. The maximum stress on the retaining sleeve is estimated by performing mechanical FEA simulations before the optimization. When the maximum acceptable stress on the sleeve is known, the minimum sleeve thickness can be calculated with the MATLAB function “fmin”. The function to minimize is the difference between the desired stress and the calculated stress on the retaining sleeve with its thickness as the optimization variable. This method can be used to save a variable and a constraint in the main optimization problem. Indeed, the minimum thickness of the retaining sleeve can be calculated separately at every optimization iteration when the maximum mechanical stress is known.

4.2. Aerodynamic Losses Analytical Model

Windage power losses are determined exactly as in [33] and considering that the airflow is turbulent, that there are no shrouds and that rotor poles are not salient. These losses can then be expressed as:

$$W_{losses}=\pi C_d{\rho }_{air}{\left(\frac{D_s}{2}\right)}^4{\mathsf{\Omega }}^{3}L$$

(13)

where ρ_air is the density of air, D_s is the diameter measured at the exterior of the retaining sleeve, L is the active length of the rotor and C_d is the skin friction coefficient. The skin friction coefficient is the value of C_d that respects the following equation [32]:

$$\frac{1}{\sqrt{C_d}}=2.04+1.768 ln\left(Re\sqrt{C_d}\right)$$

(14)

where

$$Re=\left(\frac{D_s}{2}\right)\frac{T_a\mathsf{\Omega }}{\left(\frac{\mu }{{\rho }_{air}}\right)}$$

(15)

where $Re$ is the Reynolds number T_a is the thickness of the air gap and $\mu$ is air dynamic viscosity. This equation can be solved numerically with the MATLAB function “fmin” in the same way the minimum retaining sleeve thickness equation was solved.

4.3. Copper Losses and Phase Resistance

Phase total resistance can be evaluated as following when the average coil turn length is known. The number of phases is equal to 3:

$$N_{coil}=\frac{N_{slots}{N}_{layers}}{6}$$

(16)

$$S_{cu}=\frac{\alpha S_{slot}}{N_{layer}}$$

(17)

$$R_s=\frac{N_s^2{N}_{coil}{L}_{av}}{{\sigma }_{cu}{S}_{cu}}$$

(18)

where N_coil is the number of coils per phase, S_slot is the total area of a slot minus the area used by the wedges holding the coils or the tooth tips, σ_cu is the copper conductivity, α is the copper fill coefficient and other definitions can be found in Table 1. Once the phase resistance is known, copper losses can easily be evaluated.

$$P_{copper}=3R_s{I_s}^2$$

(19)

5. Direct Optimization Algorithm Implementation in MATLAB

The motor geometry was optimized with the MATLAB function “fmincon” with the “sqp” algorithm. The function handles provided in parameters of “fmincon” started a Flux2D server that launches the python scripts which construct and simulate the motor with the finite element software. Optimization variables were scaled from 0 to 1 and the output of the constraint function was also manually scaled. The overall implementation in MATLAB is presented on Figure 7.

In summary, “fmincon” takes in arguments a function handle for the computation of the total losses “functionToMinimise”, a function handle for the computation of the constraints “constraintsFunction” and an initial guess of the optimization variables. The function “functionToMinimise” to simulate motors in finite element software, computes all the useful results and saves them for further use in “constraintsFunction”. The function “constraintsFunction” obtains the last saved results, compares them to the chosen values and returns a vector that represents constraints violations.

6. Comparative Study of Optimized Motors

To illustrate the advantages of this design procedure, this method was applied to compare machines optimized for a power of 150 kW and for different rated speeds. The aim of the study is to determine the trade-offs between power density and efficiency. It is also possible to deduce from this study, several design rules concerning the choice of the number of poles, the current density and the flux weakening operation at the nominal speed.

To carry out this comparative study, it was necessary to repeat the optimizations of the machine by changing the number of poles, the nominal speed, and the power density. Only the three-phase Halbach rotor machines with an internal rotor were considered. The number of slots per pole (spp.) and per phase in the stator is always set equal to 2.

6.1. Selected Motor Structure

Since the machine structure is imposed with a number of slots per pole and per phase equal to two, the minimum study area for the finite element simulation is always one pole even when the number of poles changes. The configuration of the stator winding can also remain the same. A winding configuration with two layers and a short pitch of 5/6 was selected (Figure 8).

The average coil length ($L_{av})$, as shown in Figure 9, was estimated considering that wires protrude by 1 cm on each side of the stator and that coil ends are 1.41 times longer than the coil span.

$$CoilSpan=\frac{2\pi r_{cslot}{N}_{slotSteps}}{N_{slots}}$$

(20)

$$L_{av}=2\left(L+0.02\right)+2\left(1.41\ast CoilSpan\right)$$

(21)

Figure 10 and

Table 4. Material properties for the finite element simulations.

Description	Symbol	Value	Description	Value
Magnets remanent induction		1.19 T	Magnets density	8300 kg/m3
Magnets relative permeability		1	Rotor back iron density	7872 kg/m3
Magnets resistivity		1.6 × 10−6 Ωm	Stator iron density	7600 kg/m3
Stator sheets thickness	d	0.2 mm	Retaining sleeve density	1500 kg/m3
Stator sheets conductivity	σ	1,776,199 S/m	Retaining sleeve relative permeability	1
Hysteresis loss coefficient	kh	183.9 WsT−2m−3	Copper density	8933 kg/m3
Excess loss coefficient	ke	0.4 W(Ts−1)−3/2m−3	Copper resistivity	2.59 × 10−8 Ωm
Stator sheets stacking factor	kf	0.98	Maximum mechanical constraint on retaining sleeve	2650 MPa

detail material properties.

The magnets are NdFeB magnets and their demagnetization curve was supposed linear in the simulations. Rotor back iron was fabricated of 1010 steel (Figure 10a) and stator sheets are NO20 (Figure 10b).

Iron loss parameters presented in Table 4 are for NO20 laminations and these parameters are kept constant even if the number of poles and the electric frequency vary. Iron losses are computed element by element by the FEA software using the Bertotti method with Equation (22).

$$d{P}_{moy}=k_h B_m^2 f k_f+\frac{1}{T}{{\displaystyle \int }}_0^T\left[\frac{\sigma d^2}{12}{\left(\frac{dB}{dt}\left(t\right)\right)}^2+k_e{\left(\frac{dB\left(t\right)}{dt}\right)}^{\frac{3}{2}}\right]k_{f}dt$$

(22)

6.2. Optimization Problem Definition

The objective function of the optimization problem consists in minimizing the total losses by imposing the nominal power and the total mass as constraints. Thus, it is possible to compare machines having the same power density (5 kW/kg, 10 kW/kg, 20 kW/kg, and 30 kW/kg) even if the rotor speed is different.

The number of poles and the rotor speed must be fixed before each optimization. It was decided not to exceed an electric frequency of 2.5 kHz given the uncertainty associated with the estimation of the magnetic losses if the same parameters are always used for the Bertotti model.

Table 5. Optimization variables.

Parameter Description	Symbol	Units
Rotor internal diameter	Dint	m
Axial magnetic circuit length	L	m
Stator yoke thickness	Ts	m
Stator slots depth	ds	m
Magnets thickness	Tmag	m
Slot opening angle	θe	rad
RMS value of the phase current	Is	A
Electrical current angle between current and no-load emf	Ψ	rad

presents the eight variables used for this optimization problem. Other geometric parameters such as the thickness of the rotor yoke rotor Tb-iron were fixed and either calculated from these variables or from analytical models such as the thickness of the sleeve. The constraints are detailed in

Table 6. Optimization constraints.

Parameter Description	Value	Units
Minimum motor power	150	kW
Maximum motor mass	5, 7.5, 10, 15	kg
Maximum induction in stator teeth	1.8	T
Maximum induction in the stator yoke	1.5	T
Maximum mechanical stress on the retaining sleeve	2650	MPa

. It should be noted that four different mass values were used according to the problem studied.

6.3. Losses in Function of Nominal Speed and Number of Poles

The results correspond to an optimization with a mesh of 4124 nodes, with 24 rotor positions per simulation, resulting in a time of 100 s ($\approx$0.02778 h) to complete each simulation. A total of 300 simulations were considered, which corresponds to a total of 8.33 h, for a personal computer (16 Gb of RAM with an i7 processor)

Optimization results were gathered to compare the performance of machines with the same power density. Thus, Figure 11a shows machines optimized for 30 kW/kg; those of Figure 11b, machines of 20 kW/kg; those of Figure 11c, machines of 10 kW/kg; and those of Figure 11d, machines of 5 kW/kg.

The total losses of optimized motors were smaller as the nominal speed increase. This can be understood by the fact that these motors had the same nominal power but not the same torque since they were not rotating at the same speed. It can be seen that the differences between the loss minima (or the efficiency maxima) are relatively small as the power density varies. However, the position of these minima is different depending on the nominal speed of the motor and its number of poles. The minimum losses were 3.5 kW above 30 krpm in the case of machines at 30 kW/kg (efficiency of 97.7%), 3 kW above 25 krpm in the case of machines at 20 kW/kg (efficiency of 98%), 2.3 kW above 15 krpm in the case of machines at 10 kW/kg (efficiency of 98.4%), 3.2 kW above 8 krpm in the case of machines at 5 kW/kg (efficiency of 97.9%).

There were several assumptions taken into account in the study, such as the consideration of an ideal current form in the machine, which neglects the presence of a large number of harmonic currents injected by the inverter. The current ripple changes the magnetomotive force and generates additional losses in the stator and rotor core, due to hysteresis and eddy current losses [34,35,36]. These additional losses would be significant in a high speed motor application but several modifications of rotor geometries are possible to minimize them such as a magnet segmentation [36].The main source of error in these calculations comes from the analytical models used for the calculation of magnetic losses and windage losses [37]. However, since the parameters of these models were kept constant in this comparative analysis, the trends that can be deduced from them are valid.

6.4. Effect of the Current Control Angle at Nominal Speed and Compromise between Copper Losses and Iron Losses

The optimization method minimizes losses in the machine by acting on the electrical angle Ψ between phase current and no-load emf. If the angle Ψ is greater than 0, the stator magnetic reaction produces a flux weakening effect. This action makes it possible to adjust the stator iron losses in relation to the Joule losses and find the best compromise on the total losses. Indeed, it can be seen in Figure 12 and Figure 13, that for a given number of poles, total losses remain constant at the same value as the electric frequency increases.

According to Bertotti model, iron losses should be almost proportional to the square of the electric frequency, but the optimization process decreases Ψ for higher speeds which reduces flux densities and iron losses in the stator. This can clearly be observed on the 6 poles, 10 kw/kg motors. Indeed, total losses reached a plateau at approximately 20 krpm which also corresponds to the point where the angle Ψ was starting to decrease. It can also be observed that when total losses remain constant, iron losses were almost equal to copper losses. For a given power density, total losses are minimized almost in the range of electric frequency regardless of the motor’s nominal speeds, as shown in Figure 13.

Table 7. Optimal electrical frequency ranges in function of targeted power densities.

Targeted Power Density	Optimal Frequency Range
kW/kg	Hz
10	800–2000
20	>1400
30	>2500

summarizes the optimal electrical frequency ranges.

For a given nominal speed and a targeted power density, the optimal frequency range can help the designer to set the number of poles of the machine. For example, according to these results, a 10 kW/kg, 14 krpm motor should have between 4 and 10 poles (Figure 13).

6.5. Cooling Effort Approximation with Equivalent Current Density Product

Once optimizations were finished, the cooling efforts were approximated with the equivalent current density product “AJ_eq” [12]. This parameter approximates the cooling effort required to cool the motor and should not be greater than 2 × 10¹² Wm⁻³Ω⁻¹ with a direct water cooling as explained in [12]. This parameter can be calculated as follows:

$$A{J}_{eq}=\frac{P_{tot} {\sigma }_{cu}}{k_{tb} 2\pi R L}$$

(23)

where P_tot are the total losses, k_tb is the end-winding coefficient, R is the inner stator radius, L is the active length of the motor and σ_cu is copper conductivity. A fixed end-winding coefficient of 1.4 was used. Figure 14 presents the equivalent density products of the motors obtained by direct optimization. Total losses were still the function to minimize.

The minimum cooling efforts of optimized motors in the function of their power density were approximately 6 × 10¹² Wm⁻³Ω⁻¹ for 30 kW/kg, 3.6 × 10¹² Wm⁻³Ω⁻¹ for 20 kW/kg, 2.8 × 10¹² Wm⁻³Ω⁻¹ for 10 kW/kg and 1.5 × 10¹² Wm⁻³Ω⁻¹ for 5 kW/kg. However, lower equivalent current density products can be obtained by setting it as the function to minimize by the optimization.

7. Conclusions

The design of a permanent magnet motor using an optimization process based on FEA ensures excellent accuracy and very good sensitivity of the numerical model. The finite elements modelling of the magnetic field allows to consider nonlinear phenomena such as magnetic saturation, the rotor movement, the presence of eddy currents and the motor control method. This reduces the number of approximations and simplifies the assumptions of the numerical model. Although the computation time was increased, we showed that an eight-variable optimization problem can be solved in less than 9 h with a typical personal computer. Under these conditions, it is easy to repeat the optimizations to compare different machine topologies for a given specification and determine the best solution.

We used this design procedure to analyze the trade-offs between power density and efficiency and determine the most suitable number of poles and electrical frequency for a 150-kW application at different rotational speeds. More than 90 machines were optimized for four different power densities (30 kW/kg, 20 kW/kg, 10 kW/kg, 5 kW/kg) and then compared. This represents approximately 800 h of computing time with a single personal computer. We found several configurations of machines with similar performances knowing that the optimization adjusts the d-axis current to achieve a flux weakening and to make a better compromise between the joule losses and the magnetic losses.

For all power density values, this comparative analysis shows that it is advantageous to increase the electrical frequency to minimize the losses but a flux weakening control must be used to balance the magnetic losses and the copper losses. The electrical frequency can be adjusted by appropriately selecting the number of motor poles relative to the rated operating speed. In this case, the motor efficiency can be over 97.7% in a 150-kW application, regardless of motor power density.

This type of design procedure is very efficient, but it is essential to consider all the phenomena and the physical couplings. In particular, a reliable thermal model must be added to take into account the cooling method in the mass and power density estimation.