Multidisciplinary Design Optimization of the Actuation System of a Hybrid Electric Aircraft Powertrain

Abstract

In the context of hybrid electric and full electric powertrains for future less-pollutant aircrafts, this paper focuses on the multidisciplinary design optimization (MDO) of the actuation system, including a surface-mounted PMSM in order to maximize the power density of the device: this study is a preliminary approach before integrating the whole powertrain. After an introduction of the MDO context, the analytical model of the electric motor is detailed. It integrates multi-physical aspects (electric, magnetic, mechanical, thermal, partial discharges and insulation, control and flight mission) and takes several heterogeneous design constraints into account. The optimization method involves a genetic algorithm allowing the reduction of the actuation weight with regard to a wide set of constraints. The results show the crucial sensitivity of the electro-thermal coupling, especially the importance of transient modes during flight sequences due to thermal capacitance effects. Another major point is related to the performance of the thermal cooling, which requires the introduction of an “internal cooling” in the stator slots in addition to the “base cooling” for stator and rotor. Gathering these analyses, the MDO leads to high power density actuators beyond 15 kW/kg with high-voltage–high-speed solutions, satisfying all design constraints (insulation, thermal, magnet demagnetization) over the flight mission.

1. Introduction

Power integration lowering both masses and volumes of powertrain devices embedded in transport applications is actually a great challenge for researchers, especially for actuation systems [1]. In ground transportation, numerous studies are focused on optimization strategies for power integration, such as in the review proposed in [2] for hybrid electric vehicles. If these challenges are huge for ground vehicles, reducing embedded weights in more electric aircrafts is a key driver for aeronautic evolution, as reviewed for example in [3]. It is especially true that typical “snowball effects” occur in aircrafts: the more embedded weight, the higher the wing surface and the more fuel burn. For example, [4] has shown that one additional ton would increase the fuel burn by 6% in the case of a regional aircraft. The ACARE (Advisory Councilor Aviation Research and Innovation in Europe) sets environmental objectives for 2050 technologies with a 75% reduction in CO₂ emissions per passenger kilometer and a 90% reduction in NOx emissions. The perceived noise emission of flying aircrafts should be reduced by 65%. These are relative to the capabilities of a typical new aircraft in 2000. More generally, the aviation industry actually faces the “revolution towards more electric aircrafts” [5,6,7,8].

Superconducting technologies in machines and power electronics may bring significant efficiency and weight reduction benefits over conventional components [9], but most of these promising technologies are currently, and for the near future, at a relatively low technology readiness level. Thus, conventional technologies are often preferred for transport applications. Several machine topologies can be selected and compared [10], especially in the automotive field for hybrid (and/or electric) applications [11,12]; in the automotive field, high specific power machines are embedded, until 4.3 kW/kg for the Tesla S60 rotating at 15,000 rpm [13,14]. Indeed, the trend towards high-speed actuation systems clearly exists in transport applications [15], which tends to reduce weight and volume. In aeronautic applications, the PMSM (permanent magnet synchronous machines) and the IM (induction machines) are the most adequate technologies. In addition to the higher efficiency of the PMSM over the IM, the PMSM also features higher specific torque, and this solution is seen today as the most suitable for weight optimization: this latter device has been selected in our study.

Optimizing the hybrid electric powertrain requires coupling a large set of domains, each involving heterogeneous phenomena and constraints in various physical fields. In that context, MDO (multidisciplinary design optimization) has become a methodological challenge itself, with several approaches and design strategies [4,16].

A typical MDO approach is proposed for optimizing the whole powertrain of future hybrid electric aircrafts in [17]. This study, partly presented in the proposed paper, is one part of a European Project in the framework of the Cleansky II EU project called “HASTECS” for “Hybrid Aircraft: Academic Research on Thermal and Electric Components and Systems.” In HASTECS, several studies [14,18,19,20,21] analyze innovative technologies (power electronics and advanced cooling, electric motors and its cooling, high voltage and partial discharges) and concepts for regional hybrid aircraft optimization in the case of a series hybrid architecture beyond 1 MW for power and beyond 1 kV for the bus voltage standard.

In his PHD Thesis, Pettes-Duler, M. [17] proposes the sizing optimization of the powertrain, integrating both the energy management strategy and the flight mission. This paper focuses on the actuation system including a PMSM, with the actuation weight being targeted as the optimization objective. As proposed in the HASTECS project, electric motors with high specific powers beyond 10 kW/kg are targeted in our optimization, with very high efficiencies (typically 97% at maximum power point). Aggressive targets have been chosen, but certain targets are already achieved in other studies [22]. In particular, Siemens [23], with the electric motor SP260D, has announced 5.2 kW/kg for a motor in flight tests. The University of Illinois [24] has designed a PMSM which would exceed a specific power of 13 kW/kg with an efficiency of 96%, showing that these targets may be reached. General Electric [25] has announced to reach a specific power beyond 10 kW/kg for a power inverter for Aircraft Hybrid-Electric Propulsion. Several tools can be used to design and optimize electric motors [26].

In Section 2, the context of the MDO process is introduced, synthesizing the optimization problem formulation.

The modeling task is one of the key issues for electric machine design: for example, [27] has recently presented an open-access electric machine design tool using MATLAB^® in order to enable an automated machine design. In our paper, a large set of heterogeneous sizing models are integrated, being strongly simplified in order to allow for solving this huge complex MDO process with acceptable computational times: these models are analytical for the electric motor [14] and its cooling [19], or based on parametric regression, as for partial discharge constraints [21].

One major contribution of this paper is related to the MDO process, especially coupling a large set of multidisciplinary constraints:

