EVTOL Tilt-Wing Aircraft Design under Uncertainty Using a Multidisciplinary Possibilistic Approach

Abstract

Recent development in Electric Vertical Take-off and Landing (eVTOL) aircraft makes it a popular design approach for urban air mobility (UAM). When designing these configurations, due to the uncertainty present in semi-empirical estimations, often used for aerodynamic characteristics during the conceptual design phase, results can only be trusted to approximately 80% accuracy. Accordingly, an optimized aircraft using semi-empirical estimations and deterministic multi-disciplinary design optimization (MDO) approaches can be at risk of not being certifiable in the detailed design phase of the life cycle. The focus of this study was to implement a robust and efficient possibility-based design optimization (PBDO) method for the MDO of an eVTOL tilt-wing aircraft in the conceptual design phase, using existing conventional designs as an initial configuration. As implemented, the optimization framework utilizes a deterministic gradient-based optimizer, run sequentially with a possibility assessment algorithm, to select an optimal design. To achieve this, the uncertainties which arise from multi-fidelity calculations, such as semi-empirical methods, are considered and used to modify the final design such that its viability is guaranteed in the detailed design phase. With respect to various requirements, including trim, stability, and control behaviors, the optimized eVTOL tilt-wing aircraft design offers the preferred results which ensure that airworthiness criteria are met whilst complying with predefined constraints. The proposed approach may be used to revise currently available light aircraft and develop eVTOL versions from the original light aircraft. The resulting aircraft is not only an optimized layout but one where the stability of the eVTOL tilt-wing aircraft has been guaranteed.

1. Introduction

Generally, two configurations with respect to cruise flight and vertical take-off and landing (VTOL) capabilities are of interest to urban air mobility (UAM) applications: rotary-wing and fixed-wing vehicles. Flying vehicles that primarily generate lift using the rotary-wings concept are naturally limited in cruise speed and tend to have reduced range capabilities. They are, however, superior when considering hover and VTOL capabilities. Comparatively, fixed-wing designs are less effective than rotors in hover flight, but they are more efficient and faster during cruise flight, providing higher attainable ranges. Currently, active research is ongoing to enhance the existing knowledge of fossil fuel-based aircraft, intending to eventually replace them with electric vertical take-off and landing (eVTOL) airplanes [1,2,3].

Currently, light electric aircraft design and development is focused on battery technology, electrical power systems, and electrical system architecture [4]. Battery system modelling in multidisciplinary analysis and optimization is a new field where insufficient research is available for fuel cell sub-system modelling. Research related to electric aircraft often focuses on the performance of the aircraft that is required to achieve the theoretical capabilities of the fuel cell system; essentially, fuel cell systems are scaled regarding the fuel cell output power to avoid considering fuel cell sub-systems [5,6,7,8,9,10]. On this topic, NASA has published a series of early studies on the feasibility of fuel cells in aircraft design [8,9,10]. Furthermore, NASA has provided a scalable high-fidelity solid oxide fuel cell (SOFC) model to design fuel cell powerplants [11]. Unfortunately, this analysis cannot be easily implemented for aircraft design purposes. In addition, some research has been accomplished to provide some early-stage methods for the conceptual design phase of electric aircraft [12,13,14,15,16,17,18]. Overall, the electric aircraft design literature that has been accomplished so far focuses on technological feasibility, sensitivity, or sizing and synthesis [19,20,21,22,23].

Computer-aided engineering CAx techniques such as CAD, CAM, CAE and so on are becoming more popular for the purpose of aircraft design and development. Different methods have been developed over the years [24,25,26,27]. Two classes of optimization algorithms can be considered when investigating aircraft optimization: stochastic algorithms, reliant upon randomization to conduct global exploration, and gradient-based algorithms, which guide the design into a local optimum using the derivatives of the objective function. The former is capable of finding the global optimum but does so at a great computational cost [28], while the latter can be guided only to a local minimum [29,30]. When optimizing the design of an aircraft, which includes very expensive analyses, gradient optimizers are typically used to ensure that an improved design can be found in a reasonable amount of time [31]. Here, Sequential Quadratic Programming (SQP) is employed, which is an approach to nonlinear constrained optimization built upon Lagrange optimization [32,33].

Inherent to any computational modelling approach is some margin of error, an interval that limits the extent to which calculated results can be trusted. This holds true for the estimations of semi-empirical methods which can be trusted to have approximately 80% precision with the initial geometry resulting from the initial airplane layout [34]. Ignoring this deviation brings a failure risk for when the project transitions from theory to reality; a risk that can only be reduced by greater fidelity estimations. Beyond this, even physical modelling is afflicted by random parameters such as the environment and human error. Sources of uncertainty can be modelled using probability theory provided there is sufficient data to determine a probability distribution shape and parameters. Non-probabilistic methods can be implemented to model uncertainty in cases with insufficient data to develop a precise determination of the probability density function (PDF) shape or parameters. Some of the methods are Interval and Fuzzy Modelling, Convex Modelling and Evidence Theory [35,36].

Possibility-Based Design Optimization (PBDO) can be employed to account for the unavoidable uncertainty that appears in design optimization. This can be implemented in cases where inadequate data is present to devise probability density functions but still produces a trusted configuration using membership functions for epistemic uncertainties [37]. Visualized in Figure 1, the PBDO method searches for configurations only in regions where a design can be trusted to be feasible. This is achieved by defining an interval function for each input, where the input is trusted to lie inside, and prohibiting this interface functions from exiting the feasible design space. Some strategies exist for PBDO including the vertex method, which solves full optimizations for all upper and lower boundaries for the uncertain parameters and variables. However, it becomes extremely expensive for cases with numerous uncertainties [38].

2. Methodology

The proposed eVTOL tilt-wing aircraft has a conventional configuration equipped with a rotating wing and horizontal tail, as well as leading-edge open rotors, four passenger seats and one pilot and is capable of reaching a maximum range of 300 km. The transition between horizontal and vertical flight is achieved by the rotation of the horizontal tail and the wing. Certification was accounted for by following the Special Condition for small-category VTOL aircraft by the EASA light aircraft certification requirements (CS-23) and small rotorcrafts (CS-27). The final geometry model of the optimized eVTOL tilt-wing aircraft is shown in Figure 2. The dimension details of the proposed eVTOL tilt-wing aircraft are presented in the Results.

In this case, the main concern was sizing the empennage and wing, where the optimization was completed by altering the defining geometry whilst satisfying the requirements of the design given by the aircraft geometry, electrical propulsion characteristics, and different flight circumstances. The involved disciplines were aerodynamics, weight and balance, propulsion, as well as stability and control. An initial sizing of the electric propulsion system was accomplished to estimate the required number of engines along with the number, mass, and volume of batteries.

2.1. MAPLA

The semi-empirical analysis program MAPLA, devised for the optimization of light, general aviation aircraft, was enhanced for use with electric vertical take-off and landing (eVTOL) vehicles. The original implementation includes five disciplines: aerodynamics, propulsion, performance, weight and balance as well as stability and control, all built with state-of-the-art analytical procedures and design data collections combined into a fully automated method. For this investigation, the electric propulsion submodule was developed including the following considerations: engine/propeller positioning and sizing, off-design analysis of blades, noise analysis, mass calculation, volume determination, and battery specifications.

In a series of previous investigations, this tool was demonstrated to be capable of modelling the characteristics of small aircraft with acceptable precision [39,40,41,42,43]. In these cases, the modelled aircraft were all using conventional configurations, a requirement imposed by the semi-empirical methods which underpin many of MAPLA’s analyses. These limitations are still present in the current version of the tool, and as such only eVTOL aircraft which resemble a conventional aircraft while in cruise can be successfully modelled; atypical or hybrid configurations cannot be analyzed with MAPLA.

For the purpose of this study, initially, a small propeller-driven aircraft was chosen and first analyzed by MAPLA. In the next step, the corresponding aerodynamic characteristics predicted by MAPLA were compared with those of the wind tunnel test results for the small aircraft [44,45] to validate the calculated results. Figure 3a,b shows the geometry of the original aircraft model and the modelled aircraft using MAPLA, respectively. The general characteristics of the airplane and properties of the investigated flight condition are presented in

Table 1. General characteristics of the aircraft used for MAPLA’s validation.

Variable	Value	Unit
Wing span	10.97	$\mathrm{m}$
Mean Aerodynamic Chord	1.511	$\mathrm{M}$
Wing Surface Area	16	${\mathrm{m}}^2$
Aspect Ratio	7.52	-
Mach	0.25	-
Mass	980	$\mathrm{k}\mathrm{g}$
CG	10	%
Altitude	0	m

.

