A Holistic Electrical Machine Design Tool for More-Electric and Hybrid-Electric Aircraft

The power density of electrical machines required for more-electric and hybrid-electric aircraft has become a critical aspect and target of optimization. This paper looks at the development aspects of a multi-domain tool to aid for system-level optimization of electrical machines within next-generation high power-density aerospace applications. The electromagnetic, thermal and mechanical aspects are wholly integrated, thus enabling the optimization including the non-active mass. The implementation and overall architecture are described with the aim of investigating the impact of coolant flowrate, coolant medium as well as efficiency on the achievable power to mass ratio.


I. INTRODUCTION
ith the globally increasingly stringent emissions legislations and fuel economy requirements, companies in the transportation sector are actively and intensely researching new technologies which often involve electrification and hence the use of electrical machines either for motoring or generation. The performance targets in this type of work are various and depend a lot on the specific industry and application. For example 'high power density' is often a key phrase to distinguish new developments. In the land transportation industry, more specifically for road transportation systems, where volume is often highly constrained, the key power density metric is the power to volume ratio or [kW/L] [1][2][3]. On the other hand, for the aerospace industry, mass minimization, rather than volume, is critical and the key power density metric is the power to mass ratio or [kW/kg] with various numbers published to show achievements of particular developments, such as a recent 5 [kW/kg] for a light electric aircraft [4].
System architects working on the concept generation and system integration of the aforementioned more electric transport architecture often face a bottleneck when it comes to the electrical machine. Whilst comprehensive libraries of say, high speed bearings are normally available either through supplier input or in-house designs, for the high-performance electrical machines targeted in such work the available data is very limited. Doing machine sizing in a manual manner for the range of options which the system architects want to investigate is too much time consuming and impractical, while narrowing down the options risks in missing the system optima altogether. This paper describes the development of such a tool for investigating options with the intent of maximizing the power density [kW/kg]. The targeted power density is determined as the developed power per total mass, including the mechanical parts required for the machine to operate. The main feature of the tool is that it considers electromagnetic, thermal and mechanical design aspects simultaneously. This allows optimizing electrical machine taking into account mechanical and thermal constraints which usually significantly affect machine design [5,6]. In the first part the methodology behind the tool development is described. The tool is then implemented and used to compare the achievable [kW/kg] for a 1MW 20krpm inner rotor PM machine configurations considering different types of coolant, and different flow-rates for various efficiency limits.