• Thermal constraints are involved in comparing steady-state and transient thermal behavior. Electro-thermal coupling is integrated into the optimization problem [28] emphasizing the “first-order influence” on the actuator performance. Several cooling systems can be assessed, as reviewed in [29].
• The actuation system being supplied by a high DC voltage bus, the Electrical Insulation System (EIS) in the PMSM [21], is also integrated through a simple surrogate model that involves insulation constraints due to partial discharge occurrence in stator slots.
• Regarding the bus voltage level and the equivalent impedance of the PMSM, the opportunity of a field-weakening strategy can also be assessed to optimize the performance according to the flight mission sequences. A supplementary constraint related to permanent magnet demagnetization is added for that purpose. This “flight mission-electric circuit-magnetic” coupling also affects the motor design and its specific power.

Several analysis and optimization results are presented in the last two sections: in Section 4, the huge influence of the transient behavior of the electro-thermal coupling is analyzed. The design choices are discussed, comparing steady-state and transient thermal models over the flight mission. Finally, in Section 5, the sensitivity of technological progress on actuation performance is analyzed, especially in terms of specific power and efficiency. This last part of the paper shows that the proposed MDO process allows reaching very high integration performance.

2. Context of the MDO Process

A complete design at the aircraft level is particularly complex because of the high number of decision variables with strong interactions between disciplines. A systemic study takes account of all (whenever possible) device couplings, far beyond summing local optimizations at the component level. It is within this framework that MDO is currently working because it allows gathering different fields around a single mathematical problem. Most often, this integrated design method is used to look at the sensitivity of the aircraft design, as well as its aerodynamic performance. For example, [27] linked the aerodynamic performance with non-linear physical phenomena occurring on the aircraft through an MDO. A robust and operational tool is presented in [30] in order to couple complex studies and highlight new aircraft concepts. Another optimal industrial trade-off for pylon design results from this: a demonstrator optimization test case has been implemented by the IRT Saint Exupery [31].

The study presented here is the preliminary step of a more complete MDO process is managed at the powertrain level in the thesis of Pettes-Duler [17]. This approach is applied to the design of a hybrid electric aircraft for regional flight. The final MDO process involves a large set of multi-physical aspects, as illustrated in Figure 1.

Before developing the system optimization at the powertrain level, a complete sensitivity analysis has been performed in [32], showing the major importance of the actuation system, especially the PMSM with regard to the powertrain weight and efficiency.

That is why a major preliminary step is to focus on the actuation part optimization, as presented in the next sections. Before presenting the models of the actuation system, Figure 2 illustrates the optimization formulation aiming at minimizing the electric motor mass.

3. Actuation System Modeling for MDO Process

In this section, we focus on the modeling of the inverter-fed surface-mounted permanent magnet synchronous machine (SM-PMSM) for high-speed applications. The theoretical backgrounds of the SM-PMSM model are based on the modeling of the open-circuit field and the armature reaction by current surface density [33]. From these theoretical backgrounds, an analytical sizing model based on loadability concepts has been derived [14,34] and has been revisited to be integrated in an MDO process. Such model seems to be the most suitable compared to finite element or permeance network models that lead to high computational times. Moreover, this model is dedicated to being further integrated in the whole powertrain optimal design, for which complexity and computation times constitute high challenges to be faced. In that context, the simplicity of analytical models and their short computational times are advantages that facilitate the optimizer’s tasks. In addition, less information is required for the system definition, unlike other models that sometimes require a large amount of input data. However, these models are less accurate compared with finite elements and they require simplification assumptions that are generally valid for a specific operating domain. Therefore, in order to optimize the complexity–validity trade-off, a validation of the analytical model has to be performed. In [33], the computation of the airgap magnetic flux density, the no-load back-emf and the torque obtained with this model have been validated by a finite element analysis.

The design model describes all the physical phenomena involved in the operation of the machine. Based on a power profile corresponding to a typical regional flight, and after defining the geometry of the machine, a magnetic model allows calculating most of its dimensions. Then, an electrical model specifies the different electrical circuit parameters of the machine which are required for coupling the motor with its power supply. At this stage, the structure of the machine is fully defined, and its mass and volume are determined. The losses in the different parts of the machine are computed for each operating point of the flight mission, enabling the electro-thermal coupling. Finally, the machine design must satisfy specific requirements in terms of thermal resistance of its various parts and the magnetic state of the magnets. Thus, a thermal model has been integrated taking into account the cooling systems as studied in [19]. This electrometrical MDO process is illustrated in Figure 3.

3.1. Geometrical Design Model

The electric machine geometry presented in Figure 4 corresponds to a SM-PMSM. It is the most promising topology in the aeronautic field, having the highest specific power among all machine types. Moreover, this topology can ensure flux density in air gap very close to sine wave, with very low saliency. Thus, with a sine wave-distributed winding, the analytical modeling is not too complex but is accurate enough for system integration [34]. From the set of sizing variables presented in

Table 1. Set of sizing variables.

Input Variables (in Blue is the Set of Decision Variables which Constitute the Degrees of Freedom for the Motor Optimization)
$T_{emot} \left(t\right)$	E-motor torque mission
${\Omega }_{emot}\left(t\right)$	E-motor rotational speed mission
${\mathit{R}}_{\mathit{a}\mathit{l}\mathit{e}\mathit{s}\mathit{a}\mathit{g}\mathit{e}}$	Bore radius of the electric motor
${\mathit{R}}_{\mathit{D}\mathit{r}\mathit{o}{\mathit{t}}_{\mathit{l}\mathit{m}}}$	Rotor diameter/rotor length ratio
${\mathit{R}}_{{\mathit{g}}_{\mathit{r}\mathit{a}\mathit{l}}}$	Air gap thickness/bore radius ratio
${\mathit{R}}_{\mathit{h}{\mathit{s}}_{\mathit{r}\mathit{a}\mathit{l}}}$	Slot height/bore radius ratio
${\mathit{R}}_{\mathit{p}{\mathit{m}}_{\mathit{r}\mathit{a}\mathit{l}}}$	Magnet thickness/bore radius ratio
${\tau }_{magnet}$	Pole pitch (= 100%)
${\tau }_{slot}$	Slot pitch (= 100% full-pitch winding)
$k_{carbon}$	Carbon fiber constant for sleeve equation
$\mathit{p}$	Number of pole pairs
$\mathit{q}$	Number of phases
$\mathit{n}\mathit{e}\mathit{p}\mathit{p}$	Number of slots per pole and per phase
${\mathit{N}}_{\mathit{c}\mathit{e}}$	Number of conductors per slot
$k_{fill}$	Fill factor in the slot
$J_a$	Permanent magnet flux density
${\mathit{B}}_{\mathit{y}\mathit{o}\mathit{k}\mathit{e}}$	Stator yoke flux density
${\mathit{B}}_{\mathit{t}\mathit{e}\mathit{e}\mathit{t}\mathit{h}}$	Stator teeth flux density
$B_{yok{e}_{rotor}}$	Rotor yoke flux density
${\mathit{V}}_{\mathit{u}\mathit{H}\mathit{V}\mathit{D}\mathit{C}}$	Ultra-high direct current bus voltage