In the following, in Figure 4a–c, the longitudinal aerodynamic characteristics of the investigated airplane are presented for the corresponding flight condition. Figure 5a,b presents the lateral-directional aerodynamic characteristics of the twin-engine propeller-driven small aircraft using MAPLA and the available wind tunnel test results [44,45]. Compared to the longitudinal characteristics, the lateral-directional results show higher deviation from the wind tunnel test results. A probable reason for this deviation could be the propeller effects and power contributions. Particularly, these effects change the flow characteristics over the empennage and change its contribution to aerodynamics.

2.2. Initial Design and Sizing of the Propulsion System

To size the electric propulsion system, the number of engines needs to be selected first. Deciding between having fewer larger engines or several smaller engines requires consideration of the following factors [46]:

• Disk loading: bigger propeller disk is better in hover due to lower disk loading.
• Ground clearances: smaller propellers would result in more ground clearance when the horizontal tail and the wing are in a horizontal position.
• Propeller-wing interaction: larger propellers would have a higher slipstream height, which results in more lift because of the propeller slipstream in the same slipstream speed.
• Propeller-propeller interaction: smaller propellers with lower slipstream heights allow easier placement of the engines such that the horizontal tail propellers do not lose their slipstream. This slipstream loss can lead to a reduction in thrust from the horizontal tail propellers and an increase in noise emissions.
• Blade rotation mechanism: using very small propellers makes it more difficult to alter the pitch of the blades.
• Safety in OEI conditions: using more and smaller engines, a failure of one would result in less effect on controllability in hover and on a loss of thrust.

The aircraft features two sets of electric engines with variable-pitch and variable-speed open rotors with six blades, and it was decided to use the same size to minimize differences in propeller loading and provide cheaper manufacturing. A great proportion of the engines are placed on the wing as it is closer to the CG location and requires higher thrust during hover. After considering the criteria for the number of engines and using MAPLA’s electric propulsion module, the desired number of engines was found to be 12 (8 for the wing and 4 for the horizontal tail) with a radius of 0.5 m and a maximum allowable rpm level of ~4800. For controllability purposes, the maximum thrust-to-weight ratio was found to be approximately 1.5. For an MTOW of 2500 kg and 12 propellers, this results in a required maximum thrust of about 3100 N per propeller with a pitch down of 45 deg. The design point for the flight condition was set to the cruise speed of 70 m/s at a 1 km altitude.

The solid-state lithium battery with a specific energy of 500 Wh/kg, a volumetric energy density of 1000 Wh/L, and a power density of 6500 W/kg was chosen. Each battery cell has an internal capacitance of 5 Ah and a nominal voltage of 3.7 V. Future-oriented battery specifications were chosen for the eVTOL tilt-wing aircraft [46]. The regular mission profile as shown in Figure 6 was defined according to the VTOL aircraft flight requirements including 1. Vertical take-off; 2. The transition from hover to climb; 3. Climb to cruising altitude; 4. Cruise; 5. Loiter in horizontal flight; 6. Descent; 7. Transition to hover; 8. Loiter in hover; and 9. Landing [46].

To calculate the energy used during flight using the method presented in [46], thrust at different flight regimes was calculated and converted to power. Accordingly, for take-off and landing, the following relation was used:

$${\mathrm{P}}_{\mathrm{T}\mathrm{L}}=\mathrm{T}{\mathrm{V}}_{\mathrm{\infty }}+\mathrm{k}\mathrm{T}\left(\frac{{-\mathrm{V}}_{\mathrm{\infty }}^2}{2}+\sqrt{\frac{\mathrm{T}}{2\mathsf{\rho }{\mathrm{A}}_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k}}}}\right)$$

(1)

where T accounts for the thrust, V_∞ is the speed perpendicular to the propeller disk, factor k corrects for the power losses, equal to 1.2 in this case, and A is the disc area.

To calculate the required thrust during take-off and landing the following equations of motion were used [46].

$$\sum {\mathrm{F}}_{\mathrm{X}}:-\mathrm{D}\cos \mathsf{\gamma }-\mathrm{L}\sin \mathsf{\gamma }+\mathrm{T}\cos {\mathsf{\theta }}_{\mathrm{T}}=\mathrm{m}{\mathrm{a}}_{\mathrm{x}}$$

(2)

$$\sum {\mathrm{F}}_{\mathrm{Y}}:-\mathrm{D}\sin \mathsf{\gamma }+\mathrm{L}\cos \mathsf{\gamma }+\mathrm{T}\sin {\mathsf{\theta }}_{\mathrm{T}}=\mathrm{m}{\mathrm{a}}_{\mathrm{y}}$$

(3)

where ${\mathsf{\theta }}_{\mathrm{T}}$ represents the rotated wing incidence with respect to the horizontal axis, x-axis, and ${\mathrm{a}}_{\mathrm{x}}$ and ${\mathrm{a}}_{\mathrm{y}}$ are the corresponding accelerations in the x-axis and y-axis, respectively. To find the fuselage drag contribution, it was assumed that the fuselage is at 90 degrees when flying vertically.

For cruise flight conditions, the suggested relation from [46] for power estimation was used

$${\mathrm{P}}_{\mathrm{c}\mathrm{r}\mathrm{u}\mathrm{i}\mathrm{s}\mathrm{e}}=\frac{{\mathrm{C}}_{\mathrm{D}}}{{\mathrm{C}}_{\mathrm{L}}}\frac{\mathrm{W}\mathrm{V}}{\mathsf{\eta }}$$

(4)

where η is the propulsive efficiency at the given speed for cruise flight conditions. Finally, the required energy to be stored by batteries for the entire mission profile can be calculated by [46,47]

$$\mathrm{E}=\sum \mathrm{P}.\mathrm{t}$$

(5)

Consequently, the energy required for a 300 km range flight with 15 min loiter (~86 min) using equations 1 through 5 was estimated to be ~1400 MJ with a maximum power requirement during take-off of ~2.1 MW.

To calculate the battery weight, the following equations according to the energy demands and maximum power were used [46]

