1 Introduction

At present, electric propulsion for aircraft is one of the major issues in aeronautic research, due to its advantages in comparison to combustion propulsion. The biggest problem of purely battery-driven aircraft is their short range. In spite of big efforts in research the energy density of batteries is still low. For this reason battery-driven aircraft cannot travel long distances compared to combustion driven airplanes. Many different alternative concepts of propulsion architectures are examined (e.g., [1]). However, those hybrid architectures are complex and therefore not easily applicable for small airplanes. Apart from that, the range of a fixed configuration can be increased by an efficient flight guidance. Therefore, the characteristics of efficient flight operation have to be found, the parameters of optimal steady-state operation must be determined and delivered to the pilot via cockpit-display.

Sachs demonstrates special elementary flight and range performance characteristics of electrically propulsed airplanes with the help of basic mathematical correlations in [2, 3]. Other research shows estimation of range and correlation with design parameters [1], using the Breguet-Formula for electric airplanes. Traub [4] examines the impact of the so called Peukert–Effekt on the maximum range and specifies corresponding optimal flight speed, the used model is extended by Avanzini et al. [5]. Donateo et al. [6] compares the estimated gross endurance with net endurance, determined by detailed simulation. Ostler et al. [7] shows an estimation of flight performance and correlating flight speeds. Efficiency criteria especially derived for electrically propulsed airplanes are presented in [8, 9]. Falck et al. [10] accomplishes a trajectory optimization for battery propulsed airplane, considering the thermal constraints of the involved subsystems of the power train. In contrast, [11] uses a detailed energy consumption model to perform a trajectory optimization and a steady-state evaluation. Steady-state phases in the determined optimal trajectory can be correlated with energy-optimal criteria.

This study describes an appropriate mathematical model, which is able to perform the optimization. After that, efficiency criteria for range efficient climb and range efficient horizontal flight are presented and discussed. To determine an optimal trajectory, an optimal control problem is formulated and solved with the help of the trajectory optimization tool FALCON.m (FSD OptimAL CONtrol toolbox for Matlab) [12]. The results of the optimization are discussed, steady-state arcs of the trajectory are correlated with efficiency criteria. Furthermore, an appropriate design of an indicator in the cockpit display is introduced, which can deliver steady state guidance parameters to the pilot. This cockpit display is demonstrated for the climb and the horizontal flight segment of the trajectory.

2 Modelling aspects

In this section the basic modelling for the later calculation is described. The example airplane, taken as the basis of this research, is an ultralight single-seater-airplane. Flight mechanics and aerodynamic are represented by a simple point-mass-model and an aerodynamic polar. The propulsion model, including propeller, electric power train and battery, is outlined by interpolated data and generic models, parametrised with real data.

2.1 Propeller

The fixed pitch propeller is modeled using approximated wind tunnel data of thrust coefficient \(C_\mathrm{T}(J_\mathrm{Prop})\) and power coefficient \(C_\mathrm{P}(J_\mathrm{Prop})\) depending on the advance ratio

$$\begin{aligned} J_\mathrm{Prop} = \frac{V\cdot 60}{D_\mathrm{Prop}\cdot N} \end{aligned}$$
(1)

where the flight speed of the airplane is represented as V, propeller diameter as \(D_\mathrm{Prop}\) and the rotational speed as N given in the unit [\(\hbox {RPM}\;{\hat{=}}\;{2\pi }/{\hbox {min}}\)]. Figure 1 shows the measured points and the approximated regression.

Fig. 1
figure 1

Values for \(C_\mathrm{T}\) and \(C_\mathrm{P}\) of fixed pitch propeller including approximation

As the speed of the propeller is higher than \(\text {Ma}\;0.3\) at high rotational speeds, thrust coefficient \(C_\mathrm{T}\) and power coefficient \(C_\mathrm{P}\) are corrected by Prandtl–Glauert factor \(\beta = \sqrt{1-\text {Ma}^2}\) [13, 14]. Propeller thrust is then calculated with the angular velocity \(\omega =N\frac{2\pi }{60}\) and the air density \(\rho\):

$$\begin{aligned} T_\mathrm{Prop} = C_\mathrm{T}(J_\mathrm{Prop}) \cdot \rho \left( \frac{\omega }{2\pi }\right) ^2 D_\mathrm{Prop}^4 \end{aligned}$$
(2)

The therefore required shaft power is calculated with:

$$\begin{aligned} P_\mathrm{Prop} = M_\mathrm{Prop} \cdot \omega _\mathrm{Prop} = C_\mathrm{P}(J_\mathrm{Prop}) \cdot \rho \frac{\omega ^3}{(2\pi )^2} D_\mathrm{Prop}^5 \end{aligned}$$
(3)

\(M_\mathrm{Prop}\) represents the propeller momentum.

2.2 Electric motor

Rotational acceleration \(\dot{\omega }\) and resulting motor current \(I_\mathrm{Mot}\) of a motor can be modeled as a function mainly of motor voltage \(U_\mathrm{Mot}\), load torque \(M_\mathrm{load}\) and moment of inertia \(H_\mathrm{Mot}\):

$$\begin{aligned}{}[\dot{\omega },I_\mathrm{Mot}] = f(U_\mathrm{Mot},M_\mathrm{load},H_\mathrm{Mot}) \end{aligned}$$
(4)

In this work, the motor is represented as inverted (also open loop) model [15,16,17]. For the inverted motor model, motor current \(I_\mathrm{Mot}\) and voltage \(U_\mathrm{Mot}\) are represented as a function of mainly rotational speed \(\omega _\mathrm{Mot}\) and torque \(M_\mathrm{Mot}\):

$$\begin{aligned}{}[I_\mathrm{Mot},U_\mathrm{Mot}] = f(\omega _\mathrm{Mot},M_\mathrm{Mot}) \end{aligned}$$
(5)

With this simplification, the very fast dynamics of the electrical power train (in comparison to flight dynamics) can be neglected and one state in the later trajectory optimization can be saved.

For an ideal motor (motor efficiency \(\eta _\mathrm{Mot}=1\)) the current \(I_{\mathrm{id},\mathrm{Mot}}\), the motor requires, depends on torque \(M_\mathrm{Mot}\):

$$\begin{aligned} I_{\mathrm{id},\mathrm{Mot}} = \frac{M_\mathrm{Mot}}{k_\mathrm{M}} \end{aligned}$$
(6)