, each dimension of the electric actuator can be calculated (see Figure 4): more information is available in [14,17].

This analytical model has been checked by a finite element software [33].

The computation of the first harmonic of the air gap flux density, $B_{ga{p}_{rms}}$, comes from geometrical decision variables. The electric angle, ${\gamma }_{elec}$, is half of the pole pitch angle, ${\tau }_{magnet}$. The root mean square value of the first harmonic flux density in the air gap is provided in

Table 2. First harmonic flux density computation process.

First Harmonic Flux Density Computation Process (in Blue is the Set of Decision Variables which Constitute the Degrees of Freedom for the Motor Optimization)
${\gamma }_{elec}=\pi \times \frac{{\tau }_{magnet}}{2}$	Electric pitch angle
$Bfon{d}_{ga{p}_{rms}}=\frac{R_{p{m}_{ral}}\times Ja}{log\left(-\left(1-R_{p{m}_{ral}}-R_{g_{ral}}\right)\right)}$	Air gap step value of the flux density
$B_{ga{p}_{rms}}=Bfon{d}_{ga{p}_{rms}}\times 2\sqrt{2}\times \frac{sin\left({\gamma }_{elec}\right)}{\pi }$	Fundamental value of the air gap flux density

.

All the other dimensions are determined by the decision variables (see

Table 3. Computation process of electric machine dimensions.

Dimensions of the SM-PMSM (in Blue is the Set of Decision Variables which Constitute the Degrees of Freedom for the Motor Optimization)
$L_{motor}=\frac{2\times R_{alesage}}{R_{Dro{t}_{Lm}}}$	Active length of the electric machine
$e{g}_{mag}=R_{alesage}\times R_{g_{ral}}$	Magnetic air gap thickness
$e_{magnet}=R_{alesage}\times R_{P{M}_{ral}}$	Magnet thickness
$h_{slot}=R_{h{s}_{ral}}\times R_{alesage}$	Slot height
$N_{enc}=2p\times q\times nepp$	Number of slots
$h_{yoke}=\frac{B_{ga{p}_{rms}}\sqrt{2}}{B_{yoke}}\times \frac{R_{alesage}}{p}$	Stator yoke height
$h_{yok{e}_{rotor}}=\frac{B_{yoke}}{B_{yok{e}_{rotor}}}\times h_{yoke}$	Rotor yoke height
$Rati{o}_{teeth}=\frac{2}{\pi }\times \frac{B_{ga{p}_{rms}}\sqrt{2}}{B_{yoke}}$	Tooth ratio (value at inner stator radius)
$e_{sleeve}=\frac{\left(R_{alesage}-e{g}_{mag}\right)}{\left(\frac{1}{k_{carbon}\times \max \left({\mathsf{\Omega }}_{emot}\left(t\right)\right)} \right)-1}$	Sleeve thickness
$e_{gap}=e{g}_{mag}-e_{sleeve}$	Mechanical air gap thickness
${\tau }_{teeth+slot}=\frac{2\pi \times R_{alesage}}{N_{enc}}$	Slot + tooth arc
$w_{slot}=\left(1-Rati{o}_{teeth}\right)\times {\tau }_{teeth+slot}$	Slot arc

).

In our model, slots are rectangular, and the slot length can be found according to

Table 4. Calculation of the slot length.

Slot Dimensions (in Blue is the Set of Decision Variables which Constitute the Degrees of Freedom for the Motor Optimization)
${\theta }_{slot}=\frac{w_{slot}}{R_{alesage}}$	Slot angle
$L_{slot}=2\times R_{alesage}\times \sin \left(\frac{{\theta }_{slot}}{2}\right)$	Slot length
$S_{slot}=L_{slot}\times h_{slot}$	Slot area
$S_{CU}{}_{total}=S_{slot}\times k_{fill}$	Useful copper area

.

The centrifugal pressure and the peripheral speed will be used as mechanical constraints. Thanks to the actuator dimensions, the centrifugal pressure exerted on the carbon sleeve from the maximum mechanical rotational speed value of the electric actuator, ${\mathsf{\Omega }}_{mech{}_{max}}$, can be estimated (see

Table 5. Sleeve mechanical constraints.

Mechanical Constraint Computations
$R_{yok{e}_{roto{r}_{OUT}}}=R_{shaft}+h_{yoke}{}_{rotor}$	Rotor yoke external radius
$R_{ma{g}_{OUT}}=R_{yok{e}_{roto{r}_{OUT}}}+e_{magnet}$	Magnet external radius
$R_{sleev{e}_{OUT}}=R_{ma{g}_{OUT}}+e_{sleeve}$	Sleeve external radius
$$\begin{gathered} P_{centrifuga{l}_{max}}=\frac{1}{3}\frac{{\mathsf{\Omega }}_{mec{h}_{max}}^2}{e_{sleeve}}\left({\rho }_{carbon}.\left(R_{sleev{e}_{OUT}^{}}^3-R_{ma{g}_{OUT}}^3\right) \\ +{\rho }_{PM}.\left(R_{ma{g}_{OUT}^{}}^3-R_{yok{e}_{roto{r}_{OUT}}}^3\right)\right) \end{gathered}$$	Maximum centrifugal pressure
$V_{periphera{l}_{max}}=R_{sleev{e}_{OUT}}\times {\mathsf{\Omega }}_{mech}{}_{max}$	Maximum peripheral speed

).