$${\mathrm{m}}_{\mathrm{b}\mathrm{a}\mathrm{t}}=\max \left(\frac{{\mathrm{E}}_{\mathrm{r}\mathrm{e}\mathrm{q}}}{{\mathrm{E}}_{\mathrm{s}\mathrm{p}}}\frac{1}{\mathrm{D}\mathrm{o}\mathrm{D}\times \mathrm{E}\mathrm{O}\mathrm{L}\mathrm{C}},\frac{{\mathrm{P}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{\mathrm{P}}_{\mathrm{d}\mathrm{e}\mathrm{n}}}\frac{1}{\mathrm{D}\mathrm{o}\mathrm{D}\times \mathrm{E}\mathrm{O}\mathrm{L}\mathrm{C}}\right)$$

(6)

where DoD is the depth of discharge and EOLC is the end-of-life capacity, and they are both varying between 0 and 1. And the volume of the battery was estimated using [46]:

$${\mathrm{v}}_{\mathrm{b}\mathrm{a}\mathrm{t}}={\mathrm{m}}_{\mathrm{b}\mathrm{a}\mathrm{t}}\times \frac{{\mathrm{E}}_{\mathrm{s}\mathrm{p}}}{{\mathrm{E}}_{\mathrm{v}\mathrm{o}\mathrm{l}}}$$

(7)

where ${\mathrm{E}}_{\mathrm{s}\mathrm{p}}$ is the specific energy in Wh/kg and ${\mathrm{E}}_{\mathrm{v}\mathrm{o}\mathrm{l}}$ stands for the volumetric energy density in Wh/m³ [46]. Consequently, using MAPLA’s electric propulsion module with the total energy and maximum power requirements, the total battery mass for energy storage was estimated to be 776 kg with a required volume of 0.4 m³.

In the following, the total voltage of batteries and the total current, assuming they are connected in series are presented [46].

$${\mathrm{V}}_{\mathrm{t}\mathrm{o}\mathrm{t}}={\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{V}}_{\mathrm{i}}$$

(8)

$$\mathrm{E}={\mathrm{C}\mathrm{V}}_{\mathrm{n}}$$

where $\mathrm{C}$ is the capacity and ${\mathrm{V}}_{\mathrm{n}}$ stands for nominal cell voltage. Also, the required number of battery cells for the electric engines obtained from [46]

$${{\mathrm{n}}_{\mathrm{c}}}_{\mathrm{m}\mathrm{o}\mathrm{t}}=\left[\frac{{\mathrm{E}}_{\mathrm{t}\mathrm{o}\mathrm{t}}\times \mathrm{\%}\mathrm{m}\mathrm{o}\mathrm{t}}{100{\mathrm{E}}_{\mathrm{c}}}\right]$$

(9)

$${{\mathrm{n}}_{\mathrm{c}}}_{\mathrm{m}\mathrm{i}\mathrm{s}\mathrm{c}}=\left[\frac{{\mathrm{E}}_{\mathrm{t}\mathrm{o}\mathrm{t}}\times (100-\mathrm{\%}\mathrm{m}\mathrm{o}\mathrm{t})}{100{\mathrm{E}}_{\mathrm{c}}}\right]$$

where ${\mathrm{n}}_{\mathrm{c}}$ describes the number of battery cells, $\mathrm{\%}\mathrm{m}\mathrm{o}\mathrm{t}$ stands for the energy percentage that goes into the motors and ${\mathrm{E}}_{\mathrm{c}}$ is the energy stored in each individual cell. However, the relation requires additional corrections as a battery package not only needs to provide energy but also requires the correct voltage. Accordingly, the number of batteries required could be calculated using the following relations considering parallel and series connections [46].

$${\mathrm{n}}_{{\mathrm{c}}_{\mathrm{s}\mathrm{e}\mathrm{r}}}=\left[\frac{{\mathrm{V}}_{\mathrm{m}\mathrm{o}\mathrm{t}}}{{\mathrm{V}}_{\mathrm{c}}}\right]$$

(10)

$${\mathrm{n}}_{{\mathrm{c}}_{\mathrm{p}\mathrm{a}\mathrm{r}}}=\left[\frac{{{\mathrm{n}}_{\mathrm{c}}}_{\mathrm{m}\mathrm{o}\mathrm{t}}}{{\mathrm{n}}_{{\mathrm{c}}_{\mathrm{s}\mathrm{e}\mathrm{r}}}}\right]$$

$${\mathrm{n}}_{{\mathrm{c}}_{\mathrm{n}\mathrm{e}\mathrm{w}}}={\mathrm{n}}_{{\mathrm{c}}_{\mathrm{p}\mathrm{a}\mathrm{r}}}\times {\mathrm{n}}_{{\mathrm{c}}_{\mathrm{s}\mathrm{e}\mathrm{r}}}$$

The total number of battery cells for the electric propulsion system and other electronic devices was estimated to be 20,763. Based on the required voltage of 500 V for motors taken from Calnetix [46], using MAPLA’s propulsion module, the number of cells in series is 135 and the required number of battery cells in parallel arrangement is 154 which would result in 20,790 cells due to rounding. As for the proposed eVTOL aircraft, 12 propellers were selected and 2 batteries were required for each engine, giving a total of 24 batteries. Hence, the total number of battery cells in parallel should increase to 164 to allow the modular set of 7 cells in parallel for 24 batteries. Consequently, the total number of cells will be 22,140 with a total mass of 827 kg and a volume of about 0.43 m³.

2.3. PBDO Method Outline

A deterministic gradient optimizer was coupled to a possibility-evaluation algorithm to carry out the design optimization under uncertainty. This gives a comparable framework to the Sequential Optimization and Reliability Assessment (SORA) method developed for use with PBDO [48,49]. Here, the Performance Measure Approach (PMA) was used to implement the possibility evaluation on each constraint. This algorithm employs a separate optimization loop to find the worst-case scenario for each individual constraint. Using this, the value of each uncertain variable is determined within its possibilistic interval and adjusted to ensure that even in extreme cases the constraints will not be violated by an implementation of the output design [37,50]. After this step, for each iteration, a deterministic optimization followed to progress the design. This method of interwoven optimization and possibility evaluation, here achieved by applying Sequential Quadratic Programming (SQP) for each, allows the design to be iteratively altered until a final configuration is selected which is guaranteed to adhere to the constraints [37].

The MDF approach is a single-level MDO procedure where the shared design and coupling variables are consistent across all disciplines, achieved by using an MDA loop. The MDA loop starts with the given design variable states, x and z, introduced by the optimizer using an initial estimation of the coupling variable vector, y. Each discipline is then computed sequentially to update its coupling variable using values from the output of the previous discipline. The loop continues until the coupling variable vector converges. The mathematical representation of a deterministic MDF approach for an optimization problem is presented as follows:

$$\underset{\mathrm{z},\mathrm{x}}{\min }\mathrm{f}\left[\mathrm{z},\mathrm{y}\left(\mathrm{x},\mathrm{y},\mathrm{z}\right),\mathrm{x}\right] \mathrm{i},\mathrm{j}=1,\dots ,\mathrm{n}$$

(11)

$$\mathrm{s}.\mathrm{t}.\mathrm{g}\left[\mathrm{z},\mathrm{y}\left(\mathrm{x},\mathrm{y},\mathrm{z}\right),\mathrm{x}\right]\le 0$$

where the shared design variables are defined by z and the local design variables are defined by x. The coupling variables, y, are found using an MDA loop each time the system and the local variables are updated.

The mathematical representation of a reliability-based MDF approach that was used in this study is shown as follows:

$$\underset{\mathrm{z},\mathrm{x}}{\min }\mathrm{f}\left[\overline{\mathrm{z}},\mathrm{y}\left(\overline{\mathrm{x}},\mathrm{y},\overline{\mathrm{z}}\right),\overline{\mathrm{x}},\overline{\mathrm{p}}\right]$$

(12)

$$\mathrm{s}.\mathrm{t}.\mathrm{P}\left(\mathrm{g}\left[\mathrm{z},\mathrm{y}\left(\mathrm{x},\mathrm{y},\mathrm{z}\right),\mathrm{x}\right]\le 0\right)\ge {\mathrm{P}}_{\mathrm{g}\mathrm{o}\mathrm{l}}$$

where P is the probability of feasibility for each constraint and p is the uncertain parameters. The optimizer alters z and x. The coupling variables, y, are estimated internally by the MDA loop and are not in the system optimizer.

The IDF method was used to reduce the iterations of the computationally costly MDA loop by eliminating the feasibility for every design assessed by the optimizer. Design feasibility, or rather the consistency of the coupling and shared variables, is only assured at the end. New auxiliary constraints are defined to reduce the discrepancy between the discipline responses and the estimated quantities at the optimum point. The mathematical representation of a deterministic IDF approach for an optimization case is presented as follows:

$$\underset{\mathrm{z},\mathrm{x}}{\min }\mathrm{f}\left[\mathrm{z},\mathrm{y}\left(\mathrm{x},{\mathrm{y}}^{\prime },\mathrm{z}\right),\mathrm{x}\right] \mathrm{i},\mathrm{j}=1,\dots ,\mathrm{n}$$

(13)

$$\mathrm{s}.\mathrm{t}.\mathrm{g}\left[\mathrm{z},\mathrm{y}\left(\mathrm{x},{\mathrm{y}}^{\prime },\mathrm{z}\right),\mathrm{x}\right]\le 0$$

$${\mathrm{y}}^{\prime }-\mathrm{y}\left(\mathrm{x},{\mathrm{y}}^{\prime },\mathrm{z}\right)=0$$

where the coupling variables, $\mathrm{y}$, are determined by calculating each discipline’s result using the coupling variable states, ${\mathrm{y}}_0$, defined by the system optimizer. The coupling variables, $\mathrm{y}$, are defined for the system optimization variables in the IDF-based PBDO method considering the auxiliary constraints. With respect to the PMA-based reliability analysis,

$$\min {\mathrm{G}}_{\mathrm{i}}\left(\mathrm{U}\right)$$

(14)

$$\mathrm{s}.\mathrm{t}.‖\mathrm{U}‖={\mathsf{\beta }}_{\mathrm{t}}$$

$$\begin{gathered} \mathrm{becomes} \\ \min {\mathrm{G}}_{\mathrm{i}}(\mathrm{U},{\mathrm{y}}^{\prime }) \end{gathered}$$

$$\mathrm{s}.\mathrm{t}.‖\mathrm{U}‖={\mathsf{\beta }}_{\mathrm{t}}$$

$${\mathrm{y}}^{\prime }-\mathrm{y}\left(\mathrm{x},{\mathrm{y}}^{\prime },\mathrm{z}\right)$$

IDF approaches are also achievable using a revised reliability analysis loop with a double loop or sequential PBDO methods.

Validation of the implemented framework was completed previously [43]. Accordingly, the following nonlinear single-discipline PBDO test example was solved with various methods [51,52]:

$$\min \mathrm{f}={\overline{\mathrm{x}}}_1+{\overline{\mathrm{x}}}_2 \mathrm{s}\mathrm{t}.\mathrm{P}\left[{\mathrm{G}}_{\mathrm{j}}\left(\mathrm{x}\right)\ge 0\right]\ge {\mathrm{P}}_{\mathrm{g}\mathrm{o}\mathrm{a}\mathrm{l},\mathrm{j}},\mathrm{j}=\mathrm{1,2},3$$

(15)

$${\mathrm{G}}_1\left(\mathrm{x}\right)=\frac{{\mathrm{x}}_1^2{\mathrm{x}}_2}{20}-1$$

$${\mathrm{G}}_2\left(\mathrm{x}\right)=\frac{{\left({\mathrm{x}}_1+{\mathrm{x}}_2-5\right)}^2}{30}+\frac{{\left({\mathrm{x}}_1-{\mathrm{x}}_2-12\right)}^2}{120}-1$$

$${\mathrm{G}}_3\left(\mathrm{x}\right)=\frac{80}{\left({\mathrm{x}}_1^2+8{\mathrm{x}}_2+5\right)}-1$$

$${\mathrm{x}}_1~\mathrm{N}\left({\mathrm{x}}_1,0.3\right),{\mathrm{x}}_2~\mathrm{N}({\mathrm{x}}_2,0.3)$$

Table 2. Performance comparison of different optimization methods [43].

Method	${\mathbf{f}}_{\mathbf{m}\mathbf{i}\mathbf{n}}$	${\mathbf{X}}_1$	${\mathbf{X}}_2$	${\mathbf{n}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l}}$
Deterministic	5.1769	3.1134	2.0636	16
Single Loop	6.6198	3.4413	3.2866	16
PMA/Sequential	6.7043	3.4506	3.2537	651
PMA/Double Loop	6.7043	3.4506	3.2537	1004
RIA/Double Loop	6.7257	3.4391	3.2866	1530