where \(k_\mathrm{M}\) represents the torque constant of the electric motor. With \(P_\mathrm{el}=P_\mathrm{mech}\), the required voltage \(U_{\mathrm{id},\mathrm{Mot}}\) can be calculated with the following equation:

$$\begin{aligned} U_{\mathrm{id},\mathrm{Mot}} = k_\mathrm{M} \cdot \omega _\mathrm{Mot} \end{aligned}$$
(7)
Fig. 2
figure 2

Equivalent circuit model of electric motor (refer to [9, 15, 18])

The motor model includes a resistor \(R_\mathrm{Mot}(T_\mathrm{Mot})\) which depends on the temperature of the motor \(T_\mathrm{Mot}\), and a current loss \(I_{\mathrm{loss},\mathrm{Mot}}(M,\omega )\) (refer to the equivalent circuit model shown in Fig. 2). The entire motor current sums up to

$$\begin{aligned} I_\mathrm{Mot} = I_{\mathrm{id},\mathrm{Mot}}(M) + I_{\mathrm{loss},\mathrm{Mot}}(M,\omega ), \end{aligned}$$
(8)

while the entire motor voltage \(U_\mathrm{Mot}\) evaluates to

$$\begin{aligned} U_\mathrm{Mot} = U_{\mathrm{id},\mathrm{Mot}}(\omega _{}) + I_\mathrm{Mot}{(M,\omega )} \cdot R_\mathrm{Mot}(T_\mathrm{Mot}) \end{aligned}$$
(9)

The amount of power loss \(P_{\mathrm{loss},I_\mathrm{loss}}\) due to \(I_{\mathrm{loss},\mathrm{Mot}}(M,\omega )\) is

$$\begin{aligned} P_{\mathrm{loss},I_{loss}} = I_{\mathrm{loss},\mathrm{Mot}} \cdot U_{\mathrm{id},\mathrm{Mot}} \end{aligned}$$
(10)

and includes losses caused by eddy currents, hysteresis, mechanical and air friction and residuals. Table 1 shows the assumed percentage of each loss of the nominal power \(P_\mathrm{nom}\) of the motor and—together with Eq. (11)—the mathematical correlation to the state variables of the motor, derived from [18,19,20,21].

Table 1 Power losses due to \(I_{0,\mathrm{Mot}}(M,\omega )\) [11]

All power losses due to \(I_{\mathrm{loss},\mathrm{Mot}}(M,\omega )\) sum up to

$$\begin{aligned} P_{\mathrm{loss},I_{\mathrm{loss}}} = \sum \limits _{i} k_i \cdot P_\mathrm{nom} \left( \frac{ \omega }{\omega _\mathrm{Des}} \right) ^{e_{\omega }} \left( \frac{ M}{M_\mathrm{Des}} \right) ^{e_M}, \end{aligned}$$
(11)

where \(e_{\omega}\) and \(e_\mathrm{M}\) represent the exponents shown in Table 1, \(\omega _\mathrm{Des}\) together with \(M_\mathrm{Des}\) are the design parameters of the motor.Footnote 1 For example, losses due to eddy current can be calculated with the exponents \(e_{\omega} =2\) and \(e_\mathrm{M}=2\). Therefore, Eq. (11) shows the quadratic relation between the amount of eddy current and angular velocity \(\omega\) and moment M.

The resistance of the motor depends on the temperature of the motor (e.g., [10, 17] for temperature modelling):

$$\begin{aligned} R_\mathrm{Mot}(T_\mathrm{Mot}) = R_{\mathrm{i},T_\mathrm{ref}} \cdot \left( 1 + (\alpha _{\mathrm{Cu}}(T_\mathrm{Mot}-T_\mathrm{ref})) \right) \end{aligned}$$
(12)

with the reference resistance \(R_{\mathrm{i},T_\mathrm{ref}}\) at the reference temperatureFootnote 2\(T_\mathrm{ref}\) and the fractional increase in resistivity with temperature \(T_\mathrm{Mot}\) which is \(\alpha _{\mathrm{Cu}} = {0.0039}\;{{1}/{\hbox {K}}}\) for copper. Most of all, power losses due to the resistance \(P_{\mathrm{loss},R_{\mathrm{i},\mathrm{Mot}}}=R_{\mathrm{i},\mathrm{Mot}}(T_\mathrm{Mot}) \; I_\mathrm{Mot}^2\), eddy currents \(P_{\mathrm{loss},\mathrm{eddy}}\) and hysteresis \(P_{\mathrm{loss},\mathrm{hyst}}\) create heat, which increases the motor winding’s temperature:

$$\begin{aligned} \dot{T}_\mathrm{Mot}= & {} \frac{\mathop {}\!\mathrm {d}T_\mathrm{Mot}}{\mathop {}\!\mathrm {d}t} = ( R_{\mathrm{i},\mathrm{Mot}}(T_\mathrm{Mot}) \; I_\mathrm{Mot}^2 \nonumber \\&+P_{\mathrm{loss},\mathrm{eddy}} + P_{\mathrm{loss},\mathrm{hyst}} - \dot{Q}_\mathrm{cool} ) \; \frac{1}{k_{T_\mathrm{Mot}}} \end{aligned}$$
(13)

Here, \(k_{T_\mathrm{Mot}}\) is the thermal conductivity of the windings

$$\begin{aligned} k_{T_\mathrm{Mot}} = m_{\mathrm{Mot},\mathrm{therm}} \cdot \overline{c}, \end{aligned}$$
(14)

where \(m_{\mathrm{Mot},\mathrm{therm}}\) is defined as the thermal active mass and \(\overline{c}\) as the average specific thermal conductivity of the respective materials. \(\dot{Q}_\mathrm{cool}\) represents the cooling flux, which counteracts the raise in temperature of the windings. For an air-cooled motor, the cooling flux is assumed by

$$\begin{aligned} \dot{Q}_\mathrm{cool} = \alpha _\mathrm{cool} \cdot S_\mathrm{cool} \cdot \left( T_\mathrm{Mot} - T_\mathrm{Air}(h) \right) \end{aligned}$$
(15)

with the heat transfer coefficient \(\alpha _\mathrm{cool}\) (assumed to be constant), the cooling surface \(S_\mathrm{cool}\) and the air temperature \(T_\mathrm{Air}(h)\), which decreases by \({6.5}\hbox { K}\) per \({1000}~{\hbox {m}}\) altitude.

Furthermore, the electric motor must be operated within its limitations. Those are