The maximum centrifugal pressure, $P_{centrifugal{}_{max}}$, and maximum peripheral speed, $V_{peripheral{}_{max}}$, constraints are set in order to design the right carbon sleeve thickness.

3.2. Partial Discharges Model

Partial discharges in electrical machines represent an important issue for systems’ reliability, especially in the context of a more electric aircraft where the combination of fast switching devices and long cables between power electronics converters can cause non-negligible overvoltage and lead to premature failures [35,36]. A complete study related to the electric insulation issue and especially the partial discharges that may appear with high-voltage-fed machines was achieved in Philippe Collin’s thesis [21]. In order to take this issue into account in an integrated MDO framework, the process is as follows: From the geometrical dimensions of the actuator, typically the number of conductors used in the windings, the partial discharges model sizes the insulation thickness required to avoid the corona phenomena appearance. A calibration abacus is used to take environmental conditions into account (temperature and pressure). In order to be conservative, the limit values have been used (limit temperature of the actuator and cruise altitude for the pressure). Then, a first set of parametric regression is used for the turn-to-turn insulation, increasing the thickness of the copper insulation, while the second parametric regression is used for the yoke-turn insulation, adding an insulation layer between the yoke and the conductors. The final part of this model checks if the windings have conveniently been integrated with all constraints. The final process is detailed in Figure 5.

Thanks to this model, the penalty in partial discharges due to the rise of temperature in the windings and pressure conditions is taken into account in this study. The electric design is then related to electromagnetic, mechanical and thermal aspects, but also to partial discharge constraints.

3.3. Performance Model based on Behn–Eschenburg Diagram

The Behn–Eschenburg linear model has been used in this study to characterize the actuator in phasor reference frame (see Figure 6). The winding configuration derives from the geometrical dimensions and allows the computation of the electric parameters (see

Table 6. Winding configuration computation process.

Electric Motor Winding Configuration (in is Blue the Set of Decision Variables which Constitute the Degrees of Freedom for the Motor Optimization)
$k_{tb}=\frac{\left(L_{motor}+{\tau }_{slot}\times \frac{\pi }{p}\times R_{alesage}+\pi \times {\tau }_{teeth+slot}\right)}{L_{motor}}$	Head winding coefficient
$kd1=\frac{\sin \left(\frac{1}{2}\times nepp\times p\times \frac{{\tau }_{teeth+slot}}{R_{alesage}}\right)}{nepp\times \sin \left(\frac{1}{2}\times p\times \frac{{\tau }_{teeth+slot}}{R_{alesage}}\right)}$	Twist factor
$kr1=\sin \left({\tau }_{slot}\times \frac{\pi }{2}\right)$	Shortening factor
$ki1=1$	Distribution factor
$k_{windings}=kd1\times kr1\times ki1$	Global winding factor
$k_{lc}=1$	AC coefficient losses

).

Three parameters are used to control the electromechanical actuator, the DC resistor, $R_{s_{DC}}$, the cyclic inductance, $L_{cS}$, and the electromotive force, $E_{rms}$.

Once again, the dimensions of the actuator are required to calculate the electric circuit parameters, as detailed in

Table 7. Electric circuit parameter computation.

Electric Circuit Parameters (in Blue is the Set of Decision Variables which Constitute the Degrees of Freedom for the Motor Optimization)
${\phi }_{rms}{}_{noload}=\frac{N_{ce}\times N_{enc}}{2.q.p}\times 2.k_{windings}\times L_{motor}\times R_{alesage}\times \sqrt{2}.B_{ga{p}_{rms}}$	RMS value of no-load flux
$L_p=\frac{4}{\pi }\times {\mu }_0\times k_{windings}^2{\left(\frac{N_{ce}\times N_{enc}}{2.q.p}\right)}^2\times \frac{L_{motor}\times R_{alesage}}{e_{magnet}+e{g}_{mag}}$	Self-inductance
$M=-\frac{1}{2}{L}_p$	Mutual inductance
$L_{cs}=L_p-M$	Cyclic inductance
$R{s}_{DC}=\frac{{\rho }_{CU}\times 2\times \frac{N_{ce}\times N_{enc}}{2.q}\times k_{tb}\times L_{motor}\times k_{lc}}{S_{C{U}_{total}}}$	DC resistance

.

The control strategy can be implemented from these parameters. A maximum torque per ampere strategy (“$I_d$ = 0”) has been performed with the capability of field-weakening over the mission. In particular, increasing speed leads to enhancing the motor voltage, which is limited by the bus voltage and the PWM source inverter: the amplitude of the motor voltage in the Park’s reference frame is limited ($Vdq\le Vd{q}_{max}$, defined in

Table 8. Initialization of the actuator control strategy.

Control Strategy
$I{q}_{mission}=\frac{T_{emot}}{p\times {\phi }_{rm{s}_{noload}}\times \sqrt{3}}$	Required torque => Iq current
$I{d}_{centre}=\frac{-{\phi }_{rm{s}_{noload}}\times \sqrt{3}}{L_{cs}}$	Center of the actuator circle
$Vd{q}_{max}=\frac{V_{uHVDC}}{2}\times m_a\times \sqrt{\frac{3}{2}}$	Stop-voltage (max inverter voltage)
$R_{mission}=\frac{Vd{q}_{max}}{L_{cs}\times {\mathsf{\Omega }}_{ele{c}_{mission}}}$	Radius of the actuator circle

). Maintaining a certain level of torque at higher speeds than the base speed induces a negative current “$I_d<0$”: this “counter-field current” reduces the air gap flux “${\mathsf{\Phi }}_g$”. This field-weakening operation operates at maximum voltage for high speeds.