shows the results of the comparison between various methods including the double-loop PMA approach, the double-loop Reliability Index Approach (RIA), and the sequential PMA method.

Of the tested configurations, the single-loop approach diverged from PMA and RIA methods for large values of ${\mathsf{\beta }}_{\mathrm{t}}$ and as such, was not selected. Overall, the PMA method was found to require fewer iterations when compared to its RIA counterparts; with the sequential PMA approach requiring fewer iterations than its double-loop implementation. For this reason, it was the approach that was ultimately selected, with the hope that fewer iterations for this problem would translate to a more rapid aircraft optimization process.

The speed and stability of this method were also assessed by running the optimization procedure through 100 starting values. It was found to, on average, have an error below ${10}^4$ %, and require less time than the double loop methods [43].

For further validation of the MDO methodology, the following multi-discipline problem was solved and compared with different methods [43,53].

$$\min \mathrm{f}=-{\left({\mathrm{x}}_1-6\right)}^3+{\mathrm{y}}_1^2-\exp \left(-\frac{{\mathrm{y}}_1}{{\mathrm{y}}_2}\right)$$

(16)

$${\mathrm{y}}_1={\mathrm{x}}_1^2+\frac{{\mathrm{y}}_2}{2}$$

$${\mathrm{g}}_1=-{\mathrm{y}}_2+\exp \left(\frac{{\mathrm{y}}_1}{{\mathrm{y}}_2}+2.2{\mathrm{x}}_1\right)$$

$${\mathrm{y}}_2={\mathrm{x}}_1+{\mathrm{x}}_2+\frac{3{\mathrm{x}}_1{\mathrm{x}}_2}{{\mathrm{y}}_1}$$

$${\mathrm{g}}_2={\mathrm{y}}_2-{\mathrm{y}}_1-{\left({\mathrm{x}}_1+1\right)}^2-{\left({\mathrm{x}}_2-4\right)}^3$$

where the uncertain variables of ${\mathrm{x}}_1$ and ${\mathrm{x}}_2$ were considered to have a variation coefficient of 0.04 defined as the standard deviation, $\sigma$, to the random variable mean value, $\mu$. After solving the problem using different methods at a reliability level of $3\sigma$ and finding the optimum points, the results are presented along with the number of function evaluations, ${\mathrm{n}}_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l}\mathrm{s}}$ in

Table 3. Comparison of different methods’ stability [43].

Method	${\mathbf{f}}_{\mathbf{m}\mathbf{i}\mathbf{n}}$	${\mathbf{X}}_1$	${\mathbf{X}}_2$	${\mathbf{n}}_{\mathbf{e}\mathbf{v}\mathbf{a}\mathbf{l}}$
MDF/Deterministic	115.797	1.648	3.004	68
MDF/Single Loop	119.942	1.656	2.913	1495
MDF/Double Loop/PMA	119.656	1.598	2.806	16,474
MDF/Sequential/PMA	119.683	1.621	2.843	7413
IDF/Double Loop/PMA	119.446	1.600	2.787	23,765

.

The results indicated that the IDF method with a sequential PMA strategy for the PBDO method achieved the most efficient results. Although MDF-based methods were also accurate, the effort needed to solve the problem and find the most optimized case with regard to the ${\mathrm{n}}_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l}\mathrm{s}}$ was considered to be excessive.

3. MDO Framework

3.1. Methods

The optimization framework is similar to the one used for Example 1 and Example 2 and the IDF/Sequential/PMA method with an SQP optimization procedure was implemented. The algorithm consists of a PBDO solver coupled with MAPLA’s different disciplines including aerodynamics, mass and balance, electric propulsion, along with stability and control.

Considering various disciplines of aerodynamics, mass and balance, electric propulsion and battery management, as well as stability and control, the XDSM block diagram in Figure 7 is described for the MDO of the eVTOL tilt-wing aircraft. Accordingly, to run the system-level optimizer, each sublevel was estimated sequentially with the result of one discipline used by another. Accordingly, the empennage and wing surface areas were optimized to achieve minimum drag based on the shape parameters presented and described in section III.B and the corresponding flight environment. The framework is described in the following and the corresponding XDSM flowchart [54] is presented in Figure 8. The steps for the approach are as follows:

• Determine the estimated errors

  (I)

  Model each entry based on the horizontal/vertical tails, wing model and battery management

  (II)

  Calculate the aerodynamics of the eVTOL tilt-wing aircraft

  (III)

  Apply the uncertainty analysis

  (IV)

  Use the best-fit PDF curve for each uncertainty

• Select the initial vector and a set of reliability indices
• Run the IDF-based deterministic optimization
• Run the MDF-based PMA reliability assessment at the current reliability index, and change variables based on the sequential method
• Check for convergence based on the current reliability goal

  (I)

  if yes, update starting vector with the current optimum point, select the next reliability goal and return to Step 3

  (II)

  if no, return to Step 3

• Advance to the next target reliability level

  (I)

  retain solution as a new starting vector

  (II)

  return to Step 3

3.2. Design Variables

Since, in this work, the focus was to optimize the wing and empennage of the eVTOL aircraft, the corresponding variables were chosen based on the semi-empirical analysis’s sensitivities considering the effectiveness of each component with respect to the semi-empirical analysis sensitivity. These variables are presented in Figure 8, Figure 9 and Figure 10. In addition to the horizontal/vertical tail and wing components, battery placement was also parametrized because of its effect on the desired center of gravity location to minimize manufacturing penalties. A list of the design variables considered in this study is supplied in

Table 4. Design variables and their descriptions for the eVTOL tilt-wing Aircraft.

Component	Variable	Details	Limits
Horizontal Tail	${\mathrm{i}}_{\mathrm{h}}$	Incidence angle of horizontal tail, deg	0 to 3
	bh	Span of horizontal tail, m	4 to 6
	${\mathrm{c}}_{{\mathrm{r}}_{\mathrm{h}}}$	Root chord of horizontal tail, m	1 to 1.45
	${\mathrm{c}}_{{\mathrm{t}}_{\mathrm{h}}}$	Tip chord of horizontal tail, m	0.5 to 1
	${{\mathsf{\Lambda }}_{\mathrm{L}\mathrm{E}}}_{\mathrm{h}}$	Horizontal tail leading edge sweep angle, deg	10 to 20
	lh	Distance, parallel to X-axis, from the the horizontal tail mean aerodynamic chord to the nose of the fuselage, m	8 to 8.7
	zh	Distance, parallel to Z-axis, from the quarter chord of the horizontal tail mean aerodynamic chord to the X-axis, positive down, m	−0.4 to −0.1
	celevator	Elevator to horizontal tail chord ratio	0.2 to 0.5