II. METHODOLOGY
In many classical optimization approaches, the optimization is first done on the electromagnetics, then a housing is designed around the optimized electromagnetic design. However, the housing can be a very significant proportion of the total mass and integrating the housing design with the overall machine kW/kg optimization has significant potential for extra power W density entitlement. Appropriate multidomain calculators, which serve as the essential building blocks with which the kW/kg optimization is performed, are required. To this end, electromagnetic, thermal, and mechanical analytical models are developed, as shown in Fig. 1.
The arbitrary SPM machines are defined in terms of their characterising geometries, constituent materials defined by their magnetic, mechanical and thermal properties, as well as the coolant properties which include the coolant temperature and flow rate. The following sections detail the multidomain calculators implemented and used within the optimization tool.
The methodology allows calculation of electromagnetic, thermal and mechanical parameters of synchronous permanent magnet machine with Halbach array. An overall structure of multidomain calculator is shown in Fig. 1. The electromagnetic calculator has a set of input parameters which describe the machine geometry such as stator inner and outer diameters, active length, magnet height, tooth width, tooth height, etc. Fig.  1. The machine active power, speed, material properties and number of turns per phase also should be specified as input parameters. The electromagnetic calculator also estimates copper and iron loss which are then used as input parameters for the thermal calculator Fig. 1. The thermal module does a calculation of temperature state of the machine winding, stator back-iron, teeth and housing Fig. 1. Coolant properties as well as flow rate and initial temperature have to be defined as input parameters. The mechanical module estimates active and passive machine mass and also calculates mechanical stress of the retention sleeve Fig. 1.  The arbitrary machine geometry is initially used to compute the no-load magnetic field according to the analytical model of machines with a Halbach array [7]. Under linear behaviour of magnetic materials, the solution for the fundamental of the radial component of flux density in the air gap for the IR topology is given in [8].
The machine is analyzed in generation mode. The d-q axis model is adopted with the rotor considered aligned to the q-axis. Imposing the d-axis current to zero, the resulting phasor diagram is shown in Fig. 3. The phase current I can then be estimated from (1): where ܲ is the electromagnetic power, m is the number of phases and ߟ is the efficiency. The reactance is calculated by means of analytical expressions [9], and the power factor is thus derived accordingly. For the loss calculation, the DC copper losses and the iron losses are considered. The flux densities for the iron loss calculation are evaluated by means of a linear magnetic circuit calculation and considered as average values on the overall structure. The armature current flux density is evaluated using the approach described in [10].
Given the specific loss of a lamination material say at a frequency of 60 Hz and at an induction of 1 T -W ி̴ ǡଵ , for any stator fundamental frequency -݂ ௦ , and iron flux-density level ‫ܤ‬ ௦ , the specific iron losses can be approximated from [11] (2): where ‫ܤ‬ ௦ is the on-load flux density.
The machine geometry can be described by a set of variables (3) -(6), Fig. 1: The variable ܹܶ (1) defines a tooth width ratio and varied in the range from 0 to 1 I. The variable ‫ܪܶ‬ (2) defines tooth height ratio and also varied in the range of 0 -1 I. The variable ݇ describes a ratio between an outer and inner stator diameters and takes values higher than 1 I. For the mechanical calculator, the input is the same arbitrary active electrical machine geometry of Fig. 2, together with the mechanical properties of the constituent materials. Based on these inputs, the mechanical calculator estimates the mass of the active materials (i.e. mass of copper, magnets and iron) and also sizes an appropriate housing around the active geometry with the intent of calculating the non-active mass. Fig. 5 details the housing sized around the Inner Rotor SPM configuration showing the axial cross-section (left) and radial cross-section (right) with the colour grey denoting the stator housing, colour green denoting the rotor shaft, colour yellow for the stator laminations, colour brown for the copper and colour red/blue for the permanent magnets. There are a total of eight inactive parts which are sized for a mechanical factor of safety of at least 1.5 and assuming the machine to be foot mounted. The aforesaid eight inactive parts are numbered in Fig. 5 and listed in Table I. The stator housing (grey) is made up of items#1-#3 which correspond to two cylindrical shells (item#1), enclosing end-flanges (item#2), and cooling channels. All the aforesaid components are made of lightweight aluminum having a density of 2700 ሾ ݇݃Ȁ݉ ଷ ] in order to minimize the mass.
The rotor assembly supporting the Halbach magnet array is also made up of three constituent parts, items#4-#6 corresponding to a rotor shaft made of high strength, magnetically permeable 17-4ph stainless steel which has a hollow cross section to minimize the mass (item#4), two rotor balance plates made of stainless steel (item#5), and a lightweight carbon-fiber sleeve (item#6) which retains the magnet under compression.

 
A minimum shaft diameter ‫ܦ‬ is calculated based on the torque transmission requirement as described by equation (7): where T is the shaft torque, and J is the hollow shaft's polar moment of inertia.
With the minimum shaft diameter determined, the bearing inner diameter can therefore be selected. For the drive-end a cylindrical bearing is selected (item#7), while for the non-driveend a pair of back-to-back angular-contact bearings are used (item#8) as shown in Fig. 5. In determining the bearing mass, linear correlations between bearing inner diameter and mass of the bearing are derived based on the available bearing data [13].

V. THERMAL CALCULATOR
For the thermal calculator, the inputs are the same arbitrary electrical machine geometries as well as the thermal conductivities of the materials used. In addition to these, the thermal calculation tool requires as inputs the coolant thermal properties such as the density, viscosity and conductivity, the inlet temperature and the flow-rates. Based on the aforesaid inputs, together with the calculated machine losses from the electromagnetic calculator, the thermal calculator estimates the temperatures at various locations within the electrical machine. The cooling technique strongly impacts the power density level [6] and for the purpose of this paper intensive, high flow-rate direct oil-cooling [14] is considered as described in the following sections since the framework of this research targets high power density.
In order to efficiently extract the heat out of the machine, it is best to locate the heat sink as closely as possible to the heat source. The minimization of the thermal resistances between the heat sources and the coolant enables efficient heat removal [15]- [16]. For the SPM machines under investigation, rectangular cooling channels are thus created in the stator slots, back iron and teeth as shown in Fig. 5, where most of the heat generated by the copper and iron losses is located.
Referring to Fig. 5, oil flows into the flooded stator chamber from the non-drive end of the machine, impinging on the endwinding surfaces through multiple jet-nozzles, as shown by the black arrows. The oil at the non-drive end chamber is forced to flow through the cooling channels in the stator core, and within slots, as shown by thearrows in Fig. 5. A thin sleeve is applied in the airgap to ensure that the rotor chamber is free of oil and this helps avoid high windage losses since the machine is running at high speed. A separate oil flow is also added to the machine housing-jacket, as shown by the yellow arrows in the aforesaid figure.
Based on the above cooling concept, a lumped parameter thermal model is developed. For convective heat transfer, the heat transfer coefficients inside the rectangular cooling channels depend on the flow patterns. The flow patterns are in turn determined by the evaluation of the Reynolds number ܴ݁, defined in [17].
When the Reynolds number ܴ݁ of the flow in cooling channel is less than 2300, the flow is said to be laminar, whilst when higher than 2300 the flow in the cooling channels is turbulent. For laminar flow, the heat transfer coefficient ℎ can be calculated from [16] For the end winding cooling, the heat transfer coefficient is estimated based on previous experimental work by the authors and measured to be around 1000 [ܹ /݉ ଶ ‫]ܭ‬ to 3000 [ܹ /݉ ଶ ‫]ܭ‬ depending on the flow rate and oil jet design [12].

VI. OPTIMIZATION MODEL
In order to determine the optimum design for the IR topology presented, a Genetic Algorithm (GA) which is embedded in the commercial optimization software modeFrontier TM is used. The optimization model, shown in Fig. 6, consists of a so-called 'scheduling project node' and a 'nested optimization procedure'. In essential terms, the scheduling project node initializes the pole-pairs and speed variables for the nested optimization procedure. For each combination of speed and pole-pair numbers defined in the 'scheduling project node', the aforesaid parameters are transferred as input parameters to the 'nested optimization procedure'. The machine geometry is characterized by seven defining variables, grouped under the heading 'Input Variables', shaded in blue in Fig. 6 and related with the geometry presented in Fig. 2. These variables are the split ratio ' ݇ ', stator inner diameter ' ܵ ூ ', aspect ratio ‫݉ݎ݂_݇'‬ ', magnet height ‫ܯ'‬ ௧ ', tooth width coefficient 'ܹܶ ', tooth height coefficient ‫ܪܶ'‬ ', and carbon fiber sleeve thickness ' ‫ܮܵ‬ ௧ '. These parameters are easily understandable by the electrical machine designer and allow to parametrize the main geometry of the machine.
The initial population is typically set to around 300. The optimization algorithm is selected to achieve fast computation and solution robustness. The final part in the solver options core is the Matlab Interface block enabling information exchange between the algorithm and the multi-domain calculators presented in the previous sections.
This interface enables the optimization to access the electromagnetic-thermal-mechanical model and return the results to the optimization algorithm. This brings the setting-up of the optimization problem to the third and final level, shaded in green in Fig. 6, in which the Output Variables and Optimization Targets are defined. The problem in hand is single-objective in nature, targeting the minimization of the machine total mass (active plus non-active parts). In achieving this target, a number of constraints are defined on the outputs. The first two constraints 'Tendsmax' and 'Tcoilmax' relate to the thermal limitations, and ensure that for any design to be considered feasible the temperature in the winding must not exceed a defined limit, which for the case in hand is set as 200 [℃] corresponding to class C insulation. Also related to the thermal domain, minimum practical cooling channel areas are defined by the parameters 'SlotChW' and 'ToothChW'. For a design to be considered feasible the power factor 'PF' and the efficiency 'eff' must also be higher than defined thresholds (in this case power factor over 0.75 and efficiency over 97%) while the on-load tooth and core flux densities are limited to up to 2.  Fig. 7. The machine is optimized for the case of 1 and 2 slots per pole per phase. The inlet oil temperature is 50°C with a flow rate of 150 l/min. It can be seen from Fig. 7 that the power density for 1 slot per pole per phase machine gradually increases from 11 kW/kg to almost 26 kW/kg for the investigated range of pole pairs. A specific optimum at 4-5 pole pairs can be observed for 2 slots per pole per phase design. This is due to the increase of slot number and taking into account that there are practical considerations which translate to minimal manufacturable channel dimensions. For the purpose of comparison and within the scope of system level studies the optimization of machine power density was also performed for the case of both oil as well as water cooling, taking into account the different electric conductivities of both fluids, for a range of flow rates from 10 l/min to 150 l/min. For the higher flow rates (50 l/min, 150 l/min) the machine power density is almost the same for both coolants. It gradually increases from 10 kW/kg for 2 pole pairs topology to around 25-26 kW/kg for higher pole number machines - Fig. 8, Fig. 9. For lower flow rates (10 l/min, 20 l/min) water cooling technique has an advantage and achieves up to 17 kW/kg and 23 kW/kg against the 11 kW/kg and 15 kW/kg for oil cooling - Fig. 8, Fig. 9. The results of 1 MW, 20000 rpm machine optimization with no temperature constraint are given in Fig. 10. For the goal of theoretical investigation the winding maximum temperature constraint of 200°C is disabled. The efficiency constraint is varied in the range of 0.95 -0.99. As expected the machine power density is very sensitive to the machine efficiency. The machines with the efficiency around 99 % have a power density rating from 8.5 kW/kg to 9.5 kW/kg while decrease of efficiency to 95% allows reaching more than 30 kW/kg. On the other hand the winding temperature would exceed reasonable limits for the machines with an efficiency rating less than 97 % especially for lower pole number designs Fig. 1. Thereby, for the considered cooling concept for megawatt class high speed machine the most reasonable efficiency range is from 97 % to 98.5% which depends on the number of pole pairs. The power density that can be achieved for the aforementioned efficiency range is from 15 kW/kg to 26 kW/kg Fig. 10, with the maximum coil temperature not exceeding 200°C for the considered efficiency range as shown Fig. 1.  This work described in details a design tool for analytical calculation of synchronous permanent magnet machine with Halbach array. This type of machine allows reaching very high power density and can be suitable for more electric and hybrid aircraft application. The tool includes electromagnetic, thermal and mechanical calculators and can be used for fast sizing and optimization. As an example, the optimization of 1 MW, 20000 rpm machine for different pole pairs and number of slots per pole per phase is given in the paper. The impact of coolant medium and efficiency constraints on machine power density and winding temperature are also presented.