From a geometrical point of view, the radius of the actuator circle is defined by the maximum available voltage, $V_{dqmax}$, and the actuator rotational speed over the flight mission, ${\mathsf{\Omega }}_{ele{c}_{mission}}$ (see Table 8). The intersection between the actuator circle and the current vector defines the operation point (OP_1,2).

When the voltage is not limited (without field-weakening), the circle must contain the operating point (blue circle in Figure 7b) and the equivalent condition is the following:

$$R_{mission}^2\ge I{q}_{mission}^2+I{d}_{centre}^2$$

(1)

The inequation becomes an equation when the voltage of the electric actuator reaches the maximum available voltage, $V_{dqmax}$. In this case, the field-weakening strategy occurs. The operating is then defined by the green circle. During that overspeed operation, the actuator usually operates at constant power, as illustrated in Figure 7a.

The blue circle shrinks to the green circle (increasing speed) and the current is shifted in phase to reach the operating point (see Figure 7b). Both operating points ($O{P}_1 \mathrm{and} O{P}_2$) are represented in the torque-speed plan (see Figure 7a). The circle characteristic is defined by:

$$R_{mission}^2=I{q}_{mission}^2+{\left(I_{d_{mission}}^{}-I{d}_{centre}^{}\right)}^2$$

(2)

This second order is derived in Figure 8 to set the current $I_{d_{mission}}$ in the case of the field-weakening operation:

$$I_{d_{mission}}^2-2\times \frac{-{\phi }_{rm{s}_{noload}}\times \sqrt{3}}{L_{cs}}\times I_{d_{mission}}+I{q}_{mission}^2+{\left(\frac{-{\phi }_{rm{s}_{noload}}\times \sqrt{3}}{L_{cs}}\right)}^2-{\left(\frac{Vd{q}_{max}}{L_{cs}\times {\mathsf{\Omega }}_{ele{c}_{mission}}}\right)}^2=0$$

(3)

Figure 8 shows the resolution process of the control strategy. The discriminant is computed to check if there are solutions. Finally, if solutions exist, the least restrictive I_d value is kept.

After this calculation process, the electromechanical actuator characteristics can be computed. The currents and voltages are derived from the d, q axis values (see

Table 9. Electric actuator characteristics.

$E_{rms{}_{mission}}={\phi }_{rms{}_{noload}}\sqrt{3} {\mathsf{\Omega }}_{ele{c}_{mission}}$	RMS value of no-load voltage
$V{d}_{mission}=R{s}_{DC}\times I{d}_{mission}-L_{cs}\times {\mathsf{\Omega }}_{ele{c}_{mission}}\times I{q}_{mission}$	d axis mission value of the voltage
$V{q}_{mission}=R{s}_{DC}\times I{q}_{mission}+L_{cs}\times {\mathsf{\Omega }}_{ele{c}_{mission}}\times I{d}_{mission}+$ $E_{rms{}_{mission}}$	q axis mission value of the voltage
$I{s}_{rm{s}_{mission}}=\frac{1}{\sqrt{3}}\sqrt{I{d}_{mission}^2+I{q}_{mission}^2}$	RMS one-phase current value
$V{s}_{rm{s}_{mission}}=\frac{1}{\sqrt{3}}\sqrt{V{d}_{mission}^2+V{q}_{mission}^2}$	RMS one-phase voltage value
${\mathsf{\Phi }}_{g_{noload}}={\phi }_{rm{s}_{noload}}\times \sqrt{3}$	No-load air gap flux
${\mathsf{\Phi }}_{g_{load}}=\sqrt{{\left({\mathsf{\Phi }}_{Rd}\right)}^2+{\left(I{q}_{mission}\times L_{cs}\right)}^2}$ ${\mathsf{\Phi }}_{Rd}={\mathsf{\Phi }}_{g_{noload}}+I{d}_{mission}\times L_{cs}$	Load air gap flux
$B_{yok{e}_{mission}}=B_{yok{e}_{noload}}\times \frac{{\mathsf{\Phi }}_{g_{field-weakening}}}{{\mathsf{\Phi }}_{g_{noload}}}$	Yoke flux density during the mission
$B_{teet{h}_{mission}}=B_{teet{h}_{noload}}\times \frac{{\mathsf{\Phi }}_{g_{field-weakening}}}{{\mathsf{\Phi }}_{g_{noload}}}$	Tooth flux density during the mission

).

Once currents and electric actuator characteristics are determined, losses can be computed from the mission profile according to

Table 10. Electric actuator losses computation process.

Electric Actuator Losses
$P_{JDC}{}_{}=3\times R{s}_{DC}\times I{s}_{rm{s}_{}}^2$	DC Joule losses in stator windings
$P_{Iron}=2\times \left(\sum f^{\alpha }{B}_{yoke}^{\beta }{M}_{stator{}_{yoke}}+\sum f^{\alpha }{B}_{teeth}^{\beta }{M}_{stator{}_{teeth}}\right)$	Iron losses in the stator yoke/teeth
$P_R{}_{}=2.C_{froul}\times \mathsf{\Omega }{ }_{ele{c}_{}}$	Friction losses in the bearings
$P_{Agap}{}_{}=k_1\times C_{f_{air}}.\pi .{\rho }_{air}\times {\mathsf{\Omega }}^3{ }_{ele{c}_{}}\times R_{alesage}^4\times k_{tb}\times L_{motor}$	Windage losses in the air gap
$P_{Arotor}{}_{}=C_{f_r}.\pi .{\rho }_{air}\times {\mathsf{\Omega }}^3{ }_{ele{c}_{}}\times R_{alesage}^5$	Windage losses in the two rotor surfaces