Vertical Tail	bv	Span of vertical tail, m	1.8 to 2.2
	${\mathrm{c}}_{{\mathrm{r}}_{\mathrm{v}}}$	Root chord of vertical tail, m	1.9 to 2.5
	${\mathrm{c}}_{{\mathrm{t}}_{\mathrm{v}}}$	Tip chord of vertical tail, m	0.8 to 1.4
	${\mathsf{\varphi }}_{\mathrm{T}\mathrm{E}}$	Vertical tail Trailing edge sweep angle, deg	10 to 20
	zv	Perpendicular distance from root chord of vertical-tail to X-axis, positive down, m	−0.35 to −0.15
	lv	Distance along X-axis from the leading edge of tip chord of vertical tail to the nose of the fuselage, m	8.5 to 9.5
	crudder	Rudder to vertical tail chord ratio, m	0.2 to 0.5
Wing	${\mathrm{i}}_{\mathrm{w}}$	Incidence angle of the wing, deg	2 to 5
	${\mathsf{\alpha }}_{\mathrm{t}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{t}}$	Incidence angle of the wing	−2 to −5
	bw	Wing span, m	10 to 12
	${\mathrm{c}}_{{\mathrm{r}}_{\mathrm{w}}}$	Root chord of the wing, m	1.8 to 2.5
	${\mathrm{c}}_{{\mathrm{t}}_{\mathrm{w}}}$	Tip chord of the wing, m	0.8 to 1.1
	${{\mathsf{\Lambda }}_{\mathrm{L}\mathrm{E}}}_{\mathrm{w}}$	Leading edge sweep angle of the wing, deg	2 to 5
	${{\mathsf{\Lambda }}_{\mathrm{T}\mathrm{E}}}_{\mathrm{w}}$	Trailing edge sweep angle of the wing, deg	−6 to −12
	lw	Distance, parallel to X-axis, from the leading edge of wing mean aerodynamic chord to the nose of fuselage, m	2.2 to 4
	zw	Distance, parallel to Z-axis, from the quarter chord of the wing mean aerodynamic chord to the X-axis, positive down, m	0.2 to 0.6
	$\mathsf{\Gamma }$	Dihedral angle, deg	5 to 10
Battery Placement	${\mathrm{l}}_{{\mathrm{v}}_{{\mathrm{b}}_1}}$	Battery set 1, distance from nose, m	0.8 to 1.3
	${\mathrm{l}}_{{\mathrm{v}}_{{\mathrm{b}}_2}}$	Battery set 2, distance from nose, m	6.4 to 6.9
	${\mathrm{N}}_{{\mathrm{b}}_1}$	Battery set 1, number of batteries	1 to 24
	${\mathrm{N}}_{{\mathrm{b}}_2}$	Battery set 2, number of batteries	1 to 24

.

The wing surface area was determined using Equation (17) considering the shape parameters depicted in Figure 8:

$${\mathrm{S}}_{\mathrm{w}}={\mathrm{c}}_{{\mathrm{r}}_{\mathrm{w}}}{\mathrm{b}}_{\mathrm{w}}-{\mathrm{c}}_{{\mathrm{e}}_{{\mathrm{w}}_{\mathrm{f}\mathrm{o}\mathrm{r}}}}\frac{{\mathrm{b}}_{\mathrm{w}}}{2}-{\mathrm{c}}_{{\mathrm{e}}_{{\mathrm{w}}_{\mathrm{b}\mathrm{a}\mathrm{c}\mathrm{k}}}}\frac{{\mathrm{b}}_{\mathrm{w}}}{2}$$

(17)

where ${\mathrm{c}}_{{\mathrm{e}}_{{\mathrm{w}}_{\mathrm{f}\mathrm{o}\mathrm{r}}}}$ is the extended panel out of the wing tip chord and is achieved in Figure 8 by the trailing edge angle of the wing.

$${\mathrm{c}}_{{\mathrm{e}}_{{\mathrm{w}}_{\mathrm{f}\mathrm{o}\mathrm{r}}}}=\frac{{\mathrm{b}}_{\mathrm{w}}}{2}\tan \left({\mathsf{\Lambda }}_{\mathrm{L}{\mathrm{E}}_{\mathrm{w}}}\right)$$

(18)

And ${\mathrm{c}}_{{\mathrm{e}}_{{\mathrm{w}}_{\mathrm{b}\mathrm{a}\mathrm{c}\mathrm{k}}}}$ can be calculated using the

$${\mathrm{c}}_{{\mathrm{e}}_{{\mathrm{w}}_{\mathrm{b}\mathrm{a}\mathrm{c}\mathrm{k}}}}=\frac{{\mathrm{b}}_{\mathrm{w}}}{2}\tan \left({\mathsf{\Lambda }}_{\mathrm{T}{\mathrm{E}}_{\mathrm{w}}}\right)$$

(19)

where ${\mathsf{\Lambda }}_{\mathrm{T}{\mathrm{E}}_{\mathrm{w}}}$ is the wing trailing edge sweep angle and can be achieved by Equation (19) as follows.

$${\mathsf{\Lambda }}_{\mathrm{T}{\mathrm{E}}_{\mathrm{w}}}=90-\arctan \left(\frac{{\mathrm{b}}_{\mathrm{w}}}{2\left({\mathrm{c}}_{{\mathrm{r}}_{\mathrm{w}}}-{\mathrm{c}}_{{\mathrm{t}}_{\mathrm{w}}}-{\mathrm{c}}_{{\mathrm{e}}_{{\mathrm{w}}_{\mathrm{f}\mathrm{o}\mathrm{r}}}}\right)}\right)$$

(20)

The surface area of the horizontal tail was estimated using Equation (21) considering the shape parameters presented in Figure 9 as follows.

$${\mathrm{S}}_{\mathrm{h}}={\mathrm{c}}_{{\mathrm{r}}_{\mathrm{h}}}{\mathrm{b}}_{\mathrm{h}}-{\mathrm{c}}_{{\mathrm{e}}_{{\mathrm{h}}_{\mathrm{f}\mathrm{o}\mathrm{r}}}}\frac{{\mathrm{b}}_{\mathrm{h}}}{2}-{\mathrm{c}}_{{\mathrm{e}}_{{\mathrm{h}}_{\mathrm{b}\mathrm{a}\mathrm{c}\mathrm{k}}}}\frac{{\mathrm{b}}_{\mathrm{h}}}{2}$$

(21)

where ${\mathrm{c}}_{{\mathrm{e}}_{{\mathrm{h}}_{\mathrm{f}\mathrm{o}\mathrm{r}}}}$ is the extended panel out of the horizontal tail tip chord and can be seen in Figure 9 by the horizontal tail trailing edge angle as follows.

$${\mathrm{c}}_{{\mathrm{e}}_{{\mathrm{h}}_{\mathrm{f}\mathrm{o}\mathrm{r}}}}=\frac{{\mathrm{b}}_{\mathrm{h}}}{2}\tan \left({\mathsf{\Lambda }}_{\mathrm{L}{\mathrm{E}}_{\mathrm{h}}}\right)$$

(22)

And ${\mathrm{c}}_{{\mathrm{e}}_{{\mathrm{h}}_{\mathrm{b}\mathrm{a}\mathrm{c}\mathrm{k}}}}$ can be calculated using the

$${\mathrm{c}}_{{\mathrm{e}}_{{\mathrm{h}}_{\mathrm{b}\mathrm{a}\mathrm{c}\mathrm{k}}}}=\frac{{\mathrm{b}}_{\mathrm{h}}}{2}\tan \left({\mathsf{\Lambda }}_{\mathrm{T}{\mathrm{E}}_{\mathrm{h}}}\right)$$

(23)

where ${\mathsf{\Lambda }}_{\mathrm{T}{\mathrm{E}}_{\mathrm{h}}}$ is the horizontal tail trailing edge sweep angle and can be achieved by Equation (24) as follows.

$${\mathsf{\Lambda }}_{\mathrm{T}{\mathrm{E}}_{\mathrm{h}}}=90-\arctan \left(\frac{{\mathrm{b}}_{\mathrm{h}}}{2\left({\mathrm{c}}_{{\mathrm{r}}_{\mathrm{h}}}-{\mathrm{c}}_{{\mathrm{t}}_{\mathrm{h}}}-{\mathrm{c}}_{{\mathrm{e}}_{{\mathrm{h}}_{\mathrm{f}\mathrm{o}\mathrm{r}}}}\right)}\right)$$

(24)

Similarly, the surface area of the vertical tail was determined using Equation (25), based on the shape parameters defined in Figure 9:

$${\mathrm{S}}_{\mathrm{v}}=\left(\left({\mathrm{c}}_{{\mathrm{r}}_{\mathrm{v}}}-{\mathrm{c}}_{{\mathrm{t}}_{\mathrm{v}}}+{\mathrm{c}}_{{\mathrm{e}}_{\mathrm{v}}}\right)\frac{{\mathrm{b}}_{\mathrm{v}}}{2}\right)+\left({\mathrm{c}}_{{\mathrm{t}}_{\mathrm{v}}}{\mathrm{b}}_{\mathrm{v}}\right)\left(\frac{{\mathrm{b}}_{\mathrm{v}}^2\mathrm{t}\mathrm{a}\mathrm{n}\left({\mathsf{\varphi }}_{\mathrm{T}\mathrm{E}}\right)}{2}\right)$$