$$\begin{aligned} \begin{aligned} \omega _\mathrm{Mot}&\le \omega _{\mathrm{max},\mathrm{Mot}} \\ M_\mathrm{Mot}&\le M_{\mathrm{max},\mathrm{Mot}} \\ I_\mathrm{Mot}&\le I_{\mathrm{max},\mathrm{Mot}} \\ U_\mathrm{Mot}&\le U_{\mathrm{max},\mathrm{Mot}} \\ P_\mathrm{Mot}&\le P_{\mathrm{max},\mathrm{Mot}} \\ T_\mathrm{Mot}&\le T_{\mathrm{max},\mathrm{Mot}} \end{aligned} \end{aligned}$$
(16)

The electric power \(P_{\mathrm{in},\mathrm{Mot}}\), the motor takes is calculated with:

$$\begin{aligned} P_{\mathrm{in},\mathrm{Mot}}=U_\mathrm{Mot}\cdot I_\mathrm{Mot} \end{aligned}$$
(17)

2.3 Inverter

The inverter transforms the voltage and the current, provided by the battery, into the voltage and the current, required by the motor. For an ideal inverter

$$\begin{aligned} P_{\mathrm{in},\mathrm{Inv}} = P_{\mathrm{out},\mathrm{Inv}} \end{aligned}$$
(18)

is valid. The presented inverter model considers power losses due to the inverter resistance \(P_{\mathrm{loss},R_\mathrm{Inv}}=R_\mathrm{Inv} \; I_\mathrm{Bat}^2\) and switching losses \(P_{\mathrm{loss},\mathrm{sw}}\). Figure 3 shows the equivalent circuit model for the inverter and the battery.

Fig. 3
figure 3

Equivalent circuit model of inverter and battery

Equation (24) leads to a quadratic equation [9]:

$$\begin{aligned} \begin{aligned}&\left( R_\mathrm{Bat} + R_\mathrm{Inv} \right) I_\mathrm{Bat}^2 - U_{0,\mathrm{Bat}} \cdot I_\mathrm{Bat} + (P_\mathrm{Mot} + P_{\mathrm{loss},\mathrm{sw}}) = 0 \end{aligned} \end{aligned}$$
(19)

\(R_\mathrm{Inv}\) and \(R_\mathrm{Bat}\) represent the electrical resistance of the inverter and the battery, \(I_\mathrm{Bat}\) is the battery current and \(U_{0,\mathrm{Bat}}\) is defined as the no-load battery voltage. Therefore, the only mathematically appropriate solution is the following:

$$\begin{aligned} \begin{aligned} I_\mathrm{Bat}&= \frac{U_{0,\mathrm{Bat}}}{2 \left( R_\mathrm{Bat} + R_\mathrm{Inv} \right) } \\&\quad -\frac{\sqrt{U_{0,\mathrm{Bat}}^2 - 4 \; (P_\mathrm{Mot}+P_{\mathrm{loss},\mathrm{sw}}) \left( R_\mathrm{Bat} + R_\mathrm{Inv} \right) }}{2 \left( R_\mathrm{Bat} + R_\mathrm{Inv} \right) } \end{aligned} \end{aligned}$$
(20)

To obtain a real solution, the discriminant of Eq. (20) must be greater than 0:

$$\begin{aligned} {U_{0,\mathrm{Bat}}^2 - 4 (P_\mathrm{Mot}+P_{\mathrm{loss},\mathrm{sw}}) \left( R_\mathrm{Bat} + R_\mathrm{Inv} \right) } \ge 0 \end{aligned}$$
(21)

Since the inverter transforms voltage by means of PWM (Pulse Width Modulation), every switching operation between \(U_\mathrm{Inv}=U_\mathrm{Bat}\) and \(U_\mathrm{Inv}={0}\) (and vice versa) creates switching losses. Figure 4 shows temporal progressions for one switching procedure, holding a finite switching time \(t_\mathrm{sw}\). During that cycle, energy loss in one switching event is assumed to be

$$\begin{aligned} W_{\mathrm{loss},\mathrm{sw}} = \int \limits _{t_\mathrm{act}}^{t_\mathrm{act}+t_\mathrm{sw}} U(t) \cdot I(t) \mathop {}\!\mathrm {d}t . \end{aligned}$$
(22)

Providing linear progressions of voltage and current during switching event [23], energy loss in one switching procedure simplifies to \(W_{\mathrm{loss},\mathrm{sw}} = \frac{1}{6}\;\Delta U \cdot \Delta I\) with the switched voltage level \(\Delta U = U_\mathrm{sw}-{0}~{\hbox {V}}\) and the switched current level \(\Delta I = I_\mathrm{sw}-{0}~{\hbox {A}}\). Furthermore it is assumed, that the inverter contains six switches, every switching cycle includes two similar switching procedures (on and off) and that \(U_\mathrm{sw}=U_\mathrm{Bat}\) and \(I_\mathrm{sw}=I_\mathrm{Bat}\). With [24] the power losses due to switching can be estimated with

$$\begin{aligned} P_{\mathrm{loss},\mathrm{sw}} = 2 \cdot I_\mathrm{Bat} \cdot U_\mathrm{Bat} \cdot f_\mathrm{sw} \cdot t_\mathrm{sw} \end{aligned}$$
(23)

with the switching frequency \(f_\mathrm{sw}\).

Fig. 4
figure 4

Temporal progressions of voltage U, current I and resulting switching power \(P_\mathrm{sw}\) during switching procedure

Losses due to the ripple of voltage and current are not taken into account here. Furthermore, limitations through the inverter are not considered, it is assumed that the inverter is well adapted to the electric motor. Cooling and temperature of the inverter are not considered.

2.4 Battery

The battery pack is modeled by an ideal voltage source and a serial resistor (e.g., [25] and the equivalent circuit model in Fig. 3). Here, the amount of the no-load voltage source \(U_{0,\mathrm{Bat}}(\mathrm{SoC})\) depends on the actual State of Charge \(\mathrm{SoC}\). The voltage drop due to the resistance of the battery pack reduces the no-load voltage:

$$\begin{aligned} U_\mathrm{Bat} = U_{0,\mathrm{Bat}}(\mathrm{SoC}) - I_\mathrm{Bat} \cdot R_\mathrm{Bat} \end{aligned}$$
(24)

No-load voltage is correlated with the actual \(\mathrm{SoC}\) with the help of measured data, depicted in Fig. 5.

