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

: Hybrid electric aero-propulsion requires high power-density electric motors. The use of a constrained optimization method with the ﬁnite element analysis (FEA) is the best way to design these motors and to ﬁnd 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 ﬁnite 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 ﬁnd the best machine conﬁgurations and establish design rules based on the speciﬁcations 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 ﬂux weakening is clearly shown.


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 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, [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 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 mag- netic 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.

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. 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.

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. 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.

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

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.

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 and Figure 4. The chosen input parameters are general so that it is possible to construct machines with different topologies. Table 1. Input parameters of the parametric construction and meshing of an electric motor.

Description
Symbol Units Value in Figure 4 Internal 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. Figures 4 and 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: where  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. Figures 4 and 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: For rounded bottom slots, θe2 respects Equation (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 For rounded bottom slots, θ e2 respects Equation (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].
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].

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. The motor circuit is studied in receptor convention. 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 I s A Electrical angle from current "I s " to no-load emf "E" Ψ rad 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: 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 Q = I s cos(Ψ) (11) 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. The motor circuit is studied in receptor convention. 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 current waveforms are sinusoidal with an RMS value of Is 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: where ωs is the electric frequency in rad/s. RMS current values in the d axis and q axis can be obtained from Is and Ψ with the following relations.

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.

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.

Retaining Sleeve Stress Analytical Model
Input parameters required to calculate the mean stress on the retaining sleeve are indicated in the Table 3. 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]: 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.

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: 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]: where where Re is the Reynolds number T a is the thickness of the air gap and µ 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.

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: R s = N 2 s N coil L av σ 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.

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 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.

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.

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

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 ( ), 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. 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. Figure 10 and Table 4 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).
+ k e dB(t) dt 3 2 k f dt (22) wires protrude by 1 cm on each side of the stator and that coil ends are 1.41 times longer than the coil span.  Figure 10 and Table 4 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

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 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. It should be noted that four different mass values were used according to the problem studied.

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 (≈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.
positions per simulation, resulting in a time of 100 s ( 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.

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 Figures 12 and 13, that for a given number of poles, total losses remain constant at the same value as the electric frequency increases.

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 Figures 12 and 13, that for a given number of poles, total losses remain constant at the same value as the electric frequency increases. 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.

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 Figures 12 and 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 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).

Cooling Effort Approximation with Equivalent Current Density Product
Once optimizations were finished, the cooling efforts were approximated with the equivalent current density product "AJeq" [12]. This parameter approximates the cooling 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 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).

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 12 Wm −3 Ω −1 with a direct water cooling as explained in [12]. This parameter can be calculated as follows: 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 12 Wm −3 Ω −1 for 30 kW/kg, 3.6 × 10 12 Wm −3 Ω −1 for 20 kW/kg, 2.8 × 10 12 Wm −3 Ω −1 for 10 kW/kg and 1.5 × 10 12 Wm −3 Ω −1 for 5 kW/kg. However, lower equivalent current density products can be obtained by setting it as the function to minimize by the optimization.
where Ptot are the total losses, ktb 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 12 Wm −3 Ω −1 for 30 kW/kg, 3.6 10 12 Wm −3 Ω −1 for 20 kW/kg, 2.8 10 12 Wm −3 Ω −1 for 10 kW/kg and 1.5 10 12 Wm −3 Ω −1 for 5 kW/kg. However, lower equivalent current density products can be obtained by setting it as the function to minimize by the optimization.

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

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.