(25)

where ${\mathrm{c}}_{{\mathrm{e}}_{\mathrm{v}}}$ is the extended panel out of the vertical tail root chord and could be achieved using the trailing edge sweep angle with the following equation.

$${\mathrm{c}}_{{\mathrm{e}}_{\mathrm{v}}}={\mathrm{b}}_{\mathrm{v}}\mathrm{t}\mathrm{a}\mathrm{n}\left({\mathsf{\varphi }}_{\mathrm{T}\mathrm{E}}\right)$$

(26)

3.3. Disciplines

3.3.1. Aerodynamics and Stability

The first discipline in this study was aerodynamics consisting of longitudinal and lateral-directional aerodynamic coefficients and stability derivatives. MAPLA’s aerodynamics module was used and was originally developed based on the work accomplished by NASA [44,45]; however, to generalize the platform and enhance the work, particularly for the purpose of small aircraft design and development, some changes were made to include the following features [39]:

• The maximum lift of a twisted wing
• The zero-lift pitching moment of the twisted wing
• The drag of the twisted wing
• The high-lift surfaces
• Lift of the horizontal tail and elevator surfaces
• The drag of the control and high-lift surfaces

The aerodynamics module results were compared against the wind tunnel test data in section II.A. In addition, the aerodynamics module results were compared with wind tunnel test data of a small aircraft at a scale of 1:3, the estimations provided by the DATCOM, and a VLM-based prediction in our previous work [39]. The results indicated that the aerodynamics module was able to determine the aerodynamic characteristics of the small aircraft with acceptable accuracy from zero-lift to stall conditions in all configurations. For the purpose of this study, aerodynamic characteristics were first estimated for each component of the tailless eVTOL tilt-wing aircraft, while running the optimization system, the horizontal/vertical tail contribution was also added to the calculations to consider the contributions of the fixed part of the empennage and high-lift surfaces. In the next step, the static stability derivatives were determined, and the characteristics were found by integrating the contributions of each component of the eVTOL tilt-wing aircraft. After that, the dynamic stability of the eVTOL tilt-wing aircraft was calculated with static stability characteristics.

3.3.2. Weight and Balance

The next discipline was weight and balance where the weight of each part was estimated using statistical estimations [55,56]. The weights were determined in accordance with the MIL-STD-1374 through a “Summary Group Weight Statement” for light aircraft. Both the mass and center of gravity of each part were estimated using statistical relations to provide the total empty weight of the eVTOL tilt-wing aircraft [57,58]. Finally, the moments of inertia calculations were carried out for ${\mathrm{I}}_{\mathrm{x}\mathrm{x}}$, ${\mathrm{I}}_{\mathrm{y}\mathrm{y}}$ and ${\mathrm{I}}_{\mathrm{z}\mathrm{z}}$ using the following relations [55]:

$${\mathrm{I}}_{\mathrm{x}\mathrm{x}}=\frac{{\left({\mathrm{b}}^2·\mathrm{W}·\mathrm{R}\left(1\right)\right)}^2}{4}$$

(27)

$${\mathrm{I}}_{\mathrm{y}\mathrm{y}}=\frac{{\left({\mathrm{l}}_{\mathrm{f}}^2·\mathrm{W}·\mathrm{R}\left(2\right)\right)}^2}{4}$$

(28)

$${\mathrm{I}}_{\mathrm{z}\mathrm{z}}=\frac{{\left({(\frac{\mathrm{b}}{2}+\frac{{\mathrm{l}}_{\mathrm{f}}}{2})}^2·\mathrm{W}·\mathrm{R}\left(3\right)\right)}^2}{4}$$

(29)

$$\mathrm{R}=[0.25,0.38,0.39]$$

(30)

where b is the wingspan in m, W is the total weight in kg and ${\mathrm{l}}_{\mathrm{f}}$ is the length of the fuselage in m.

3.3.3. Handling Quality

The next discipline was control where the handling quality of the aircraft was evaluated using MAPLA’s Flying Quality subprogram. It enables the assessment of the trim characteristics of the aircraft in different flight conditions and aircraft configurations, hence presenting the handling quality level for all longitudinal and lateral-directional modes in accordance with the regulations. By altering mass, the moment of inertia, geometry, aerodynamic coefficients, stability derivatives, the centre of mass, and the handling qualities of the eVTOL tilt-wing airplane were assessed. Hence, having the required aerodynamic characteristics, trim conditions were applied to generate the data for the trimmed flight. Accordingly, elevator deflection, aircraft angle of attack and power values were estimated considering the trim flight as follows [42,43]:

$${\mathrm{C}}_{\mathrm{L}}={\mathrm{C}}_{{\mathrm{L}}_1}=\frac{\mathrm{W}}{0.5·\mathsf{\rho }·{\mathrm{V}}^2·{\mathrm{S}}_{\mathrm{w}}}$$

(31)

$${\mathrm{C}}_{\mathrm{m}}=0$$

(32)

The perturbed equations of motion of the eVTOL tilt-wing aircraft were modelled based on the small perturbation equations. Having the transfer function for each case, the corresponding damping ratio and natural frequency were calculated and stored to be used for the handling quality analysis. Then, using airworthiness and aviation regulation requirements for small aircraft, particularly from CS-23, CS-27, CS-29, or CS-25 and available information from military regulations for aircraft dynamic modes, the aircraft handling quality level for different longitudinal and lateral-directional modes were estimated and ranked from 1 to 4. Level 1 denotes the best flying quality, while Level 4 represents poor flying quality for the aircraft mode [42,43].

3.3.4. Battery System

Finally, the battery system sizing and placement were estimated to find the optimum location scenario to keep the aft and forward CG positions in the desired locations for manufacturing purposes. Hence, considering the number of battery cells as 22,140 in total for a range of 300 km with 15 min loiter (~86 min) based on the results from the initial propulsion system sizing presented in section II.B, the arrangement of the battery package optimized. Accordingly, the battery sets were divided into two compartments where an optimization process was conducted to find the number of cells and the distance of each one of the compartments from the nose of the airplane that could compromise the required number of cells in total while maintaining the most aft and forward CG positions (e.g., 27% and 10%, respectively) in the desired manufacturing position of the base aircraft.

3.4. Objective Function

An objective function to guide the design was required for this investigation, selected to achieve the minimum drag force of the components while satisfying the design constraints of the total aircraft. Accordingly, the Stability and Control module of MAPLA was employed while running the optimization algorithm to find the trim condition in the cruise and the corresponding drag force. The contributions of each component to the total drag of the eVTOL tilt-wing aircraft are summarized here:

• Zero-lift drag contribution of the wing, horizontal tail, and vertical tail.
• Zero-lift drag contribution of the fuselage and nacelles.
• Zero-lift interference drag contribution of the components.
• Drag contribution of the wing and empennage package at the angle of attack.
• Drag contribution of the fuselage and nacelle at the angle of attack.
• Wing-fuselage interference drag contribution at the angle of attack.

Accordingly, the main objective was defined to minimize the total drag of the eVTOL tilt-wing aircraft using Equation (33).

Main objective:

$$\min {\mathrm{C}}_{{\mathrm{D}}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}=\min (\mathrm{C}{\mathrm{D}}_{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{g}}+\mathrm{C}{\mathrm{D}}_{\mathrm{e}\mathrm{m}\mathrm{p}\mathrm{e}\mathrm{n}\mathrm{n}\mathrm{a}\mathrm{g}\mathrm{e}}+{\mathrm{C}\mathrm{D}}_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\mathrm{s}})+(\mathrm{C}{\mathrm{D}}_{\mathrm{f}\mathrm{u}\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{a}\mathrm{g}\mathrm{e}}+\mathrm{C}{\mathrm{D}}_{\mathrm{n}\mathrm{a}\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\mathrm{e}\mathrm{s}})$$

(33)

3.5. Constraints

A series of constraints were imposed on the optimization problem. Starting with the manufacturing constraints were applied to define the maximum and minimum acceptable values for each variable. Considering historical data and initial weight and balance estimation, the empty mass of the eVTOL tilt-wing airplane was restricted to less than 1600 kg. According to the Special Condition for small-category VTOL aircraft by EASA light aircraft certification requirements and CS-27, small rotorcrafts, the weight of the eVTOL must be less than 3175 Kg [59]. Static and dynamic characteristics requirements were also set as constraints for the optimization model. Moreover, the effectiveness of the control surfaces during cruise and climb-out flights was determined to guarantee adequate pitch authority as accomplished in slow flights.