Fig. 5
figure 5

Decrease of normalized no-load voltage \(\overline{U}_{0,\mathrm{Bat}}(\mathrm{SoC})\) related to state of charge \(\mathrm{SoC}\)

The no-load voltage \(U_{0,\mathrm{Bat}}(\mathrm{SoC})\) is normalised to \(\overline{U}_{0,\mathrm{Bat}}(\mathrm{SoC})\)\(= \frac{U_{0,\mathrm{Bat}}(\mathrm{SoC})}{U_\mathrm{nom}}\), with the nominal voltage of the battery pack. In simulation, the state of charge decreases with

$$\begin{aligned} \dot{\mathrm{SoC}}=\frac{\mathop {}\!\mathrm {d}\mathrm{SoC}}{\mathop {}\!\mathrm {d}t} = \frac{I_{\mathrm{eff},\mathrm{Bat}}}{C_{\mathrm{nom},\mathrm{Bat}}} \end{aligned}$$
(25)

(e.g., [26]) in respect to the nominal capacity of the battery \(C_{\mathrm{nom},\mathrm{Bat}}\). The effective battery \(I_{\mathrm{eff},\mathrm{Bat}}\) current is slightly higher because of the Peukert–Effekt [27, 28]

$$\begin{aligned} I_{\mathrm{eff},\mathrm{Bat}} = I_\mathrm{Bat} \cdot \left( \frac{I_\mathrm{Bat}}{I_{\mathrm{nom},\mathrm{Bat}}} \right) ^{f_{\mathrm{peu}}-1} \end{aligned}$$
(26)

with the nominal battery current \(I_{\mathrm{nom},\mathrm{Bat}}\) and the Peukert-Factor \(f_{\mathrm{peu}}=1.05\) for Li-ion-batteries [28]. Recuperation by windmilling is not considered in this work. Hence, battery currents \(<0\) are not taken in account. Transient effects like considered e.g., in [10] are neither regarded.

2.5 Aerodynamic and flight mechanics

Aerodynamic forces lift L and drag D are determined by an aerodynamic polar \(C_\mathrm{D}=f(C_\mathrm{L})\) of the airplane derived from measurements. The atmospheric model is based on [29]. For all results, wind is not considered and only straight ahead flight, without turns, is regarded. Therefore, both lateral coordinates x and y are substituted to the distance s. In trajectory optimization, the equations of motion are represented by

$$\begin{aligned} \dot{s}&= V \cdot \cos \gamma \end{aligned}$$
(27)
$$\begin{aligned} \dot{h}&= V \cdot \sin \gamma \end{aligned}$$
(28)
$$\begin{aligned} \dot{V}&= \frac{{T}-{D}}{m}-g \cdot \sin {\gamma } \end{aligned}$$
(29)
$$\begin{aligned} \dot{\gamma }&= \frac{L}{m\cdot V}- \frac{g\cdot \cos \gamma }{V} \end{aligned}$$
(30)

with the altitude h, the velocity V, the flight path angle \(\gamma\), the acceleration due to gravitation g and the weight of the airplane m (cf. e.g., [30, 31]). Therefore, it is assumed that thrust is co-linear and lift is perpendicular to the flight direction (V).

Considering a steady-state calculation without acceleration \(\dot{V} = {0}\) and with \(\dot{\gamma } = {0}\), the flight path angle is calculated by

$$\begin{aligned} \sin \gamma = \frac{T-D}{m\cdot g} \end{aligned}$$
(31)

Wind is not considered in this work.

3 Optimization methods

In this work, two Optimization methods are applied. To obtain an optimal trajectory from the model, a typical optimal control problem is formulated and solved. With the help of a steady-state evaluation of the model [with Eq. (31)] steady-state model outputs (e.g., \(\gamma\), \(I_{\mathrm{eff},\mathrm{Bat}}\), etc.) can be determined. This second approach is performed to calculate the values of the efficiency criteria. Both approaches are explained next.

3.1 Optimal control problem

The current optimization problem can be modeled as a classical optimal control problem [11]. For a given distance \(s_\mathrm{f} = {220}\,{\hbox {km}}\), the trajectory should reach the final boundary condition with a maximum amount of State of Charge \(\mathrm{SoC}\):

Find

$$\begin{aligned} \max \; J = \mathrm{SoC}(t_\mathrm{f}) \end{aligned}$$
(32)

with respect to the dynamic constraints

$$\begin{aligned} \dot{\mathbf {x}} = f(\mathbf {x}(t),\mathbf {u}(t),\mathbf {p}), \end{aligned}$$
(33)

with the state vector

$$\begin{aligned} \mathbf {x} = \begin{pmatrix} {s}\\ {h}\\ {V}\\ {\gamma }\\ {T_\mathrm{Mot}}\\ {\mathrm{SoC}} \end{pmatrix} \end{aligned}$$
(34)

and the state derivatives from Eqs. (27), (28), (29), (30), (13) and (25). The control vector contains the lift coefficient \(C_L\) and rotational speed of the propeller N:

$$\begin{aligned} {\mathbf {u}} = \begin{pmatrix} {C_L}\\ {N} \end{pmatrix} \end{aligned}$$
(35)

As the trajectory optimization will be performed in 3 phases, the parameter vector contains two transition times (\(t_{12}\) and \(t_{23}\)) and the final time \(t_\mathrm{f}\) at \(s_\mathrm{f}\):

$$\begin{aligned} \mathbf {p} = \begin{pmatrix} t_{12} \\ t_{23} \\ t_\mathrm{f} \end{pmatrix} \end{aligned}$$
(36)

Besides, the initial conditions

$$\begin{aligned} \begin{aligned} s_{t_0}&= {0}\hbox { m}, \\ h_{t_0}&= {500}\hbox { m},\\ \mathrm{SoC}_{t_0}&= {95}\,{\%}, \end{aligned} \end{aligned}$$
(37)

the final condition

$$\begin{aligned} \begin{aligned} s_{t_\mathrm{f}}&= s_\mathrm{f} \\ h_{t_\mathrm{f}}&= {500}\hbox { m} \end{aligned} \end{aligned}$$
(38)

and the path constraints