.

The losses are the inputs of the thermal model used for determining the temperature constraints in each part of the e-machine. Note that the electro-thermal coupling in mono-directional as temperature variations do not influence motor parameters in our model.

3.4. Thermal Model Using a Nodal Network

In the thermal modeling, two levels of a direct cooling device deeply studied in [19] have been implemented in the optimization process:

• “Base cooling” (first level), outside of the stator through a water jacket and inside of the rotor through a shaft cooling system.
• “Internal cooling” integrated in stator slots (“winding channel”) of the electric motor: this second level involves a more efficient cooling.

Both systems are water cooling circuits and are linked to a heat exchanger allowing to dissipate the heat flux in the plane environment (see Figure 9 and Figure 10).

The thermal model of the electric actuator and its cooling system is based on a lumped parameter thermal model: a nodal network is implemented in order to compute specific temperatures inside the actuator with a conductance matrix linking the different nodes. The capacitive transient effects have been integrated to filter the temperature transients inside the electric motor and to provide the temperature evolution during the flight mission. To do so, a thermal capacity matrix, G, of the motor system has been considered, with each capacity corresponding to the thermal capacitances of a given part of the system. Variations in input data such as the external temperature and the heat flow profile (motor losses profiles) are also considered. Thanks to this transient effect, the thermal limits can be released to achieve more compact and optimal sizing.

The thermal balance equation is as follows: $C.\frac{dT}{dt}=GT+\mathsf{\Psi }$, where C and G are respectively the matrices representing capacity and conductance effects, and T and Ψ are the vectors representing temperatures and heat generation in the thermal system (see Table 11).

The model is based on several simplifying assumptions: the thermal model is homogeneous in the axisymmetric rotor section, the radiation heat transfer is neglected and thermo-physical properties do not depend on the temperatures.

An implicit Euler method is used to solve these equations.

Table 11. Flux and conductance expressions for each transfer mode.

Flux and Conductance Expressions
Heat Transfer Mode	Conductance	Heat Flow Expression
Conduction (axial and ortho-radial)	$G_{axial}^{cond}\left(i,j\right)=\lambda .\frac{S_{ij}}{L_{ij}}$	$G_{axial}^{cond}\left(i,j\right)\left(T_j-T_i\right)$
Conduction (radial)	$G_{radial}^{cond}\left(i,j\right)=\frac{2\pi \lambda h}{\ln \left(\frac{r_j}{r_i}\right)}$	$G_{radial}^{cond}\left(i,j\right)\left(T_j-T_i\right)$
Convection	$G_S^{conv}\left(i\right)=H{S}_i$	$G_S^{conv}\left(i\right)\left(T_f-T_i\right)$
Fluidic flow	$G_{}^{fluid}\left(i,j\right)=\dot{m}{C}_p$	$G_{}^{fluid}\left(i,j\right)\left(T_j-T_i\right)$

Table 11. Flux and conductance expressions for each transfer mode.

Flux and Conductance Expressions
Heat Transfer Mode	Conductance	Heat Flow Expression
Conduction (axial and ortho-radial)	$G_{axial}^{cond}\left(i,j\right)=\lambda .\frac{S_{ij}}{L_{ij}}$	$G_{axial}^{cond}\left(i,j\right)\left(T_j-T_i\right)$
Conduction (radial)	$G_{radial}^{cond}\left(i,j\right)=\frac{2\pi \lambda h}{\ln \left(\frac{r_j}{r_i}\right)}$	$G_{radial}^{cond}\left(i,j\right)\left(T_j-T_i\right)$
Convection	$G_S^{conv}\left(i\right)=H{S}_i$	$G_S^{conv}\left(i\right)\left(T_f-T_i\right)$
Fluidic flow	$G_{}^{fluid}\left(i,j\right)=\dot{m}{C}_p$	$G_{}^{fluid}\left(i,j\right)\left(T_j-T_i\right)$

where:

• $S_{ij}$ is the surface of heat transfer between volumes represented by nodes $i$ and $j$,
• $S_i$ is the surface exposed to convection heat transfer,
• $T_f$ is the average temperature of the fluid surrounding surface $S_i$,
• $r_i,r_j$ are the radius of nodes $i$ and $j$ ($\mathrm{with} r_j\ge r_i$),
• $L_{ij}$ is the distance between nodes $i$ and $j$,
• $h$ is the height of the cylindrical object,
• $H$ is the heat transfer coefficient,
• $\dot{m}$ is the mass flow rate,
• $\lambda$ is the thermal conductivity,
• $C_p$ is the thermal capacitance.

Each temperature results from the thermal balance equation. Three of them have been chosen as constraints: the end-winding, the stator yoke and the magnet temperatures. Thanks to these constraints, the electric motor optimization can be performed with the mass of all components (active part + cooling system) as the objective function.

4. The Cooling System: A Leading Challenge

In aeronautics, “mass is the enemy number one”. Cooling plays a key role in the sizing of electromechanical components and constitutes a leading optimization constraint.

In Section 4.1, the difference between two thermal models has been firstly analyzed: a steady-state thermal model (R nodal network) and a transient model involving thermal capacitance effects (R,C nodal network). This comparative study highlights the importance of the coupling between the flight mission and the thermal modes which are simulated with the transient model (R,C) but not in the steady-state one. In Section 4.1, a “base cooling” system is considered involving two subsystems: a water jacket for stator external cooling and a liquid-cooled shaft system for the rotor.

In Section 4.2, the performance of the actuator optimization with this “base cooling” will be compared with a second level of thermal subsystem adding a stator slot “internal cooling” with the previous “base cooling”. The “internal cooling” device is directly integrated inside stator slots to be close to the heat dissipated by copper losses in windings.