Regulation compliance for VTOL aircraft according to the Means of Compliance (MOC) with the Special Condition VTOL aircraft from EASA [60], should use regulations from CS-23, CS-27, CS-29, or CS-25 unless specifically indicated by the MOC VTOL. Accordingly, for the case of the proposed eVTOL tilt-wing aircraft regulations are considered from CS-23 for light aircraft.

A negative value for ${{\mathrm{C}}_{\mathrm{m}}}_{\mathsf{\alpha }}$ was required to achieve longitudinal static stability. Typically, ${{\mathrm{C}}_{\mathrm{m}}}_{\mathsf{\alpha }}$ for light aircraft should be around −0.3 to −1.5 1/rad. Regarding longitudinal dynamic stability, a negative value for the real part of the roots of the longitudinal characteristics was desired. ${{\mathrm{C}}_{\mathrm{m}}}_{\mathrm{q}}$ was a key player for this criterion such that a negative value has a strong stabilizing outcome. Usually, ${{\mathrm{C}}_{\mathrm{m}}}_{\mathrm{q}}$ for most small aircraft is around −5 to −30 1/rad [61]. Furthermore, a static margin of 5–10% was constrained to longitudinal stability.

A negative value for ${\mathrm{C}}_{{\mathrm{l}}_{\mathsf{\beta }}}$ was necessary to maintain roll stability. In addition, ${\mathrm{C}}_{{\mathrm{n}}_{\mathsf{\beta }}}$ must be positive to have directional static stability. The typical value of ${\mathrm{C}}_{{\mathrm{n}}_{\mathsf{\beta }}}$ for most aircraft is between +0.1 and +0.4 1/rad. A negative value for the real part of the roots of the lateral-directional characteristic equation is required to maintain directional dynamic stability. ${\mathrm{C}}_{{\mathrm{n}}_{\mathrm{r}}}$ is a major contributor to this criterion such that a negative value has a robust stabilizing effect. The typical value for most aircraft is between −0.1 and −1 1/rad [62]. For a full-thrust climb scenario in OEI flight circumstances, yaw static stability must be ensured. The constraint ensures that the rudder and vertical tail have adequate authority [63,64,65,66].

Airworthiness requirements for small airplanes were considered, principally from the CS-23 standard for aircraft dynamic modes and defined as constraints for the optimization model. All constraints are listed in

Table 5. Constraints used for the optimization of the eVTOL tilt-wing aircraft.

Constraint	Description	Limits
${\mathbf{W}}_{\mathbf{e}}$	Empty weight, kg	${\mathrm{W}}_{\mathrm{e}}$ < 1700
${\mathbf{C}}_{{\mathbf{m}}_{\mathsf{\alpha }}}$	Pitching moment coefficient, a/rad	$-1.5<{\mathrm{C}}_{{\mathrm{m}}_{\mathsf{\alpha }}}<-0.3$
${\mathbf{C}}_{{\mathbf{n}}_{\mathsf{\beta }}}$	Weathercock stability coefficient, 1/rad	$0.1<{\mathrm{C}}_{{\mathrm{n}}_{\mathsf{\beta }}}<0.4$
${\mathbf{C}}_{{\mathbf{l}}_{\mathsf{\beta }}}$	Effective dihedral coefficient, 1/rad	${\mathrm{C}}_{{\mathrm{l}}_{\mathsf{\beta }}}<0$
${\mathbf{C}}_{{\mathbf{m}}_{\mathbf{q}}}$	Pitching moment coefficient due to pitch rate, 1/rad	$-30<{\mathrm{C}}_{{\mathrm{m}}_{\mathrm{q}}}<-5$
${\mathbf{C}}_{{\mathbf{n}}_{\mathbf{r}}}$	Damping in yaw derivative, 1/rad	$-1<{\mathrm{C}}_{{\mathrm{n}}_{\mathrm{r}}}<-0.1$
$\mathbf{C}{\mathbf{G}}_{\mathbf{a}\mathbf{f}\mathbf{t}}$	Centre of gravity, maximum afterwards, %	<27%
$\mathbf{C}{\mathbf{G}}_{\mathbf{f}\mathbf{o}\mathbf{r}}$	Centre of gravity, maximum forwards, %	>10%
SM	Static margin	$5\mathrm{\%}\le \mathrm{S}\mathrm{M}\le 10\mathrm{\%}$
${\mathsf{\delta }}_{{\mathbf{r}}_{\mathbf{m}\mathbf{a}\mathbf{x}}}$	Max rudder deflection, deg	$\pm 25$
${\mathsf{\delta }}_{{\mathbf{e}}_{\mathbf{m}\mathbf{a}\mathbf{x}}}$	Max elevator deflection, deg	$\pm 25$
$\mathbf{H}{\mathbf{Q}}_{\mathbf{P}}$	Handling quality, phugoid mode	$\mathrm{H}{\mathrm{Q}}_{\mathrm{P}}=1$
$\mathbf{H}{\mathbf{Q}}_{\mathbf{S}\mathbf{P}}$	Handling quality, short period mode	$\mathrm{H}{\mathrm{Q}}_{\mathrm{S}\mathrm{P}}=1$
$\mathbf{H}{\mathbf{Q}}_{\mathbf{D}}$	Handling quality, Dutch roll mode	$\mathrm{H}{\mathrm{Q}}_{\mathrm{D}}=1$
$\mathbf{H}{\mathbf{Q}}_{\mathbf{R}}$	Handling quality, roll mode	$\mathrm{H}{\mathrm{Q}}_{\mathrm{R}}=1$
$\mathbf{H}{\mathbf{Q}}_{\mathbf{S}}$	Handling quality, spiral mode	$\mathrm{H}{\mathrm{Q}}_{\mathrm{S}}=1$

.

3.6. Source of Uncertainty

In this example, instead of using a uniformly distributed random parameter for the error term of aero uncertainty, a nonuniform approach was used where each aerodynamic characteristic presented in

Table 6. Nonuniform error terms used for the uncertainty of the aerodynamic characteristics [44,45].

Aerodynamic Characteristics	Deviation
$\mathbf{C}{\mathbf{L}}_{\mathsf{\alpha }}$	$5\mathrm{\%}$
$\mathbf{C}{\mathbf{D}}_{\mathsf{\alpha }}$	$10\mathrm{\%}$
$\mathbf{C}{\mathbf{m}}_{\mathsf{\alpha }}$	$15\mathrm{\%}$
$\mathbf{C}{\mathbf{m}}_{\mathbf{q}}$	$30\mathrm{\%}$
$\mathbf{C}{\mathbf{m}}_{\dot{\mathsf{\alpha }}}$	$30\mathrm{\%}$
$\mathbf{C}{\mathbf{L}}_{\mathbf{q}}$	NA, assumed 20%
$\mathbf{C}{\mathbf{L}}_{\dot{\mathsf{\alpha }}}$	NA, assumed 20%
${\mathbf{C}}_{{\mathbf{n}}_{\mathsf{\beta }}}$	$20\mathrm{\%}$
${\mathbf{C}}_{{\mathbf{l}}_{\mathsf{\beta }}}$	$5\mathrm{\%}$
${\mathbf{C}}_{{\mathbf{l}}_{\mathbf{p}}}$	NA, assumed 20%
${\mathbf{C}}_{{\mathbf{n}}_{\mathbf{p}}}$	NA, assumed 20%
${\mathbf{C}}_{{\mathbf{l}}_{\mathbf{r}}}$	$10\mathrm{\%}$
${\mathbf{C}}_{{\mathbf{n}}_{\mathbf{r}}}$	$10\mathrm{\%}$

had a specific error term according to the similar aircraft model used as the base for the design of the new aircraft [44,45]. Accordingly, the error term is considered to be distributed between the maximum and minimum detected value. The error factor was defined as the ratio of estimated aerodynamic results to the observed data from the original database, ${\mathsf{\epsilon }}_{\mathrm{a}\mathrm{e}\mathrm{r}\mathrm{o}}$. The aerodynamics results estimated using MAPLA’s Aerodynamics module were then scaled by this ratio. The error factor determined by Equation (34) was set to be nonuniformly distributed according to the deviation results shown in Table 6.

$${\mathsf{{\rm E}}}_{\mathrm{a}\mathrm{e}\mathrm{r}\mathrm{o}}=\frac{\mathrm{a}\mathrm{e}\mathrm{r}{\mathrm{o}}_{{\mathrm{c}\mathrm{h}\mathrm{a}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{s}}_{\mathrm{p}}}}{\mathrm{a}\mathrm{e}\mathrm{r}{\mathrm{o}}_{\mathrm{c}\mathrm{h}\mathrm{a}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{s}}}$$