$$\begin{aligned} V_\mathrm{IAS,min}\;&\leq & V_\mathrm{IAS}\;&\leq\; V_\mathrm{IAS,max}\\ 0,9\;&\leq & n_\mathrm{Z}\;&\leq\; 1,1\\ && P_\mathrm{Mot}\;&\leq\; P_\mathrm{Mot,max}\\ && I_\mathrm{Mot}\;&\leq\; I_\mathrm{Mot,max}\\ 0\;&\leq & T_\mathrm{Mot}\;&\leq\; T_\mathrm{Mot,max}\\ 0\;&\leq & h(t)\;&\leq\; {3000}\hbox{ m}\\ & & I_\mathrm{Mot}\;&\leq\; I_\mathrm{Mot,max}\\ 0\;&\leq & I_\mathrm{Bat}\;&\leq\; I_\mathrm{eff,Bat,max}\\ && \texttt{Disc}\;&\leq\; 1\\ 0\;&\leq & \Delta h(s)&\\ \end{aligned}$$
(39)

are not to be violated. The load factor is defined as the ratio \(n_\mathrm{Z} = \frac{L}{m\;g}\). It is limited to smoothen transitions between flight phases. The path constraint Disc represents the discriminant from Eqs. (20) and (21) as a purely mathematical constraint, while \(\Delta h(s) = h(s)-h_\mathrm{terr}(s)\) enforces the trajectory on altitudes higher than a fictional terrain.Footnote 3 In phase 2, the flight path angle is further set to \(\gamma = 0\) to enforce a horizontal flight segment.

The trajectory optimization is performed using the Matlab toolbox FALCON.m [12]. This toolbox helps the user to define the optimal control problem in Matlab and solves the problem very efficiently. FALCON.m applies trapezoidal collocation as transcription and discretization scheme. The resulting classical parameter optimization problem is passed to IPOPT (Interior Point OPTimizer), which solves the problem and returns the result to FALCON.m back again.

3.2 Steady state evaluation

The steady state evaluation is performed using evaluation points on a grid. Therefore, Matlab provides fast “element wise operations” [32] to compute large grids of evaluation points. The steady state model inputs (airspeed \(V_{\mathrm{TAS}}\) and rotational speed N) are arranged in respectively one matrix, so that every combination of input value exists once. With the mathematical model described above [especially Eq. (31)], the model outputs can be calculated at every evaluation point. The model outputs contain the battery current \(I_{\mathrm{eff},\mathrm{Bat}}\) and the flight path angle \(\gamma\) corresponding to the evaluation point. Those values in model output grids who violate a constraint [ref. to Eqs. (16), (39)] are replaced with an error symbol \(\zeta\) and are therefore not considered in the following analysis.Footnote 4

Fig. 6
figure 6

Example for a criterion depending on two parameters P1 and P2

In the next step, efficiency criteria as \(\max \left( \frac{V}{I_{\mathrm{eff},\mathrm{Bat}}} \right)\) (described in later Sect. 4) are computed. Figure 6 shows the values of an exemplary criterion which depends on two parameters P1 and P2 in 3D. By projecting the criterion surface in Fig. 6 into the P1 and P2 plane, the criterion’s silhouette depicts the relation of the respective parameter and the criterion (Fig. 7).

Fig. 7
figure 7

Projection of criterion into the P1 and P2-plane

The black lines in Fig. 7 show the maximum values of the criterion depending in corresponding parameter. The maximum values of the criterion are marked with \(\diamondsuit\). Furthermore, both \(\Delta\) marks depict the range in which the “loss” of criterion’s value is \(<{5}\,{\%}\).

Considering the horizontal flight criterion’s values, first combinations of steady-state controls \(V_{\mathrm{IAS}}\) and N must be found in the calculation grid, that cause a flight path angle \(\gamma \approx {0}^{\circ }\). Therefore, \(\min |\gamma |\) is searched with the method presented in Sect. 3.2. By means of the acquired indices fulfilling the constraint \(\min |\gamma |\), the criterion’s values can be analysed subsequently.

4 Flight performance preliminaries

In the following section basic efficiency criteria are shortly presented and discussed. All criteria can only be applied on steady-state arcs of flight without acceleration.

4.1 Criterion of range optimal horizontal flight

To implement a most-efficient horizontal steady-state cruise flight, the distance travelled \(\Delta s\) has to be maximized with a smallest possible amount of \(\Delta \mathrm{SoC}\):

$$\begin{aligned} \text {max }\left( \frac{\Delta s}{\Delta \mathrm{SoC}} \right) _{\gamma =0} \end{aligned}$$
(40)

For a small \(\mathop {}\!\mathrm {d}t\), the criterion reads [with Eq. (25)]:

$$\begin{aligned} \text {max }\left( \frac{\Delta s}{\Delta \mathrm{SoC}} \frac{\frac{\mathop {}\!\mathrm {d}}{\mathop {}\!\mathrm {d}t}}{\frac{\mathop {}\!\mathrm {d}}{\mathop {}\!\mathrm {d}t}} \right) _{\gamma =0} = \text {max }\left( \frac{V}{I_{\mathrm{eff},\mathrm{Bat}}} \right) _{\gamma =0} \end{aligned}$$
(41)

In Fig. 8 maximum values of the criterion related to the flight altitude h are depicted. The corresponding steady-state flight guidance parameters airspeed \(V_{\mathrm{IAS}}\) and rotational speed N of the propeller are drawn next to it. While the rotational speed N increases with the altitude, the optimal airspeed \(V_{\mathrm{IAS}}\) is quasi constant. The values of the criterion itself decreases slightly with increasing altitude (approx. \({1.5}{\%}\)). An optimal horizontal cruise flight is slightly more efficient at lower altitudes.

Fig. 8
figure 8

Evaluation of the criterion of range optimal horizontal flight

The dashed lines in Fig. 8 show the influence of guidance parameter variation on the value of the criterion.Footnote 5 The inner lines mark a loss of \({2,5}{\%}\), whereas the outer lines one of \({5}{\%}\) of the maximum of the criterion (these marks can be seen in the cockpit display later).

Nevertheless, the maximum amount of the criterion is nearly constant over the altitude. This can be interpreted as a quasi linear dependence between required power in horizontal flight and the consumption equivalent \(I_{\mathrm{eff},\mathrm{Bat}}\). Figure 9 depicts the dependence derived from the presented model within the interval of required power for horizontal flight between \(h={0}\hbox { m}\) and \(h={3000}\hbox { m}\). Obviously, the quadratic approximation is very similar to the linear approximation. This demonstrates the described quasi linear dependence and implies a nearly constant efficiency of the power train within the power spectrum required for horizontal flight between \(h={0}\hbox { m}\) and \(h={3000}\hbox { m}\) (similar approach can be found in [33]).