4.1. Optimization Process

The clearing procedure [37] was used for optimizing the PMSM mass with regard to the design constraints. Clearing is a niching elitist genetic algorithm which usually outperforms standard genetic algorithms on difficult problems with multiple non-linear constraints and multimodal features [38]. All constraints were scaled and integrated into the objective function with penalty coefficients. The population size and the number of generations were respectively set to 100 and 200. Classical values for crossover and mutation rates were used (i.e., p_c = 1 and p_m = 1%). For each optimization case, multiple runs were carried out in order to take the stochastic nature of the algorithm into account and to ensure the reproducibility of results.

The global set of decision variables is provided in

Table 12. Set of decision variables.

Decision Variables	Name of the Variable	Lower Bound	Upper Bound
$V_{uHVDC}\left(\mathrm{V}\right)$	Ultra-high direct current voltage	540	2040
$R_{alesage}\left(\mathrm{m}\right)$	Inner radius of the stator	0.05	0.25
$R_{drot{}_{lm}}\left(\%\right)$	Ratio between rotor diameter and active length	50	125
$R_{\mathrm{h}s{}_{ral}}\left(\%\right)$	Ratio between stator slot and inner radius	10	150
$R_{g{}_{ral}}\left(\%\right)$	Ratio between the air gap thickness and the inner radius of the stator	1	10
$R_{pm{}_{ral}}\left(\%\right)$	Ratio between the magnet thickness and the inner radius of the stator	5	50
$B_{yoke{}_{max}}\left(\mathrm{T}\right)$	Maximum yoke flux density	1	1.53
$B_{teet\mathrm{h}{}_{max}}\left(\mathrm{T}\right)$	Maximum teeth flux density	1	1.53
$N_{ce}$ (-)	Number of conductors per slot	1	4
$n_{epp}$ (-)	Number of slots per poles and per phases	1	3
$p$ (-)	Number of pole pairs	1	7

.

Most of the decision variables are geometrical parameters. From this set, the geometrical dimensions of the electric motor are defined. The electric circuit parameters are derived to compute the actuator losses. The profile losses are used in the cooling model to estimate the temperatures inside the actuator. Finally, the partial discharges model checks the integration of the windings into the slot to avoid the phenomena.

The constraints of the optimization problem are listed in

Table 13. List of constraints used in the optimization process.

List of Constraints
$R_{shaft}\ge R_{shaft}{}_{min}$
$e_g\ge e_{g_{min}}$
$V_{peripheral}\le V_{periphera{l}_{max}}$
$P_{centrifugal}\le P_{centrifugal{}_{max}}$
$Missionperformed?\left(Flag\right)$
$demagnetizingcurrent\left(Flag\right)$
$T_{stato{r}_{yoke}}\le T_{stato{r}_{yoke{}_{max}}}$
$T_{copper}\le T_{coppe{r}_{max}}$
$T_{magnet}\le T_{magnet{}_{max}}$
$k_{fill}ok?\left(Flag\right)$
$D_{copper}\ge D_{copper{}_{min}}$
$Windingsintegrated?\left(Flag\right)$

. The first six constraints are checked after the electric motor model and the next three, which are related to the cooling model, are verified by simulating the flight mission. Finally, the last three constraints are calculated from the partial discharges model. When one constraint is not fulfilled, a penalized value (function of the number of constraints checked or not) is returned to the optimizer to facilitate the optimization convergence in a continuous way. Once all constraints are satisfied, the mass of the electric machine including its cooling is returned (being the objective function) to the optimizer. An optimal solution is found after the launch of 200 independent runs. This optimization process highlights the interaction between the sizing of an electromechanical component and its cooling.

4.2. Thermal Modeling: The Importance of Transient Modes Coupled with the Flight Mission

The reference flight mission profile integrated into the optimization is depicted in Figure 11. In this figure, shaft power and corresponding rotational speed are represented in per unit versus time during all flight phases (Taxi, Take off, Climb, Cruise, Approach and Landing, Taxi). The electric motor sizing and optimization are then performed by integrating the flight cycle in the loop according the MDO process described in Figure 2. All losses over the flight mission were computed according to the control strategy presented in Section 3.3.

Three different optimizations are carried out in the following sections: two are related to the e-motor with the base cooling (2025 target) presented in Section 3.4, using a steady-state or transient thermal model. The third optimization is applied on the e-motor with internal cooling (2035 target) using a transient thermal model.

A first level of thermal modeling has been used to optimize the electric motor mass: “the steady-state model” (in blue, left, Figure 12). A second level of “transient thermal model” involving transient modes is compared (in red, right, Figure 12). The improvements between both models are spectacular: the specific power has been multiplied by three (2.6 vs. 7.5 kW/kg), meaning that transient phases are highly sensitive regarding the motor sizing. In fact, the yokes based on magnetic sheets involve significant thermal capacitances which filter the temperature variations. As a result, the trade-off between specific power and energy efficiency is influenced by the thermal behavior: increasing the specific power of the electric actuator by considering thermal transient modes tends to decrease its efficiency. In that case, with the e-motor being less sensitive to losses, the latter are increased by lowering the mass.

A limit temperature of 220 °C has been used for the end-windings and the stator yoke in the optimization based on the “steady-state thermal model”; with this model, thermal limits have been released knowing that thermal capacitance effects would filter the thermal transients. For the “transient thermal model”-based optimization, a lower value (180 °C) has been taken to be conservative and safer. In the same way, two limit values have been used for the magnets: 200 °C for the “steady-state model”-based optimization and 150 °C for the transient.

Regarding the steady-state thermal model-based optimization, the maximum of the temperature values is reached during the take-off because of the high-power demand (see Figure 13). This steady-state temperature profile is really penalizing: electric and magnetic loads have to be reduced to satisfy thermal constraints during take-off, reducing losses (better efficiency) but increasing the actuator mass.