(34)

$$\mathrm{a}\mathrm{e}\mathrm{r}{\mathrm{o}}_{\mathrm{c}\mathrm{h}\mathrm{a}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{s}}={\mathsf{\epsilon }}_{\mathrm{a}\mathrm{e}\mathrm{r}\mathrm{o}}\mathrm{a}\mathrm{e}\mathrm{r}{\mathrm{o}}_{{\mathrm{c}\mathrm{h}\mathrm{a}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{s}}_{\mathrm{p}}}$$

The aircraft models from the compiled database were used and the calculated empty mass of each one was compared with the original results from the database. Equation (35) was used to define the mass error factor. Hence, empty-mass calculations were determined by MAPLA first and scaled by the error factor to consider the uncertainty in the estimation of the eVTOL tilt-wing aircraft mass.

$${\mathsf{{\rm E}}}_{{\mathrm{M}}_{\mathrm{e}}}=\frac{{\mathrm{M}}_{{\mathrm{e}}_{\mathrm{p}}}}{{\mathrm{M}}_{{\mathrm{e}}_{\mathrm{d}\mathrm{b}}}}$$

(35)

$${\mathrm{M}}_{\mathrm{e}}=\frac{{\mathrm{M}}_{{\mathrm{e}}_{\mathrm{p}}}}{{\mathsf{\epsilon }}_{{\mathrm{M}}_{\mathrm{e}}}}$$

4. Results

After running the optimization process considering both the constraints and uncertainties, the resulting characteristics of the aircraft are presented in Figure 11, generated by the Geometry module of MAPLA. The resulting variables are discussed in

Table 7. Design variables for optimization of the eVTOL tilt-wing aircraft.

Component	Variable	Description	Result
Horizontal Tail	${\mathrm{i}}_{\mathrm{h}}$	Incidence angle of horizontal tail, deg	1.92
	bh	Span of horizontal tail, m	4.95
	${\mathrm{c}}_{{\mathrm{r}}_{\mathrm{h}}}$	Root chord of horizontal tail, m	1.29
	${\mathrm{c}}_{{\mathrm{t}}_{\mathrm{h}}}$	Tip chord of horizontal tail, m	0.82
	${{\mathsf{\Lambda }}_{\mathrm{L}\mathrm{E}}}_{\mathrm{h}}$	Horizontal tail leading edge sweep angle, deg	12.2
	lh	Distance, parallel to X-axis, from the the horizontal tail mean aerodynamic chord to the nose of the fuselage, m	8.32
	zh	Distance, parallel to Z-axis, from the quarter chord of the horizontal tail mean aerodynamic chord to the X-axis, positive down, m	−0.278
	celevator	Elevator to horizontal tail chord ratio	0.37
Vertical Tail	bv	Span of vertical tail, m	1.847
	${\mathrm{c}}_{{\mathrm{r}}_{\mathrm{v}}}$	Root chord of vertical tail, m	1.955
	${\mathrm{c}}_{{\mathrm{t}}_{\mathrm{v}}}$	Tip chord of vertical tail, m	0.874
	${\mathsf{\varphi }}_{\mathrm{T}\mathrm{E}}$	Vertical tail Trailing edge sweep angle, deg	17.15
	zv	Perpendicular distance from root chord of vertical-tail to X-axis, positive down, m	−0.23
	lv	Distance along X-axis from the leading edge of tip chord of vertical tail to the nose of the fuselage, m	9.08
	crudder	Rudder to vertical tail chord ratio, m	0.418
Wing	${\mathrm{i}}_{\mathrm{w}}$	Incidence angle of the wing, deg	2.74
	${\mathsf{\alpha }}_{\mathrm{t}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{t}}$	Incidence angle of the wing	−3.15
	bw	Wing span, m	11.95
	${\mathrm{c}}_{{\mathrm{r}}_{\mathrm{w}}}$	Root chord of the wing, m	2.143
	${\mathrm{c}}_{{\mathrm{t}}_{\mathrm{w}}}$	Tip chord of the wing, m	0.9
	${{\mathsf{\Lambda }}_{\mathrm{L}\mathrm{E}}}_{\mathrm{w}}$	Leading edge sweep angle of the wing, deg	3.2
	${{\mathsf{\Lambda }}_{\mathrm{T}\mathrm{E}}}_{\mathrm{w}}$	Trailing edge sweep angle of the wing, deg	−9.5
	lw	Distance, parallel to X-axis, from the leading edge of wing mean aerodynamic chord to the nose of fuselage, m	2.76
	zw	Distance, parallel to Z-axis, from the quarter chord of the wing mean aerodynamic chord to the X-axis, positive down, m	0.205
	$\mathsf{\Gamma }$	Dihedral angle, deg	7.5
Battery Placement	${\mathrm{l}}_{{\mathrm{v}}_{{\mathrm{b}}_1}}$	Battery set 1, distance from nose, m	0.9
	${\mathrm{l}}_{{\mathrm{v}}_{{\mathrm{b}}_2}}$	Battery set 2, distance from nose, m	6.5
	${\mathrm{N}}_{{\mathrm{b}}_1}$	Battery set 1, number of batteries	16
	${\mathrm{N}}_{{\mathrm{b}}_2}$	Battery set 2, number of batteries	8

.

The results for the aerodynamic characteristics of the optimized eVTOL aircraft, by using the proposed multidisciplinary possibilistic approach, are presented for cruise flight conditions at an airspeed of 70 m/s. Figure 12 depicts the results for the longitudinal static characteristics for the optimized eVTOL aircraft configuration.

Figure 13 presents the results for the lateral-directional static characteristics of the optimized eVTOL tilt-wing aircraft in cruise flight conditions.

Figure 14 shows the results for the longitudinal dynamic characteristics of the optimized eVTOL tilt-wing aircraft in cruise flight conditions.

In the following, the dynamic stability results for the optimized eVTOL tilt-wing aircraft in horizontal flight are presented. As can be seen, all periodical modes were stable for the eVTOL tilt-wing aircraft at the tested speeds and flight conditions. The stability characteristics were defined as constraints for the optimization model. The constraints were implemented considering a penalty for any scenarios out of the boundaries. It should be noted that all results have been considered for the maximum weight and the aftmost CG location as the worst-case scenario.

With respect to the CS-23 standards [67], the phugoid mode was noted as not strong. In Figure 15, the phugoid mode damping ratio considering the MIL airworthiness requirements for this aircraft type is shown [68,69]. As was visibly apparent, the damping characteristics completely satisfied the Level 1 handling quality for all considered speeds. With respect to Figure 16 presents the time to double for the phugoid amplitude. As can be seen, the oscillations were stable and negative all the time.

For Short Period oscillations, according to the CS-23 criteria, the characteristics were very significant. With respect to Figure 17 which shows the short period undamped natural frequency against different acceleration sensitivity values, the results met the Level 1 requirements according to MIL requirements [70]. Furthermore, Figure 18 presents the damping ratio values for the short period mode, which nicely satisfies Level 1 of handling quality according to the CS-23 criteria [34].

With respect to CS-23.181, the Dutch roll oscillations between the stalling speed and the maximum speed must be damped to 1/10 amplitude in seven cycles. Figure 19 demonstrates that this requirement has also been met.

With respect to Figure 20, the roll mode was effectively damped. Also, as shown in Figure 21 where the spiral mode characteristics are presented, the values for time to double show that the spiral mode is damped according to the requirements defined in CS-23 and MIL-F-8785-C [34,67]. Complying with both procedures, the spiral mode shows Level 1 of handling quality.

5. Conclusions

The projected multidisciplinary possibilistic method effectively presented the capability to optimize the eVTOL tilt-wing airplane by changing the geometry of the empennage package as well as wing and battery placement in the early design phases where several uncertainties exist because of low fidelity calculations. In this example, the disciplines of aerodynamics, propulsion, weight and balance, and stability and control were considered. With respect to different design characteristics, including longitudinal/lateral/directional trim, stability and control characteristics, and other criteria, the resulting eVTOL tilt-wing aircraft provided the optimized design with the desired characteristics to ensure that airworthiness requirements were met, and predefined constraints were complied with. The proposed approach may be used as a basis for light aircraft companies to upgrade their current products to the eVTOL versions. The resulting eVTOL design will be not only an optimized version but one where the stability of the aircraft is guaranteed. With this, the manufacturers of small companies can enhance the design and development of new aircraft without some of the costs associated with the older development methods, primarily because many of the required considerations in the conceptual design phase have been handled by the proposed approach.