Fig. 9
figure 9

Consumption equivalent \(I_{\mathrm{eff},\mathrm{Bat}}\) versus required power \(P_\mathrm{hor}\) in horizontal flight between \(h={0}\hbox { m}\) and \(h={3000}\hbox { m}\)

4.2 Criterion of range optimal climb

To realize a preferably efficient steady-state flight, but without constraining the flight path angle \(\gamma\) to zero, it is possible to gain altitude first and to increase the travelled distance by gliding without energy input afterwards. Figure 10 shows both phases. It is assumed, that the glide phase II is performed within the best glide ratio \(E_\mathrm{max}=\tan (\gamma _\mathrm{gl})\). The entire travelled distance sums up with both phases horizontal parts to

$$\begin{aligned} \Delta s_\mathrm{tot} = \Delta s_{\mathrm{cl},\mathrm{horiz}} + \Delta s_{\mathrm{gl},\mathrm{horiz}} \end{aligned}$$
(42)

From

$$\begin{aligned} \Delta s_{\mathrm{cl},\mathrm{horiz}}=\Delta s_\mathrm{cl}\cdot \cos (\gamma _\mathrm{cl}) \end{aligned}$$
(43)

and

$$\begin{aligned} \Delta s_{\mathrm{gl},\mathrm{horiz}}=\Delta s_\mathrm{cl}\cdot \sin \gamma _\mathrm{cl}\cdot E_\mathrm{max} \end{aligned}$$
(44)

the criterion of Range Optimal Climb can be derived for a small \(\mathop {}\!\mathrm {d}t\) [cf. Eq. (41)]:

$$\begin{aligned} \text {max }\left( \left( \cos \gamma + \sin \gamma \cdot E_\mathrm{max} \right) \frac{V}{I_{\mathrm{eff},\mathrm{Bat}}}\right) \end{aligned}$$
(45)

As shown above for the criterion of Range Optimal Horizontal Flight, Fig. 11 depicts the maximum values of the criterion of Range Optimal Climb related to the flight altitude h. The corresponding steady-state flight guidance parameters airspeed \(V_{\mathrm{IAS}} = \sqrt{\frac{\rho (h)}{\rho _0}}V_{\mathrm{TAS}}\) and rotational speed N of the propeller are shown next to it. The maximum values for both efficiency criteria are plotted for comparison, including the corresponding guidance parameters. For the sake of simplicity, the sensitivity of \({2.5}{\%}\) was omitted here. It can be derived, that a flight within the criterion of a Range Optimal Climb is more efficient than a flight within the criterion of a Range Optimal Horizontal Flight, because the values for the criterion of Range Optimal Climb are slightly larger than the values for the criterion of Range Optimal Horizontal Flight. The airspeed \(V_{\mathrm{IAS}}\) and the rotational speed N is located slightly higher for the criterion of a Range Optimal Climb and the sensitivity of \({5}{\%}\) includes a larger interval of guidance parameters.

Fig. 10
figure 10

Criterion of range optimal climb [11]

Fig. 11
figure 11

Evaluation of the criterion of range optimal climb.

The lines of optimal airspeed \(V_{\mathrm{IAS}}\) and rotational speed N proceed less smooth for the criterion of Range Optimal Horizontal Flight than for the criterion of a Range Optimal Climb. The reason for this is the resolution of the grid points and the evaluation for horizontal flight described in Sect. 3.2. Unlike the criterion of a Range Optimal Climb, the criterion of a Range Optimal Horizontal Flight requires the minimization of \(|\gamma |\) first, which causes an additional constraint. Hence, the finer the resolution the smoother are the lines.

5 Results

In the following section, the major results of the trajectory optimization are presented, like the energy optimal trajectory (cf. Fig. 12). Thereby, steady-state arcs are correlated with efficiency criteria and the cockpit display concept is demonstrated.

Fig. 12
figure 12

Optimal trajectory for \(s_\mathrm{f} = {220}\hbox { km}\)

Figure 12 shows the histories of the states (h, \(V_{\mathrm{TAS}}\), \(\gamma\), \(T_\mathrm{Mot}\) and \(\mathrm{SoC}\)), the controls (\(C_\mathrm{L}\) and N) and the constraints (\(V_{\mathrm{IAS}}\), \(J_\mathrm{Prop}\), \(n_\mathrm{Z}\), \(P_\mathrm{Mot}\), \(I_\mathrm{Mot}\), \(I_{\mathrm{eff},\mathrm{Bat}}\) and \(\texttt {Disc}\)) plotted over the flown distance s. All dashed lines represent minimum or maximum values of the path constraints described in (39). The path constraint due to terrain is shown in the progression of state h.

Because of the terrain constraint, the trajectory is structured in 3 phases (consistent with the three phases described in Sect. 3.1), a climb, a horizontal and a descend segment. In phase 2, the constraint \(\gamma = 0\) forces the trajectory to a perfect horizontal flight. Without a terrain constraint the solution of the optimal control problem would be a trivial nearly straight flight towards the final point at \(s_\mathrm{f}\). The reason for this is the general dependence of altitude and the efficiency shown in Fig. 8. Therefore, a flight in lower altitude is slightly more efficient than in higher altitude. Due to the correlation of the altitude h and the air density \(\rho\), the optimal flight speed \(V_{\mathrm{TAS}}\) generally rises with altitude. In contrast, airspeed \(V_{\mathrm{IAS}}\) is more or less constant during the steady-state arcs.

During the first half of the climb phase, the motor temperature \(T_\mathrm{Mot}\) increases because of the high motor power \(P_\mathrm{Mot}\) required in the climb segment. At the same time, the motor cooling effect increases, because the air temperature \(T_\mathrm{Air}\) decreases with the altitude. The motor temperature constraint is not violated during the whole trajectory. Hence, the cooling system of the motor is sufficient. Only in the transient phases of the trajectory, the constraint in \(n_\mathrm{Z}\) is reached. During the residual trajectory, all values are within the constraints. The mathematical constraint \(\texttt {Disc}>0\) is not violated at all, either. The trajectory ends with a remaining \(\mathrm{SoC}(t_\mathrm{f}) = {30,1}{\%}\).

5.1 Climb segment