On the contrary, simulating transient modes involves filtering effects during the flight mission: the maximum values of temperatures are now switched at the top of the climb, because of the transient phase, which allows delaying the temperature rise (see Figure 14).

Furthermore, the transient behavior during a flight mission is also influenced by the control strategy. In order to be able to make a performance comparison, the torque-speed characteristics of both actuators are displayed in Figure 15 (in blue for the e-motor obtained with a steady-state (SS) thermal model, in red for that associated with the transient state (TS) thermal model and in grey for the situation of the “mission points” for each flight sequence). The actuator characteristic is limited by a maximum current (maximum torque), while its speed limit is set by a limit voltage (Vdq_max) related to the power electronic supply (related to the bus voltage, V_uHVDC, see Table 8). When this limit voltage is reached, the control strategy is as follows:

• Reduce the Iq current of the actuator, allowing it to go into a higher-speed phase without shifting the current (Id = 0).
• Add negative current in the d-axis of the actuator (Id < 0) to reduce the flux density (field-weakening) in the air gap, thus operating at constant power by increasing speed.

Both optimized machines (with both thermal models) use field-weakening capabilities. The field-weakening strategy optimizes the loss balance, reducing iron losses when the Joule losses are at their maximum (see Figure 16). Field weakening is exploited in both optimizations (SS and TS) at exactly the same time during take-off and almost during the whole duration of the climb.

The loss profile highlights the difference between both e-motor optimal designs. Figure 17 completes the analysis: the voltage reaches a maximum value (field-weakening) during the take-off/climb phase, while the current is increased. Regarding the absolute values, the current and the voltage have similar values for both optimizations, but the main difference resides in the current density, which is multiplied by three for the TS optimization.

With lower thermal constraints (as in the case of the transient thermal model), the number of pole pairs (see Figure 18) can be increased (also increasing iron losses) to reduce the actuator size. The field-weakening strategy allows counterbalancing the iron loss increase by lowering them during takeoff and climb sequences when the power demand is maximum (see Figure 16). Finally, thanks to this “thermal-control-mission” coupling, the specific power has been multiplied by three between both optimizations.

4.3. Sensitivity of the Cooling Performance

In order to once again increase the specific power of the electromechanical actuator, a second cooling solution has been implemented, adding an “internal cooling” (directly into the stator slots) to the “base cooling” (see Figure 9). As mentioned previously, the objective function is the electric motor mass, including the cooling system. The same set of decision variables has been used with the same list of constraints: only the cooling system has been changed. Let us recall that the “base cooling” is composed of a stator water jacket and a rotor shaft cooling.

Adding the “internal cooling”, the Joule losses are directly dissipated inside slots by the cooling channel. It also involves differences in terms of thermal dynamics: losses of the electric actuator are important with this high-performance cooling regarding the volume of the machine, and the transient phase is faster than before (see Figure 19). Once again, the stator yoke and the end-windings are the limiting temperatures and the magnet temperatures have been reduced because of the efficiency of the new cooling system. Enhancing the performance of the cooling device allows to significantly increase the specific power of the electric motor by a ratio higher than two (see Figure 20).

Regarding the decision variables, the inner stator radius has been reduced and compensated by a higher current density (related to the ratio between the slot height and the inner stator radius) thanks to the internal cooling system (see Figure 21). Bigger magnets ($R_{P{M}_{ral}}$) are obtained to increase the air gap flux density. The DC bus voltage is indirectly lowered with the internal cooling system optimization because of the integration of the cooling channel into the slots. Indeed, the directly cooled slots are so small that any space gain is important for the optimizer, leading to the reduction of the insulation thickness and of the insulation layer used for the copper; with this second optimization, the filling factor is close to 70%. The same trend is observed for the number of slots per pole and per phase ($N_{epp}$): the slots area cannot be too small while the internal cooling channel is integrated, which lowers $N_{epp}$ values. Increasing the number of pole pairs considerably reduces the yoke thickness and other sizes of the actuator. However, greater iron losses are expected, which is managed by the high-performance cooling. The same explanation can be used for the flux densities in the teeth: the cooling channels are in touch with the windings and the stator teeth, so the flux density in the teeth can be maximized. The yoke flux density has been decreased to minimize the iron losses in the yoke.

5. Conclusions

In this paper, a MDO was applied to a PMSM actuation system dedicated to a regional aircraft powertrain. This process involves a large set of multi-field couplings between electrical, geometrical, mechanical and thermal aspects, integrating the control strategy over the flight mission and also taking numerous optimization constraints into account. Among the environment constraints, partial discharges induce insulation issues for high-voltage electric actuators. The proposed study constitutes a preliminary step before the optimization of the whole powertrain embedded inside future hybrid electric aircrafts. With “mass being the enemy number one” in aeronautics, and while the electric motor is one of the most sensitive elements among all powertrain devices, we have focused this study on the mass reduction of the electric motor and its cooling. An analytical modeling of the electromechanical and thermal parts allows optimizing the actuator system. This model is a good trade-off between accuracy (at the system design level) and efficiency in terms of computational time to consider a future optimal integrated design at the whole powertrain level.

The main results of this MDO approach are related to the influence of thermal aspects, which constitute the leading constraint to enhance specific powers for electro-mechanical actuation. First, the influence of thermal transients during the flight mission has been highlighted, showing that main (copper and iron) losses can be balanced over the flight mission thanks to an optimal control involving field-weakening capabilities: this is a typical example of triple systemic coupling between the thermal model, the control strategy and the flight mission cycle. Second, it has been demonstrated that enhancing the performance of the thermal system by adding an internal cooling device directly inside stator slots strongly increases the power density of the electric motor. With the most performant cooling and with exploiting transient control of thermal modes over the flight mission, the specific power of the actuator system with its cooling may be beyond 15 kW/kg, which is a promising result for the field of hybrid electric propulsion in aeronautics.