During the climb segment of the optimal trajectory, the flight path angle \(\gamma (t)\) is constant. The Airspeed \(V_{\mathrm{IAS}}\) is constant (just as lift coefficient \(C_L\)), while \(V_{\mathrm{TAS}}\) increases with the rising altitude. The motor moment \(M_\mathrm{Mot}\) (not depicted) is constant, thus the rising motor power \(P_\mathrm{Mot}\) is caused by the increasing rotational speed N in the climb segment. Hence, battery current \(I_{\mathrm{eff},\mathrm{Bat}}\) rises as well.

Fig. 13
figure 13

Criterion of range optimal climb depicted over \(V_{\mathrm{IAS}}\) in the climb segment

With the approach described in Sect. 3.2, Fig. 13 was generated for the criterion of Range Optimal Climb. It shows the criterion’s maximum values in relation to the airspeed \(V_{\mathrm{IAS}}\). Analogously, Fig. 14 shows the same for the second steady state guidance parameter, rotational speed N. Both figures are generated for the altitude \(h = {1883}\hbox { m}\) in the climb segment of the optimal trajectory.

For comparison of the operation points of the trajectory and the efficiency criteria, the figures include the actual operating point (marked with \(\pentagon\)) in \(h = {1883}\hbox { m}\). Additionally, sensitivities of \({2.5}{\%}\) and \({5}{\%}\) (ref. to Fig. 7) of deviation of the maximum value of the criterion are plotted with \(\triangledown\) and \(\Delta\). The optimal airspeed regarding the efficiency criterion of Range Optimal Climb is \(V_{\diamondsuit ,\mathrm{IAS}}={39.9}\;{{\mathrm{m}}}/{\mathrm{s}}\), while the operation point in the trajectory is \(V_{\pentagon ,\mathrm{IAS}}={39.7}\;{\mathrm{m}}/\mathrm{s}\) (difference: \({0.5}{\%}\)). Similarly, the optimal rotational speed is \(N_{\diamondsuit }={2778}\hbox { RPM}\), while the operation point in the trajectory is \(N_{\pentagon}={2728}\;{\mathrm{RPM}}\) (difference: \({1.8}{\%}\)). The deviation for the rotational speed is slightly bigger than for the optimal airspeed, but the “loss” of criterion only amounts to \({0.5}{\%}\).

Fig. 14
figure 14

Criterion of range optimal climb depicted over N in the climb segment

The upper sensitivity for \({5}{\%}\) does not appear in Fig. 14, because the related rotational speed is located at \(>N_\mathrm{max}(={3000}\hbox { RPM})\). Considering the small differences between steady-state evaluated value \(\diamondsuit\) and operation point in the trajectory \(\pentagon\), both can be correlated. That means, for a Range Optimal Climb, there exists a steady-state efficiency criterion and corresponding guidance parameters can be determined without trajectory optimization.

5.2 Cockpit display in the climb segment

A pilot who tries to guide his airplane efficiently in any climb segment, should try to control his airplane within the airspeed \(V_{\mathrm{IAS}}\) and the rotational speed N, which are derived from the efficiency criterion of Range Optimal Climb. Therefore, a new cockpit display is introduced here (Fig. 15). Both indicators—for airspeed and rotational speed—contain additional small green and blue beams and a white line. The white line shows the optimal operation parameter derived from the criterion. The green and blue beams limit the sectors within the “loss” of \({2.5}{\%}\) and \({5}{\%}\) of criterion’s maximum. With those beams, the pilot is shown the sensitivity of the corresponding steady-state guidance parameter. The close location of upper bounds of \({2.5}{\%}\) and \({5}{\%}\) described above is easy to see in Fig. 15 Both actual values of steady-state parameters are very close to the maximum white line marks of the display.

Fig. 15
figure 15

Cockpit display with indicators for the criterion of range optimal climb taken from climb segment in the optimal trajectory in \(h_\mathrm{ISA}={1883}\hbox { m}\) (display is further developed using basic code by [34])

5.3 Horizontal segment

In the horizontal flight nearly every parameter in Fig. 12 is constant. Solely, the \(\mathrm{SoC}\) decreases persistently. Since the battery no-load voltage \(U_{0,\mathrm{Bat}}(\mathrm{SoC})\) decreases with descending \(\mathrm{SoC}\) (Fig. 5) and the required power is constant, the battery current \(I_\mathrm{Bat}\) increases slightly.

Fig. 16
figure 16

Criterion of range optimal horizontal flight depicted over \(V_{\mathrm{IAS}}\) in the horizontal segment

Fig. 17
figure 17

Criterion of range optimal horizontal flight depicted over N in the horizontal segment

Figures 16 and 17 show—in the same manner as above—the criterion of a Range Optimal Horizontal Flight over airspeed \(V_{\mathrm{IAS}}\) and rotational speed N. As shown above, both figures include the actual operating point (marked with \(\pentagon\)) in the trajectory. Additionally, sensitivities of \({2.5}{\%}\) and \({5}{\%}\) of deviation of the maximum value of the criterion are plotted with \(\triangledown\) and \(\triangle\). The operating point in the trajectory is even closer to the maximum of the criterion of a Range Optimal Horizontal Flight compared to above. The optimal airspeed is \(V_{\diamondsuit ,\mathrm{IAS}}={39.0}\;{\hbox {m}}/{\mathrm{s}}\), while the operation point in the trajectory is \(V_{\pentagon ,\mathrm{IAS}}={39.1}\;{\mathrm{m}}/\mathrm{s}\) (difference: \({0.3}{\%}\)). The optimal rotational speed \(N_{\diamondsuit }={2622}\hbox { RPM}\), while the operation point in the trajectory is \(N_{\pentagon }={2627}\hbox { RPM}\) (difference: \({0.2}{\%}\)). In comparison to the climb segment of the trajectory, both guidance parameters correlate very well to the parameters taken from the criterion. Here the “loss” of criterion between steady-state evaluated value \(\diamondsuit\) and the actual operation point \(\pentagon\) in the trajectory only amounts to \({0.2}{\%}\). Therefore, the steady-state criterion and operation points in the trajectory can be correlated even better.

In Fig. 17 a special phenomenon can be seen: for rotational speeds below approx. \({2450}\hbox { RPM}\), there is an ambiguity in the evaluation of the model.

Figure 18 describes how this ambiguity occurs.Footnote 6 Neglecting the angle between thrust and velocity vector, in steady-state horizontal flight thrust T equals aerodynamic drag D:

$$\begin{aligned} T = D \end{aligned}$$
(46)

Thus, the horizontal operation points are always located at intersections of drag and thrust curve. For rotational speeds \(>{2450}\hbox { RPM}\), the maximum thrust curve reaches the power limit \(P_{\mathrm{max},\mathrm{Mot}}\) of the motor model. Therefore, the rotational speed at the intersection in area II is slightly lower than in area I. As a consequence, ambiguous values of the criterion do not exist. Maximum thrust curves for rotational speeds \(<{2450}\hbox { RPM}\) intersect the drag curve two times before reaching the maximum power limit \(P_{\mathrm{max},\mathrm{Mot}}\). Here, two horizontal operation points exist for each rotational speed, the ambiguity occurs in Fig. 17.

Fig. 18
figure 18

Phenomenon of ambiguity in horizontal flight

The described phenomenon appears to be caused by a double solution comparable to a root search. For further development of the steady-state evaluation of the model, a mechanism which captures the non useful solution is proposed. Since the type of the model is analytical/numerical mixed, the non useful solution can only be identified numerically. An additional examination of every evaluation point and its neighbouring points could lead to the desired result. Alternatively, one could calculate the derivative of the criterion with respect to rotational speed and the airspeed and use this information to find the desired point.

5.4 Cockpit display in the horizontal segment

Figure 19 shows the cockpit display for an operating point taken from the horizontal segment of the trajectory (red mark).

Fig. 19
figure 19

Cockpit display with indicators for the criterion of range optimal horizontal flight taken from horizontal segment in the optimal trajectory in \(h_\mathrm{ISA}={2550}\hbox { m}\)

It is obvious, that the operation point in the trajectory is located exactly at the white lines of the optimum of the criterion of the Range Optimal Horizontal Flight.

For rotational speeds \(>{2450}\hbox { RPM}\), there is only one operation point (see above). Every airspeed \(V_{\mathrm{IAS}}\) requires one unique rotational speed N to perform a horizontal flight without acceleration. Only one of both steady-state parameters must be displayed in the cockpit. However, for the sake of completeness both parameters are displayed.

5.5 Descend segment

Also the descend segment of the optimized trajectory shows nearly constant parameters. The airspeed is \(V_{\mathrm{IAS}}\,={38.9}\;\mathrm{m}/\mathrm{s}\), the flight path angle is \(\gamma ={-1.7}^{\circ }\). The motor power \(P_\mathrm{Mot}\) and the battery current \(I_{\mathrm{eff},\mathrm{Bat}}\) are \(>~0\). Hence, the descend segment is performed with an amount of energy input, because the flight path angle is \(\gamma > \arctan (E_\mathrm{max})\). Hence, the descend segment is not a pure glide flight.

Fig. 20
figure 20

Aerodynamic polar without and with the influence of the propeller

Figure 20 shows the aerodynamic polar (thin line) without propeller. It includes an additional polar with the effect of the propeller, when \(I_{\mathrm{eff},\mathrm{Bat}}\overset{!}{=}{0}\) is claimed (bold line). The corresponding best glide ratio of both polars are marked with \(\diamondsuit\). Because of the additional drag produced by the propeller, the polar assuming \(I_{\mathrm{eff},\mathrm{Bat}}\overset{!}{=}{0}\) is located a bit more right at higher values of \(C_\mathrm{D}\). Therefore, the propeller is driven by a windmilling-effect. Figure 20 additionally includes the actual operation point in the descend segment of the trajectory (marked with \(\pentagon\)). This operation point shows nearly the same lift coefficient like the mark \(\diamondsuit\) for the best glide ratio with respect to \(I_{\mathrm{eff},\mathrm{Bat}}\overset{!}{=}{0}\) (the difference between \(C_{\mathrm{L},\diamondsuit }\) and \(C_{\mathrm{L},\pentagon }\) is \({2}{\%}\)). The effective drag coefficient \(C_{\mathrm{D},\mathrm{eff},\pentagon }\) of the actual operation point is located left of both polars, because there is an energy input in the system during the descend segment. It evaluates to:

$$\begin{aligned} C_{\mathrm{D},\mathrm{eff}} = \frac{2 \cdot (-T_{\mathrm{desc}}+D_{\mathrm{desc}})}{\rho \cdot V_{\mathrm{desc}}^2\cdot S_{\text {ref}}} \end{aligned}$$
(47)

Since \(C_{\mathrm{L},\diamondsuit }\approx C_{\mathrm{L},\pentagon }\), the optimal airspeed \(V_{\mathrm{IAS}}\) corresponds to the lift coefficient of best glide regarding \(I_{\mathrm{eff},\mathrm{Bat}}\overset{!}{=}0\). A descent flight whose flight path angle is smaller than the glide angle within \(E_\mathrm{max}\) should be guided with in the airspeed of best glide ratio \(V_{\mathrm{IAS},E_\mathrm{max}}\) with respect to \(I_{\mathrm{eff},\mathrm{Bat}}\overset{!}{=}{0}\).

6 Summary and outlook

This study has developed a mathematical model for a small airplane including its battery-electrical power train. This model allows to perform both a trajectory optimization and a steady-state evaluation. Steady state efficiency criteria for electrical aircraft have been set up for a Range Optimal Horizontal Flight and for a Range Optimal Climb. Their evaluation has shown, that a horizontal flight is slightly more efficient in lower altitudes. Furthermore, the values for the criterion of a Range Optimal Climb are always slightly bigger than for a Range Optimal Horizontal Flight. The problem of trajectory optimization has been defined and then solved using the Matlab toolbox FALCON.m. Its results have been discussed and analysed. Hereby, for the climb and the horizontal segment, the efficiency criteria of Range Optimal Climb and Range Optimal Horizontal Flight could be correlated with the operation points in the trajectory. A special phenomenon of the criterion’s values in horizontal flight has been explained.

Future research will refine the mathematical model further. For example in the battery model, mid-term voltage behaviour on the battery current could be examined. Effects on the criteria caused by wind will be analysed, as well as the additional influence of a variable pitch propeller model. Moreover, thermal behaviour should be added to the models of inverter and battery. This will bring new states and constraints to the optimization. Maximum temperature of the inverter probably can be reached very fast in the trajectory, because the thermal mass of the inverter is small. Thus, an active cooling for the inverter could be necessary. Depending on the aircraft design, cooling of the battery may be complex.

One of the fundamental questions is if these results can be transferred to bigger aircraft.