A Multidisciplinary Possibilistic Approach to Size the Empennage of Multi-Engine Propeller-Driven Light Aircraft

Abstract

In considering aircraft design, it is very important to effectively size the tail configuration for stability and control. Multidisciplinary design optimization (MDO) focuses on the use of numerical optimization in the design of systems with multiple subsystems or disciplines of consideration. However, MDO uses deterministic calculations, and does not consider the uncertainties that arise from the employed analyses, including errors due to linearization and simplification. For problems with inadequate input data, the possibility-based design optimization (PBDO) scheme can be implemented in its stead to achieve reliable designs using membership functions for epistemic uncertainties. A multidisciplinary, possibilistic approach is presented to define the sizing of the empennage configuration of a twin-engine propeller-driven aircraft by changing shape parameters while satisfying the design requirements given the tailless aircraft configuration, the flight conditions, and various uncertainties. The corresponding disciplines are aerodynamics, stability and control, propulsion and weight and balance. Herein, different design requirements are considered including longitudinal/lateral/directional trim and stability characteristics, manufacturing and controllability criteria, handling qualities, operational requirements, airworthiness and survivability. The resulting aerodynamic characteristics and flight dynamic stability outcomes show that the optimized tail configuration for the proposed aircraft fully complied with airworthiness requirements and predefined constraints while considering several uncertainties due to the use of early-stage statistical estimations. The proposed approach can be used to enhance the preliminary design of multi-engine propeller-driven light aircraft where only low-fidelity, statistical estimations are available. The resulting output is not only an optimized aircraft configuration, but one where the stability of the design has been ensured. In this work, the aerodynamic characteristics have been determined using a validated semi-empirical program called MAPLA, developed for light aircraft designs and development in the preliminary design phase. Furthermore, the optimization framework consists of a deterministic optimizer that runs sequentially with a possibility assessment algorithm.

1. Introduction

Traditional civil aircraft design and development is mainly focused on providing vehicles that meet performance requirements while minimizing operational costs. This meant that a sequential strategy was required to address changing assumptions as the design matured, with no promise for an ideal design as the designer was faced with a complex multidisciplinary procedure [1]. Earlier designers were forced to rely on their creative and qualitative knowledge due to the lack of computational analysis and computer aided tools in the beginning of the conceptual design stage. As such, it was very difficult to predict the exact penalties of altering certain variables, with core decisions being made through the use of existing information and repetitive, time-consuming studies in a loop until an acceptable configuration was found. As time passed, the market demanded products that could be as low-price as possible to be affordable by individuals while providing a safe, reliable, and efficient flight. Consequently, current aircraft design is extremely multifaceted, accounting for the interaction of multiple disciplines such as aerodynamics, performance, stability and control, propulsion, structures, avionics, and even environmental impacts [2]. To realize this time-consuming process, a massive investment in different analyses is required to modernize the design process [3]. Even though many different tools, such as the Finite Element Method (FEM) and Computational Fluid Dynamics (CFD), have been introduced to accelerate the aircraft design process and enhance the accuracy of estimations, principally in the preliminary and detailed design stages, the basic idea remains unchanged: to provide a cost-effective, high-quality vehicle that meets the required criteria of the desired mission. Hence, the aircraft design process still relies heavily on engineering information, historical statistics and low fidelity studies [4,5,6,7,8,9].

Modern aircraft design and development is a complex and involved process, requiring a decision-making process that combines multidisciplinary sequential calculations with the goal of balancing cost requirements and mission capabilities, all while various system criteria are considered. Thanks to computing developments, specifically in the areas of simulation tools and design procedures, more knowledge is now available in the earlier life-cycle stages, including the conceptual and preliminary design. As such, fewer penalties are associated with varying design variables in the modern process, providing a more efficient baseline product prior to the manufacturing phase. This was largely enabled by the now-abundant options for high-fidelity simulation tools [10,11,12,13].

The aircraft design and development process has extensively used different optimization techniques at different stages, but regardless of the iterative pathway, the result will only be as accurate as the original analysis methods they were built upon. Overall, both analytical and numerical optimization techniques, in spite of being incredibly powerful problem-solving tools, face several limitations when used for real engineering problems. Simply put, due to the linearization of calculations, simplification of analysis, and implementation of statistical and semi-empirical methods, some errors appear that spread over the optimization procedure. As such, optimizations in the early design stages face several restrictions due to the lack of accurate modelling and high-fidelity analysis, and the possibility emerges that the optimized design using low-fidelity estimations may fail to meet the required mission performance characteristics [14,15]. Consequently, it is imperative to understand the kind of error an aircraft development may face in the early stages of the design and development process, especially for guaranteeing the mission requirements while balancing the compromises made for effectiveness and affordability. Furthermore, knowledge about a specific type of aircraft typically enhances this optimization process and aids in developing viable designs, using experience gained through previous projects and available information about the specifications of existing designs in the same category [16,17,18,19].

Multidisciplinary design optimization (MDO) focuses on the use of numerical optimization in the design of systems with multiple subsystems or disciplines of consideration. The main reason for using MDO is that the performance of a system is related not only to each individual discipline, but also to their interactions. With the help of the MDO approach in the early design stages, one can simultaneously enhance the design while decreasing both the time and cost of the entire design cycle [20,21]. MDO has been a promising procedure for more than three decades, but it is still not as widely implemented in all the design phases as was initially expected [22,23,24,25,26,27,28,29,30,31]. Thus far, several methods have been proposed for MDO including: Multiple-discipline-feasible (MDF), Individual discipline feasible (IDF), All-at-once (AAO), and Simultaneous Analysis and Design (SAND) [20,21,22,23,24,25,26,27,28,29,30,31,32]. Some MDO structures have also been proposed based on the IDF and MDF methods. The IDF-based methods include Collaborative Optimization (CO), Enhanced Collaborative Optimization (ECO), Quasiseparable Decomposition (QSD), Exact and Inexact Penalty Decomposition (IPD/EPD) and Analytical Target Cascading (ATC). The MDF-based methods include Concurrent subspace optimization (CSSO), Bi-level integrated system synthesis (BLISS), Multidisciplinary design optimization based on independent subspaces (MDOIS) and Asymmetric Subspace Optimization (ASO) [20] MDO can improve the pace of conceptual design for an aircraft, quickly recognizing the optimal case through the use of simplified analyses that are available at this stage [33,34,35,36]. At its core, however, MDO relies upon deterministic calculations, and does not consider the uncertainties which arise from the defining assumptions of the employed analyses, including errors due to linearization and simplification.

The defining assumptions used for physical modelling open the process up to a risk of failure when a project transitions from paper to reality; a risk that can only be minimized with higher fidelity estimation, but never completely removed. Several works have been done to study the behavior of aleatory and epistemic uncertainty. Among those methods, some have been used for design optimization problems including Interval Analysis, Fuzzy Numbers, and Probability Theory. Generally, interval analysis or fuzzy numbers are preferred for cases where inadequate information is available to precisely determine the probability density function (PDF) [37,38,39,40].

Interval modeling is commonly used as a non-probabilistic representation of uncertainty [41,42]. The general concept is that some value $\mathrm{Y}$ cannot be determined with perfect accuracy but is known to lie between the boundaries $\mathrm{A}$ and $\mathrm{B}$. Accordingly, any mathematical operation which is applied to $\mathrm{Y}$ can alternately be used on the interval of $\left[\mathrm{A},\mathrm{B}\right]$ instead. This produces an output interval, within which lies the solution associated with $\mathrm{Y}$. Interval analysis does not provide any notion of where the solution specifically is, only that it is within the boundaries. The fuzzy numbers concept extends the idea of interval analysis by adding a function which defines the degree of membership within the interval. Figure 1 shows a linear membership function. Interval analysis can additionally be expressed as a specific subset of fuzzy models with a binary membership function where 1 is membership, and 0 is non-membership [43,44,45].

Reliability Based Design Optimization (RBDO) is a methodology to consider probabilistic variables and parameters, which can be used as a method to consider uncertainty in design optimization problems [46]. Recently, RBDO has found frequent use in stochastic design optimization exercises where the input data are adequate to generate a precise numerical distribution. Comparatively, for problems with inadequate input data, the possibility-based design optimization (PBDO) method can be implemented in its stead to achieve reliable designs using membership functions for epistemic uncertainties [47].

For an aircraft design, it is very important to efficiently size the empennage configuration for stability and control [48]. Conventional tail configurations include both a horizontal tail and vertical tail, which during expected operation will experience only a portion of the total lift each can generate. This, specifically, results from the expectation that tail stall situations should never occur when the aircraft is in use [49]. To enforce this, most light aircraft feature a larger tail section than would be theoretically required for stable flight. Particularly, vertical tail surfaces are typically oversized to provide enough control in the single-engine failure scenario while being able to satisfy the requirements of dynamic motion flying qualities [50]. In addition, the safe-flight boundaries of the asymmetrically loaded airplane are of interest for propeller-driven airplanes within the terminal flight phases [51,52].

The main goal of this study is to implement a robust and efficient RBDO method for the multi-disciplinary design optimization of a twin-engine aircraft in the conceptual design stage. In doing so, the uncertainties which arose from using a multi-fidelity analysis, including semi-empirical estimations, would be adequately captured and ensured that the design would be viable in the detailed stages of the life cycle. In this specific case, the focus is on the empennage configuration, where the optimization proceeded by modifying the defining shape parameters whilst satisfying the design requirements imposed by the tailless aircraft configuration and the expected flight conditions. The involved disciplines are aerodynamics, stability and control, propulsion, as well as weight and balance.

In Section 2, the methodology of the proposed approach is presented with a focus placed on the tail and its importance for aircraft trim, stability and control (TSC). After that, the multidisciplinary tool used in this study and the optimization algorithm are introduced. Section 3 presents the MDO framework proposed for this study including the XDSM flowchart, design variables, cost objectives, constraints, and the sources for uncertainty. It also includes the corresponding results of the optimized aircraft including its resulting geometry, the aerodynamics characteristics and the flight dynamic stability. Finally, a short conclusion is provided in Section 5.

2. Methodology

The empennage in a conventional aircraft is responsible for TSC in the longitudinal and lateral-directional axes [53]. The horizontal stabilizer is responsible for the longitudinal TSC, while the vertical stabilizer is responsible for the lateral-directional TSC [48]. The design criteria that have to be guaranteed by the design process of an empennage are as follows: (1) longitudinal trim, (2) directional trim, (3) lateral trim, (4) longitudinal stability, (5) directional stability, (6) lateral stability, (7) manufacturability and controllability, (8) handling qualities, (9) operational requirements, (10) airworthiness, and (11) survivability. Usually, designing for stability is done at the expense of controllability, wherein a perfectly stable design would, by definition, be uncontrollable. Depending on the design needs, a difference balance may be struck between these two opposing priorities. For example, a commercial passenger aircraft will put a greater focus on stability, while a fighter jet will strongly favour control. The degree to which the stability and controllability of an airplane are defined and prioritized is available within civil aviation standards such as 14 CFR part 25, Military Standards (MIL-STDs), or other flight standards [53].

2.1. Aircraft TSC Requirements

When not done through an automated process, an empennage is first designed by considering the trim requirements, after which a stability and control analysis is used to modify it. When trimmed, an airplane does not rotate about the center of gravity (CG) and either moves only in a defined direction or circular path. Mathematically, the airplane can be said to be trimmed when the summation of forces and moments about the lateral ($\mathrm{x}$), longitudinal ($\mathrm{y}$), and vertical ($\mathrm{z}$) axes equal zero. After considering trim, the stability must be analyzed. This is the tendency of the airplane, as a response to some small disturbance, to return to the initial trimmed condition. It is broken down into two categories; static and dynamic, and considers the behaviour of the aircraft without any human or computer interference. Static stability considers only the initial response of the aircraft to some small and instantaneous perturbation, while dynamic stability analyzes the response over time. Additionally, control is defined as the ability of the airplane to change between different trim conditions. In a conventional design, the elevators are used to enable longitudinal control, ailerons are used for lateral control, and the rudder is used for directional control [53].

2.1.1. Longitudinal TSC

As shown in Equations (1) and (2), the airplane is considered to be in longitudinal trim when the summation of all forces in the $\mathrm{x}$-direction, including drag and thrust, is zero, and the summation of all moments about the $\mathrm{y}$-axis, including pitching moment, is zero [53]:

$$\sum {\mathrm{F}}_x=0$$

(1)

$$\sum {\mathrm{M}}_{\mathrm{CG}}=0$$

(2)

In this study, four disciplines were considered to size the aircraft empennage: aerodynamics, stability and control, propulsion, and weight and balance. For longitudinal TSC, using the MAPLA software package [54,55,56,57], vertical forces were modelled for different aerodynamic and stability characteristics including propwash. As shown in Figure 2, the tailless aircraft configuration creates a lifting force, ${\mathrm{Lift}}_{\mathrm{TL}}$, which is located at the aerodynamic center of ${\mathrm{AC}}_{\mathrm{TL}}$. As the aircraft CG is usually located aft of the aerodynamic center of the main wing, this lifting vector typically generates a positive moment about the CG which causes the aircraft nose to pitch up. For stable flight, the net moment must become zero by employing an upward force behind the CG. The horizontal stabilizer is used to provide this lifting force [48,49,50,51,52,53].

The use of the trim equation leads to Equation (3) as follows [53]:

$$\sum {\mathrm{M}}_{\mathrm{CG}}={\mathrm{M}}_{\mathrm{TL}}+{\mathrm{M}}_{\mathrm{h}}=0$$

(3)

where ${\mathrm{M}}_{\mathrm{T}\mathrm{L}}$ is the pitching moment generated by the tailless aircraft, which was calculated using the MAPLA software to model the wing, fuselage, nacelle, and propeller effects. ${\mathrm{M}}_{\mathrm{h}}$ is the horizontal tail contribution to pitching moment, which accounts for the propeller effects and downwash contribution on the horizontal tail. Equation (3) was nondimensionalized as follows:

$${\mathrm{C}}_{{\mathrm{m}}_{\mathrm{CG}}}={\mathrm{C}}_{{\mathrm{m}}_{\mathrm{TL}}}-{\mathsf{\eta }}_{\mathrm{h}}{\mathrm{V}}_{\mathrm{H}}{\mathrm{C}}_{{\mathrm{L}}_{{\mathrm{h}}_{\mathrm{p}}}}=0$$

(4)

where ${\mathrm{C}}_{\mathrm{m}}{}_{\mathrm{T}\mathrm{L}}$ is the tailless aircraft pitching moment coefficient and ${\mathrm{C}}_{{\mathrm{L}}_{{\mathrm{h}}_{\mathrm{p}}}}$ is lift coefficient of the horizontal tail. ${\mathrm{\eta }}_{\mathrm{h}}$ is the tail efficiency factor and is defined as a square of the effective airspeed at the horizontal tail over the aircraft airspeed. ${\mathrm{V}}_{\mathrm{H}}$ is the volume coefficient of the horizontal tail and can be expressed using Equation (5) as follows:

$${\mathrm{V}}_{\mathrm{H}}=\frac{{\mathrm{S}}_{\mathrm{h}}{\mathrm{l}}_{\mathrm{h}}}{\mathrm{S}\overline{\mathrm{c}}}$$

(5)

where ${\mathrm{S}}_{\mathrm{h}}$ is the horizontal tail area, ${\mathrm{l}}_{\mathrm{h}}$ is the distance from the $\mathrm{A}{\mathrm{C}}_{\mathrm{h}}$ to the CG, and $\mathrm{S}$ is the surface area of wing [53].

As longitudinal stability deals with the angular motion of the aircraft, the pertinent dynamic characteristic is the change of the pitching moment, $\mathrm{M}$, with respect to the angle of attack, $\mathsf{\alpha }$. Hence, as shown in Equation (6), the primary stability derivative used to model the static longitudinal stability is the rate of change of the pitching moment coefficient with respect to the change in $\mathsf{\alpha }$, denoted here as ${\mathrm{C}}_{\mathrm{m}}{}_{\mathsf{\alpha }}$. Moreover, with respect to Equation (7), the primary stability derivative that affects the dynamic longitudinal stability is the rate of change of ${\mathrm{C}}_{\mathrm{m}}$ with respect to the change in pitch rate, $\mathrm{q}$, presented here as ${\mathrm{C}}_{\mathrm{m}}{}_{\mathrm{q}}$ [53]:

$${\mathrm{C}}_{\mathrm{m}}{}_{\mathsf{\alpha }}=\frac{\partial {\mathrm{C}}_{\mathrm{m}}}{\partial \mathsf{\alpha }}$$

(6)

$${\mathrm{C}}_{\mathrm{m}}{}_{\mathrm{q}}=\frac{\partial {\mathrm{C}}_{\mathrm{m}}}{\partial \mathrm{q}}$$

(7)

2.1.2. Directional TSC

With respect to Equations (8) and (9), the airplane is defined to be directionally trimmed when the summation of all forces in the $\mathrm{y}$-direction and all moments about the $\mathrm{z}$-axis are equal to zero, including the effects of side forces and the yawing moment [53]:

$$\sum {\mathrm{F}}_{\mathrm{y}}=0$$

(8)

$$\sum {\mathrm{N}}_{\mathrm{CG}}=0$$

(9)

To model and ensure lateral-directional TSC, the rightward and leftward forces, along with propeller effects, were estimated using MAPLA. In considering directional behaviour, equilibrium refers to a rotational equilibrium for yawing motion, where yaw is defined to occur when the nose of the airplane is pointed away from the flight path within the $\mathrm{xy}$-plane. As shown in Figure 3, the sideslip angle $\mathsf{\beta }$ is used to represent this deviation from the flight path. This can be induced by a crosswind and to offset this angle, a moment $\mathrm{N}$ must be generated by the vertical tail. Unlike longitudinal trim where the airplane will be in equilibrium at a range of $\mathsf{\alpha }$, directional equilibrium is often at no sideslip, where $\mathsf{\beta }$ equals zero [48,49,50,51,52,53].

To maintain directional trim, a force in the $\mathrm{y}$-axis must be developed by the vertical tail, which generates a yawing moment. The resulting sum of moments can be expressed as:

$$\sum {\mathrm{N}}_{\mathrm{CG}}={\mathrm{N}}_{\mathrm{TL}}+{\mathrm{N}}_{\mathrm{v}}=0$$

(10)

where ${\mathrm{N}}_{\mathrm{TL}}$ is the tailless aircraft yawing moment including wing, fuselage and nacelle contribution calculated using the MAPLA software package. ${\mathrm{N}}_{\mathrm{V}}$ is the vertical tail contribution to the total yawing moment, with propeller effects and sidewash contribution accounted for. Equation (10) can be nondimensionalized and expressed as follows:

$${\mathrm{C}}_{\mathrm{n}}={\mathrm{C}}_{\mathrm{n}}{}_{\mathrm{TL}}+{\mathrm{V}}_{\mathrm{V}}{\mathrm{C}}_{{\mathrm{L}}_{{\mathrm{v}}_{\mathrm{p}}}}=0$$

(11)

where ${\mathrm{C}}_{\mathrm{n}}{}_{\mathrm{TL}}$ is the tailless aircraft yawing moment coefficient, and ${\mathrm{C}}_{{\mathrm{L}}_{{\mathrm{v}}_{\mathrm{p}}}}$ is a lift coefficient of the vertical tail where both consider propeller effects. ${\mathrm{V}}_{\mathrm{V}}$ is the volume coefficient of the vertical tail and can be calculated using Equation (10) as follows:

$${\mathrm{V}}_{\mathrm{V}}=\frac{{\mathrm{S}}_{\mathrm{v}}{\mathrm{l}}_{\mathrm{v}}}{\mathrm{S}\overline{\mathrm{c}}}$$

(12)

where ${\mathrm{S}}_{\mathrm{v}}$ is the vertical tail surface area, and ${\mathrm{l}}_{\mathrm{v}}$ is the distance from aerodynamic center of the vertical tail to the aircraft CG [53].

For multi-engine aircraft, the vertical tail must be designed to stabilize the yawing moment produced during one engine inoperative (OEI) operation. The vertical tail must also be sized for the spin recovery. In the case of multi-engine airplanes, however, spin recovery is not specifically considered as the aircraft are not likely to enter this situation [48]. Hence, Equation (10) can be expressed as follows for the OEI condition:

$${\mathrm{C}}_{\mathrm{n}}={\mathrm{C}}_{\mathrm{n}}{}_{\mathrm{TL}}+{\mathrm{C}}_{{\mathrm{T}}_{\mathrm{R}}}{\mathrm{Y}}_{\mathrm{T}}+{\mathrm{V}}_{\mathrm{v}}{\mathrm{C}}_{{\mathrm{L}}_{{\mathrm{v}}_{\mathrm{p}}}}=0$$

(13)

where ${\mathrm{C}}_{{\mathrm{T}}_{\mathrm{R}}}$ is the right engine thrust coefficient and ${\mathrm{Y}}_{\mathrm{T}}$ is the distance between the thrust line and the aircraft CG within the $\mathrm{xy}$-plane [53].

As the directional stability investigates the aircraft’s yawing motion, the corresponding pertinent dynamic characteristic is the change of the yawing moment, $\mathrm{N}$, with respect to the $\mathsf{\beta }$. Hence, with respect to Equation (14), the primary stability derivative to model static directional stability is the rate of change of the yawing moment coefficient, ${\mathrm{C}}_{\mathrm{n}}$, with respect to the change in the sideslip angle, denoted here as ${\mathrm{C}}_{{\mathrm{n}}_{\mathsf{\beta }}}$. Furthermore, with respect to Equation (15), the primary stability derivative that informs the dynamic directional stability is the rate of change of ${\mathrm{C}}_{\mathrm{n}}$ with respect to the change in the yaw rate, $\mathrm{r}$, and denoted here as ${\mathrm{C}}_{{\mathrm{n}}_{\mathrm{r}}}$ [48,49,50,51,52,53]:

$${\mathrm{C}}_{\mathrm{n}}{}_{\mathsf{\beta }}=\frac{\partial {\mathrm{C}}_{\mathrm{n}}}{\partial \mathsf{\beta }}$$

(14)

$${\mathrm{C}}_{\mathrm{n}}{}_{\mathrm{r}}=\frac{\partial {\mathrm{C}}_{\mathrm{n}}}{\partial \mathrm{r}}$$

(15)

2.1.3. Lateral TSC

With respect to Equations (16) and (17), the airplane is said to be in a state of lateral trim when the summation of all forces in the $\mathrm{z}$-direction and all moments about the $\mathrm{x}$-axis is equal to zero. These considerations include the lift, weight, and rolling moment experienced by the aircraft [53]:

$$\sum {\mathrm{F}}_{\mathrm{z}}=0$$

(16)

$$\sum {\mathrm{L}}_{\mathrm{CG}}=0$$

(17)

In Figure 4, the front view of an aircraft is shown where the vertical tail is generating a counteracting moment as a response to the tailless aircraft developing a rolling moment from a nonzero sideslip angle. An additional moment is generated by the rotating propellers, which is also counteracted by the vertical tail. As with directional trim, lateral equilibrium is usually present at a condition with no sideslip, where $\mathsf{\beta }=0$ [48,49,50,51,52,53].

To maintain lateral trim, the vertical tail must generate a moment about the $\mathrm{x}$-axis (${\mathrm{L}}_{\mathrm{v}}$), which is equal in magnitude to the rolling moment generated by the tailless body. Thus, the lateral trim equations can be expressed as follows:

$$\sum {\mathrm{L}}_{\mathrm{CG}}={\mathrm{L}}_{\mathrm{TL}}+{\mathrm{L}}_{\mathrm{V}}$$

(18)

where ${\mathrm{L}}_{\mathrm{TL}}$ is the tailless aircraft rolling moment calculated by MAPLA, which includes the wing, fuselage, nacelles, and propeller contributions. ${\mathrm{L}}_{\mathrm{V}}$ is the vertical tail contribution to the total rolling moment, which also contains propeller and sidewash effects. Equation (18) can be nondimensionalized and expressed as follows:

$${\mathrm{C}}_{\mathrm{l}}={\mathrm{C}}_{\mathrm{l}}{}_{\mathrm{TL}}+\frac{{\mathrm{S}}_{\mathrm{V}}{\mathrm{z}}_{\mathrm{V}}}{\mathrm{S}\overline{\mathrm{c}}}{\mathrm{C}}_{{\mathrm{L}}_{{\mathrm{V}}_{\mathrm{p}}}}=0$$

(19)

where ${\mathrm{z}}_{\mathrm{V}}$ is the distance from the AC of the vertical tail to the aircraft CG [53].

As lateral stability is concerned with the aircraft’s rolling motion, the corresponding pertinent dynamic characteristic is the change of the rolling moment, $\mathrm{l}$, with respect to the sideslip angle β. With respect to Equation (20), the primary stability derivative to calculate the static lateral stability is the rate of change of yawing moment coefficient, ${\mathrm{C}}_{\mathrm{l}}$, with respect to the change in the sideslip angle, denoted here as ${\mathrm{C}}_{{\mathrm{l}}_{\mathsf{\beta }}}$ [48,49,50,51,52,53]:

$${\mathrm{C}}_{\mathrm{l}}{}_{\mathsf{\beta }}=\frac{\partial {\mathrm{C}}_{\mathrm{l}}}{\partial \mathsf{\beta }}$$

(20)

More information is available in the flight dynamics books by Roskam [58], Nelson [59] or Etkin [60].

2.2. MAPLA

The enhanced semi-empirical multidisciplinary analysis program MAPLA was used to model the relevant disciplines. This tool is specifically built for design optimization of light, general aviation, propeller-driven aircraft, and includes analyses for four disciplines: aerodynamics, propulsion, performance, as well as stability and control. MAPLA includes state-of-the-art analytical procedures and design data collections, which have been combined and modified in a unique method, which is fully automated. Previous investigations demonstrated that MAPLA was able to determine the aerodynamic characteristics of light aircraft with an acceptable accuracy in various configurations [54,55,56,57]. While this program performs well in modelling light general aviation aircraft with conventional configuration, which was the focus of its development, the module’s use of semi-empirical methods has the unfortunate consequence of it not being generally applicable. For the purposes of this investigation, where a conventional light aircraft was targeted, MAPLA was more than sufficient.

2.3. PBDO Method Outline

The optimization algorithm contains a deterministic optimizer which works sequentially with a possibility–evaluation algorithm. This framework is comparable to the Sequential Optimization and Reliability Assessment (SORA) procedure developed for RBDO [61,62]. The Performance Measure Approach (PMA) was used to implement the possibility evaluation on each constraint, which was followed by a deterministic optimization. The PMA algorithm consists of a separate optimization loop which finds a worse-case scenario for each constraint separately. The position of each uncertain variable and parameter within their intervals are determined such that that the constraint is maximized. This guarantees that constraints will be satisfied even in the extreme worst-case where each individual source of uncertainty drives the designs towards constraining limits [63,64]. A Sequential Quadratic Programming (SQP) procedure is used in both the deterministic optimization and the PMA-based possibility evaluation. The method accomplishes deterministic optimizations in sequence with possibility evaluations, altering the constraints until convergence can be guaranteed [65].

To validate the MDO framework, the following nonlinear single-discipline RBDO test problem was solved with different methods, after which the results were compared [65,66]:

$$\begin{gathered} \min \mathrm{f}={\overline{\mathrm{x}}}_1+{\overline{\mathrm{x}}}_2\mathrm{st}. \mathrm{P}\left[{\mathrm{G}}_{\mathrm{j}}\left(\mathrm{x}\right)\ge 0\right]\ge {\mathrm{P}}_{\mathrm{goal}, \mathrm{j}}, \mathrm{j}=1,2,3 \\ {\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}\left({\mathrm{x}}_2, 0.3\right) \end{gathered}$$

(21)

The problem was solved for the constraint boundaries of ${\mathrm{G}}_1$, ${\mathrm{G}}_2$, and${ \mathrm{G}}_3$ with different reliability indices ${\mathsf{\beta }}_{\mathrm{r}}$ of 1 to 5 with the double loop method using Reliability Index Approach (RIA) and PMA, the sequential method with PMA, and the single-loop method [65]. The algorithm performance was compared at ${\mathsf{\beta }}_{\mathrm{r}}=3$ for all methods and is shown in

Table 1. Performance comparison of different optimization methods.

Method	${\mathit{f}}_{\mathit{m}\mathit{i}\mathit{n}}$	${\mathit{x}}_1$	${\mathit{x}}_2$	${\mathit{n}}_{\mathit{e}\mathit{v}\mathit{a}\mathit{l}\mathit{s}}$
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

as follows:

It was found that the double loop and sequential methods all yielded approximately equal results but required varying numbers of evaluation to arrive at these results. As shown in Table 1, the PMA-based methods solve required fewer evaluations and by extension fewer computation resources. Of the two PMA methods, the sequential implementation was noted to require less function evaluations than the double loop approach. For this problem, the single loop approach was found to need the smaller number of evaluations but demonstrated a divergence from the PMA and RIA methods when ${\mathrm{\beta }}_{\mathrm{t}}$ was increased.

The stability, speed, and accuracy of each method, shown in

Table 2. Stability comparison of different methods.

Method	No. of Failed Runs	Average Error %	Median Error %	Average Time (s)	Average Evals
Single Loop	0	1.93	1.25	0.3343	145
PMA/Sequential	5	2.54 × 10–5	2.54 × 10–5	0.5805	1648
PMA/Double Loop	0	3.77	2.97 × 10–5	2.7277	3862
RIA/Double Loop	71	6.11	0.319	4.1316	13276

, were also evaluated by running the optimization over a grid of 100 points between [0;0] and [10;10].

3. MDO Framework

3.1. Methods

MDO methods are classified into tw7o different categories: single-level and multilevel methods. Of the seven methods mentioned prior, MDF, IDF, AAO, SAND and MDOIS are categorized as single-level methods which implement a single optimizer and directly use a non-hierarchical structure. The remaining methods are the multilevel methods: CSSO, BLISS and CO. For the single-level methods, with the exception of the MDOIS approach, coupling variables are updated by the individual disciplines within a global iteration. The MDOIS approach only updates these coupling values after all disciplines have finished computation and a global multidisciplinary analysis is run. In considering the multilevel methods, the core difference is that the architectures include a hierarchical structure with a separate optimizer at each level [32].

In this paper, the RyeMDO IDF/Sequential/PMA framework with an SQP optimization algorithm has been implemented. A framework was developed which implements a PBDO algorithm to find the optimum empennage design configuration for the twin-engine propeller driven aircraft. The algorithm consists of a PBDO solver coupled with MAPLA to provide the four analyses. In the previous work [65], this framework was found to be reliable and less dependent on constraint functions and coarse objectives than the double loop procedures. Moreover, this framework was shown to have a smaller dependency on the starting vector selection compared with the single loop procedure, whilst being more effective than all methods save for the single loop approach in the considered cases. The sequential approach was chosen as it provides more benefits when solving for multi-level reliability [65].

The sequential framework resolves a sequence of deterministic optimizations first, then carries out reliability evaluations until convergence is achieved. This results in the intermediate solutions gradually passing through lower indices until the target level is reached [65]. In this work, the low reliability index was applied at the first step and then gradually higher indices were used. This allows the approach to extract intermediate answers and plot the connection between the optimal design and nominated target reliability index without significantly increasing the number of function evaluations. The objective function was defined to minimize the empennage package weight for the specified flight conditions while lowering the reliability index of the designs.

In this paper, for the reliability assessment sub-problems, the MDF approach was selected instead of the IDF method. This was justified on the grounds that, for problems with a large number of variables and fewer uncertainties, the dimensionality of the reliability analysis is significantly larger in the IDF framework. Moreover, the auxiliary equality constraints had to be considered for coupling variable consistency, producing significant complexity in the outcome that must be solved for every constraint in each iteration of the sequential procedure. By using the MDF approach for reliability assessments, the dimensionality of the PMA optimization was reduced and no auxiliary constraints had to be considered. In previous work [65], it was shown that the PMA reliability evaluations converged sooner for the MDF approach and with higher reliability.

With respect to the four different disciplines considered in this study (aerodynamics, propulsion, mass and balance, as well as stability and control), the XDSM flowchart in Figure 5 is presented for the multi-disciplinary optimization procedure of the twin-engine propeller-driven aircraft. To run the system-level optimizer, each sublevel was analyzed sequentially with the output of one discipline providing input to another. Using this approach, the empennage surface area was optimized based on the geometric parameters and the corresponding flight condition achieved by each discipline. The method is outlined as follows and the corresponding block diagram is shown in Figure 5:

1.

Find the analysis error distributions

(I)

Model each entry based on the aircraft empennage model

(II)

Determine the entire aircraft aerodynamic characteristics

(III)

Implement the uncertainty analysis

(IV)

Implement a best fit PDF curve for each individual source of uncertainty

2.

Choose the starting vector and a list of reliability indices

3.

Run the IDF-based deterministic optimization

4.

Run the MDF-based PMA reliability evaluation at the current reliability index, and alter variables according to the sequential procedure

5.

Check for convergence with current reliability goal

(I)

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

(II)

if no, return to 3

6.

Advance to next target reliability level

(I)

retain solution as new starting vector

(II)

return to 3

3.2. Disciplines

Four disciplines of analysis were considered in this investigation. The first of which was the aerodynamic analysis, which consists of longitudinal and lateral-directional aerodynamic coefficients and stability derivatives. Power-off aerodynamic coefficients are initially estimated for each section of the tailless aircraft, including the wing, fuselage and nacelle. The empennage package is subsequently introduced to consider the horizontal tail, vertical tail and high-lift surface contributions. Following this, the static stability derivatives are calculated and the total power-off characteristics are obtained through combining the contributions of each aircraft part. Afterwards, the dynamic stability characteristics of the aircraft are estimated by using the stability derivatives. Finally, the power-on parameters and propeller effects were added to the estimations through the use of semiempirical methods.

After the aerodynamic analysis concludes, the weight of each component is calculated through a statistical method based on the MIL-STD-1374 via a “Summary Group Weight Statement” [18,68,69]. Both the mass and CG of individual components are calculated by statistical equations to provide the overall empty weight of the airplane [18,68,69].

Aircraft handling quality is evaluated using the flying quality subprogram of MAPLA. This module enables the estimation of the aircraft trim characteristics in all flight and aircraft configurations, thereby providing the corresponding handling quality level for different longitudinal and lateral-directional modes in accordance with the standards of this aircraft type.

3.3. Design Variables

As this study’s main concern was to size the empennage of a twin-engine propeller-driven light aircraft, the corresponding design variables were selected based on the effectiveness of each component with respect to the semi-empirical method estimations. These variables are shown in Figure 6. In addition to the tail components, the engine placement was also parametrized due to the impact on rudder design. The engine position is represented as the lateral position on the wing, measured as a $\mathrm{y}$-direction distance from the $\mathrm{x}$-axis. The full list of design variables is shown in

Table 3. Aircraft design variables and their descriptions for empennage sizing.

Component	Variable	Description	Limits
Horizontal Tail	${\mathrm{i}}_{\mathrm{h}}$	Horizontal tail incidence angle, deg	0 to 3
	bh	Horizontal tail span, m	3.5 to 5.1
	${\mathrm{c}}_{{\mathrm{r}}_{\mathrm{h}}}$	Horizontal tail root chord, m	1 to 1.45
	${\mathrm{c}}_{{\mathrm{t}}_{\mathrm{h}}}$	Horizontal tail tip chord, m	0.5 to 1
	${\mathrm{\Lambda }}_{\mathrm{L}\mathrm{E}}{}_{\mathrm{h}}$	Leading edge sweep angle of horizontal tail, deg	10 to 20
	lh	Distance, parallel to X-body axis, from the nose of fuselage to the horizontal tail mean aerodynamic chord, m	8 to 8.7
	zh	Distance, parallel to Z-body axis, from the X-body axis to the quarter chord of the horizontal tail mean aerodynamic chord, positive down, m	–0.4 to –0.1
	celevator	Ratio of elevator chord to horizontal tail chord	0.2 to 0.5
Vertical Tail	bv	Vertical tail span, m	1.8 to 2.2
	${\mathrm{c}}_{{\mathrm{r}}_{\mathrm{v}}}$	Vertical tail root chord, m	1.9 to 2.5
	${\mathrm{c}}_{{\mathrm{t}}_{\mathrm{v}}}$	Vertical tail tip chord, m	0.8 to 1.4
	${\mathrm{\varphi }}_{\mathrm{T}\mathrm{E}}$	Trailing edge sweep angle of vertical tail, deg	10 to 20
	zv	Perpendicular distance from X-body axis to root chord of vertical-tail, positive down, m	–0.35 to –0.15
	lv	Distance along X-body axis from the nose of fuselage to leading edge of tip chord of vertical tail, m	8.5 to 9.5
	crudder	Ratio of rudder chord to vertical tail chord, m	0.2 to 0.5
Engine	YT	Lateral Distance from X-axis to thrust line, m	1.6 to 1.9

.

3.4. Objective Functions

Two scenarios are considered for the objective function. For the first scenario, the weight optimization was investigated and, for the second scenario, the optimization targeted the minimum drag in a trimmed flight condition. In the following, both scenarios were described and the corresponding results are shown in Section 4.

For the first scenario, the objective function for the optimization process is to minimize the weight of the horizontal and vertical stabilizers by varying the design variables shown in Table 1. The resulting equation has a direct relation with the surface areas of the horizontal and vertical tail, which allows the optimization scheme to produce a design with minimal area.

The horizontal tail weight, ${\mathrm{W}}_{\mathrm{h}}$, was estimated with a statistical relationship provided for the general aviation aircraft type using Equation (22) as follows [17]:

$${\mathrm{W}}_{\mathrm{h}}=0.016{\left({\mathrm{N}}_{\mathrm{z}}{\mathrm{W}}_{{\mathrm{d}}_{\mathrm{g}}}\right)}^{0.414}{\mathrm{q}}^{0.0168}{\mathrm{S}}_{\mathrm{h}}^{0.896}{\left(\frac{100{\mathrm{t}}_{{\mathrm{c}}_{\mathrm{h}}}}{{\cos \mathsf{\Lambda }}_{{\frac{ \overline{\mathrm{c}}}{4}}_{\mathrm{h}}}}\right)}^{-0.12}{\left(\frac{{\mathrm{AR}}_{\mathrm{h}}}{{\cos }^2{\mathsf{\Lambda }}_{{\frac{ \overline{\mathrm{c}}}{4}}_{\mathrm{h}}}}\right)}^{0.043}{\mathsf{\lambda }}_{\mathrm{h}}^{-0.02}$$

(22)

where ${\mathrm{N}}_{\mathrm{z}}$ is the ultimate load factor and is equal to 1.5 times of aircraft load factor limit, $\mathrm{q}$ is the dynamic pressure ratio, ${\mathrm{W}}_{{\mathrm{d}}_{\mathrm{g}}}$ is the flight design gross weight, ${\mathrm{t}}_{{\mathrm{c}}_{\mathrm{h}}}$ is the horizontal tail airfoil section thickness ratio, ${\mathsf{\Lambda }}_{{\frac{ \overline{\mathrm{c}}}{4}}_{\mathrm{h}}}$ is the sweep angle of the horizontal tail at ¼ of the mean aerodynamic chord, ${\mathrm{AR}}_{\mathrm{h}}$ is the aspect ratio of the horizontal tail and ${\mathsf{\lambda }}_{\mathrm{h}}$ is the taper ratio of the horizontal tail.

The surface area of the horizontal tail was calculated using Equation (23) based on the geometrical parameters provided in Figure 6 as follows:

$${\mathrm{S}}_{\mathrm{h}}={\mathrm{c}}_{{\mathrm{r}}_{\mathrm{h}}}{\mathrm{b}}_{\mathrm{h}}-{\mathrm{c}}_{{\mathrm{e}}_{{\mathrm{h}}_{\mathrm{for}}}}\frac{{\mathrm{b}}_{\mathrm{h}}}{2}-{\mathrm{c}}_{{\mathrm{e}}_{{\mathrm{h}}_{\mathrm{back}}}}\frac{{\mathrm{b}}_{\mathrm{h}}}{2}$$

(23)

where ${\mathrm{c}}_{{\mathrm{e}}_{{\mathrm{h}}_{\mathrm{for}}}}$ is the extended panel out of the tip chord of the horizontal tail and can be found in Figure 6 using the trailing edge angle of the horizontal tail as follows:

$${\mathrm{c}}_{{\mathrm{e}}_{{\mathrm{h}}_{\mathrm{for}}}}=\frac{{\mathrm{b}}_{\mathrm{h}}}{2}\tan \left({\mathsf{\Lambda }}_{{\mathrm{LE}}_{\mathrm{h}}}\right)$$

(24)

In addition, ${\mathrm{c}}_{{\mathrm{e}}_{{\mathrm{h}}_{\mathrm{back}}}}$ can be calculated using the

$${\mathrm{c}}_{{\mathrm{e}}_{{\mathrm{h}}_{\mathrm{back}}}}=\frac{{\mathrm{b}}_{\mathrm{h}}}{2}\tan \left({\mathsf{\Lambda }}_{{\mathrm{TE}}_{\mathrm{h}}}\right)$$

(25)

where ${\mathsf{\Lambda }}_{{\mathrm{TE}}_{\mathrm{h}}}$ is the trailing edge sweep angle of the horizontal tail and can be found using Equation (26) as follows:

$${\mathsf{\Lambda }}_{{\mathrm{TE}}_{\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{for}}}}\right)}\right)$$

(26)

Following this, the vertical tail weight, ${\mathrm{W}}_{\mathrm{v}}$, was calculated using the following statistical equation [18]:

$${\mathrm{W}}_{\mathrm{v}}=0.073\left(1+0.2\frac{{\mathrm{z}}_{\mathrm{h}}}{{\mathrm{z}}_{\mathrm{v}}}\right){\left({\mathrm{N}}_{\mathrm{z}}{\mathrm{W}}_{{\mathrm{d}}_{\mathrm{g}}}\right)}^{0.376}{\mathrm{q}}^{0.122}{\mathrm{S}}_{\mathrm{v}}^{0.3873}{\left(\frac{100{\mathrm{t}}_{{\mathrm{c}}_{\mathrm{v}}}}{{\cos \mathsf{\Lambda }}_{{\frac{ \overline{\mathrm{c}}}{4}}_{\mathrm{v}}}}\right)}^{-0.49}{\left(\frac{{\mathrm{AR}}_{\mathrm{v}}}{{\cos }^2{\mathsf{\Lambda }}_{{\frac{ \overline{\mathrm{c}}}{4}}_{\mathrm{v}}}}\right)}^{0.375}{\mathsf{\lambda }}_{\mathrm{v}}^{0.039}$$

(27)

where ${\mathrm{t}}_{{\mathrm{c}}_{\mathrm{v}}}$ is the vertical tail airfoil section thickness ratio, ${\mathsf{\Lambda }}_{{\frac{ \overline{\mathrm{c}}}{4}}_{\mathrm{v}}}$ is the sweep angle of the vertical tail at one quarter of the mean aerodynamic chord, ${\mathrm{AR}}_{\mathrm{v}}$ is the aspect ratio of the vertical tail and ${\mathsf{\lambda }}_{\mathrm{v}}$ is the taper ratio of the vertical tail.

Similarly, the surface area of the vertical tail was computed using the following equation, which was again based on the geometrical parameters provided in Figure 6:

$${\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\tan \left({\mathsf{\varphi }}_{\mathrm{TE}}\right)}{2}\right)$$

(28)

where ${\mathrm{c}}_{{\mathrm{e}}_{\mathrm{v}}}$ is the extended panel out of the root chord of the vertical tail and was found with the trailing edge sweep angle as follows:

$${\mathrm{c}}_{{\mathrm{e}}_{\mathrm{v}}}={\mathrm{b}}_{\mathrm{v}}\tan \left({\mathsf{\varphi }}_{\mathrm{TE}}\right)$$

(29)

For the second scenario, MAPLA’s stability and control module was used during the optimization process to find the trim condition and its corresponding drag force, which was used as an objective function to be minimized.

3.5. Constraints

First and foremost, manufacturing constraints were employed in defining the lower and upper limits for each design variable. The total empty mass of the aircraft was considered to be below 2000 kg according to statistical data and an initial weight and balance analysis. Static and dynamic criteria were also set as the constraints in the optimization process. The effectiveness of the control surfaces during cruise and climb-out were evaluated to ensure acceptable pitch authority as achieved in slow flight and that yaw authority was maintained in single-engine flight.

To have longitudinal static stability, ${\mathrm{C}}_{\mathrm{m}}{}_{\mathsf{\alpha }}$ must be negative. The typical value for most airplanes is between –0.3 to –1.5 1/rad. For longitudinal dynamic stability, the real part of the roots of the longitudinal characteristic equation must be negative. A major contributor to this criterion is ${\mathrm{C}}_{\mathrm{m}}{}_{\mathrm{q}}$, such that a negative value has a robust stabilizing effect. The typical value of ${\mathrm{C}}_{\mathrm{m}}{}_{\mathrm{q}}$ for most airplanes is between –5 to –30 1/rad [53]. Additionally, longitudinal stability was constrained to a static margin of 5%.

To achieve roll stability, the value for ${\mathrm{C}}_{{\mathrm{l}}_{\mathsf{\beta }}}$ must be negative [48]. Furthermore, to maintain directional static stability, ${\mathrm{C}}_{{\mathrm{n}}_{\mathsf{\beta }}}$ must be positive. The typical value for most airplanes is around +0.1 to +0.4 1/rad. For directional dynamic stability, the real part of the roots of the lateral-directional characteristic equation must be negative. A major contributor to this criterion is ${\mathrm{C}}_{\mathrm{n}}{}_{\mathrm{r}}$, such that a negative value has a robust stabilizing effect. The typical value for most airplanes is about –0.1 to –1 1/rad [53]. Yaw static stability must be ensured for a full-thrust climb scenario with one failed engine. The constraint ensures that the rudder and vertical tail have acceptable authority. All these constraints are considered in the optimization algorithm [70,71,72,73].

In addition, airworthiness requirements for light aircraft, particularly from CS-23 and available information from military standards for dynamic modes, were set as the constraints in the optimization process. The corresponding constraints for each mode were discussed in the Results section. The list of all constraints is provided in

Table 4. The list of constraints considered for the size of the empennage model.

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

as follows.

3.6. Sources of Uncertainty

The uncertainty due to the use of semi-empirical approaches in the aerodynamics and stability disciplines are somewhat difficult to measure. However, recent computing power enhancements can provide an answer with ~80% accuracy using the initial structure given by the initial aircraft layout [74]. From this, the error term was assumed to be uniformly distributed between the highest and lowest observed error. The error term was defined as the ratio of predicted characteristics to the observed characteristics from the database, ${\mathsf{\epsilon }}_{\mathrm{aero}}$. The characteristics from the aerodynamics discipline was then scaled by this ratio to correct for this error as shown in Equation (30). The error term was taken to be a uniformly distributed random parameter:

$$\begin{gathered} {\mathsf{\epsilon }}_{\mathrm{aero}}=\frac{{\mathrm{aero}}_{{\mathrm{characteristics}}_{\mathrm{p}}}}{{\mathrm{aero}}_{\mathrm{characteristics}}} \\ {\mathrm{aero}}_{\mathrm{characteristics}}={\mathsf{\epsilon }}_{\mathrm{aero}}{\mathrm{aero}}_{{\mathrm{characteristics}}_{\mathrm{p}}} \end{gathered}$$

(30)

The designs from the compiled aircraft database were modeled and the predicted empty-mass of each design was compared with the observed value from the database. The histogram of the errors indicated that the empty-mass error term could be approximated by a normal distribution curve. The mass error term was defined as shown in Equation (31). Empty-mass calculations carried out by the optimizer were scaled by the error term to account for the uncertainty in the prediction of the aircraft empty mass:

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

(31)

Furthermore, a previous study [56] showed that the required rudder deflection in OEI scenarios could be signified due to the asymmetric blade effect of propellers. It was shown that a maximum of 4 degrees of rudder deflection was expected to be considered in the design of the rudder surface for light aircraft. Accordingly, rudder deflection error was approximated by considering the asymmetric blade thrust effect. The propeller effect error term was defined in Equation (32) as follows. Maximum rudder deflection calculations carried out during the optimization process were scaled by the error term to account for the uncertainty in the prediction of the required rudder deflection:

$${\mathsf{\delta }}_{{\mathrm{r}}_{\max }}={\mathsf{\epsilon }}_{\mathrm{asym}}{\mathsf{\delta }}_{{\mathrm{r}}_{{\max }_{\mathrm{p}}}}$$

(32)

4. Results

In this section, the resulting geometry of the optimized empennage configuration, as well as the aerodynamic characteristics and dynamic stability of the twin engine propeller driven aircraft, are shown for both aforementioned objective scenarios.

4.1. Resulting Geometry

After running the optimization algorithm, the resulting geometry for the empennage configuration, considering the aforementioned objective scenarios, constraints and uncertainties, is shown in Figure 7a and the complete aircraft can be seen in Figure 7b, produced with MAPLA’s geometry module. The corresponding variables results were also presented in

Table 5. Optimized design variables for empennage sizing of the twin-engine aircraft.

Variable	Optimized Value with Mass Optimization	Optimized Value with Drag Optimization
${\mathrm{i}}_{\mathrm{h}}$ (deg)	1.95	1.96
${\mathrm{b}}_{\mathrm{h}}$ (m)	4.46	4.51
${\mathrm{c}}_{{\mathrm{r}}_{\mathrm{h}}}$ (m)	1.27	1.31
${\mathrm{c}}_{{\mathrm{t}}_{\mathrm{h}}}$ (m)	0.863	0.837
${\mathrm{\Lambda }}_{\mathrm{L}\mathrm{E}}{}_{\mathrm{h}}$ (deg)	12.18	11.23
${\mathrm{l}}_{\mathrm{h}}$ (m)	8.338	8.327
${\mathrm{z}}_{\mathrm{h}}$ (m)	–0.275	–0.281
${\mathrm{c}}_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{r}}$ (m)	0.37	0.37
${\mathrm{b}}_{\mathrm{v}}$ (m)	1.855	1.83
${\mathrm{c}}_{{\mathrm{r}}_{\mathrm{v}}}$ (m)	1.969	1.966
${\mathrm{c}}_{{\mathrm{t}}_{\mathrm{v}}}\left(\mathrm{m}\right)$	0.895	0.837
${\mathrm{\varphi }}_{\mathrm{T}\mathrm{E}}$ (deg)	17.095	16.94
${\mathrm{z}}_{\mathrm{v}}$ (m)	–0.221	–0.225
${\mathrm{l}}_{\mathrm{v}}$ (m)	9.09	9.05
${\mathrm{c}}_{\mathrm{r}\mathrm{u}\mathrm{d}\mathrm{d}\mathrm{e}\mathrm{r}}$ (m)	0.4135	0.415
${\mathrm{Y}}_{\mathrm{T}}$ (m)	1.706	1.701
${\mathrm{S}}_{\mathrm{h}}\left({\mathrm{m}}^2\right)$	4.755	4.84
${\mathrm{S}}_{\mathrm{v}}\left({\mathrm{m}}^2\right)$	2.656	2.565

. As can be seen, the resulting aircraft empennage sizing using both scenarios of mass minimization and the drag optimization in trim flight were similar.

4.2. Aerodynamics Characteristics of the Optimized Aircraft

MAPLA’s aerodynamics characteristics were compared with experimental data for a twin-engine light general aviation aircraft provided by NASA with similar characteristics to the one studied here [75,76]. Moreover, MAPLA has been used for design optimization of similar single and twin-engine propeller-driven light aircraft (e.g., two-, four- and six-seat propeller driven aircraft) [54,55,56,57]. The results demonstrated that MAPLA was able to determine the aerodynamic characteristics of the light aircraft with an acceptable accuracy. Aerodynamic characteristics of the aforementioned twin-engine aircraft are shown in the cruise flight condition, using the optimized tail model developed by implementing the proposed multidisciplinary possibilistic approach. The dotted-line results presented in the following figures are representative of the case with the drag optimization scenario.

Figure 8 shows the longitudinal static characteristics for both the tailless and whole aircraft configurations. As demonstrated in Figure 8c, the total aircraft configuration provided a value of approximately –1 1/rad for ${\mathrm{C}}_{{\mathrm{m}}_{\mathsf{\alpha }}}$ using the optimized tail configuration, which matches best within the typical value range discussed previously of $–0.3\le {\mathrm{C}}_{{\mathrm{m}}_{\mathsf{\alpha }}}\le –1.51/\mathrm{rad}$. In addition, a small positive value for the ${\mathrm{C}}_{{\mathrm{m}}_0}$ was expected, which in this case was equal to 0.05. As can be seen from the results shown in Figure 8, both optimization results show the same behavior.

Figure 9 shows the lateral-directional static characteristics for both the tailless and total aircraft configurations. As can be seen in Figure 9c, the complete aircraft configuration provided a value of about 0.12 1/rad for ${\mathrm{C}}_{{\mathrm{n}}_{\mathrm{\beta }}}$ using the optimized tail configuration which aligns with the typical value range discussed earlier of $0.1\le {\mathrm{C}}_{{\mathrm{n}}_{\mathrm{\beta }}}\le 0.41/\mathrm{rad}$. In addition, for roll stability, the value for ${\mathrm{C}}_{{\mathrm{l}}_{\mathrm{\beta }}}$ must be negative, and, for this aircraft, it was found to be approximately $–0.11/\mathrm{rad},$ which was considered acceptable. Similar to the longitudinal characteristics, the results shown in Figure 9 reveal the same behavior for both optimization scenarios with the mass optimization and minimum drag optimization in the trim flight condition.

Figure 10 presents the longitudinal dynamic characteristics for both tailless and total aircraft configurations. As can be seen in Figure 10c, the total aircraft configuration provided a value of approximately –18 1/rad for ${\mathrm{C}}_{{\mathrm{m}}_{\mathrm{q}}}$ using the optimized tail configuration, which complies with the typical value range discussed earlier, $–5\le {\mathrm{C}}_{{\mathrm{m}}_{\mathrm{q}}}\le –301/\mathrm{rad}$. Additionally, the results clearly show the proper contribution of the tail configuration in the dynamic stability of the proposed twin-engine aircraft design.

Figure 11 shows the lateral-directional dynamic characteristics for both the tailless and total aircraft configurations. As can be seen in Figure 11d, the total aircraft configuration provided a value of approximately 0.1 1/rad for ${\mathrm{C}}_{{\mathrm{n}}_{\mathrm{r}}}$ using the optimized tail configuration, which also lies within the typical value range discussed earlier, $-0.1\le {\mathrm{C}}_{{\mathrm{n}}_{\mathrm{r}}}\le -11/\mathrm{rad}$. Similar to the static aerodynamic characteristics, the results shown for the dynamic characteristics in Figure 10 and Figure 11 reveal the same behavior for both optimization scenarios with the mass optimization and minimum drag optimization in the trim flight condition.

4.3. Flight Dynamic Stability of the Optimized Aircraft

In this section, the results of the dynamic stability analysis obtained for the optimized aircraft are discussed in greater detail to investigate the dynamic stability requirements of the proposed twin-engine propeller driven aircraft. Nearly all periodical modes were stable at the tested airspeeds and flight altitudes. The stability requirements were defined as constraints for the optimization algorithm. The constraints were enforced using an indirect method which imposes a penalty for any cases beyond the boundaries. The penalty does not ensure that constraints are beyond the limits; thus, in a few cases, the boundaries were slightly crossed, but remained within an expected tolerance. In the following, the dynamic mode characteristics are discussed in detail with respect to the airworthiness criteria. The results were obtained for aftmost CG location at the maximum weight.

With respect to the CS-23 airworthiness requirements [74], the phugoid mode was noted to not be particularly strong. In Figure 12, the phugoid mode damping ratio against the MIL criteria has been shown [75,76]. The damping characteristics satisfy the Level 1 acceptance for all range of airspeeds from stall to the maximum cruise speed. Figure 13 shows time to double for the phugoid amplitude. The negative values show that the oscillations are stable.

Figure 14 presents the short period undamped natural frequency versus acceleration sensitivity. As can be seen, the results were inside the envelope of the first acceptance level according to MIL criteria [74]. Figure 15 further shows the damping ratio level for the short period mode, which nicely satisfies the first acceptance level of the CS-23 requirement to be strongly damped. Reprinted from [74].

According to CS-23.181, the Dutch roll requirements were given precisely as: “Any combined lateral-directional oscillations (“Dutch roll”) occurring between the stalling speed and the maximum allowable speed appropriate to the configuration of the aeroplane must be damped to 1/10 amplitude in 7 cycles…”. Figure 16 shows that this criterion has also satisfied the design.

Figure 17 shows the time to double the roll mode amplitude. With respect to Figure 17, the roll mode was strongly damped, which is typical for the most airplanes. The spiral mode characteristics are presented in Figure 18. As can be seen, values for time to double show that the spiral mode is also strongly damped according to CS-23 and MIL-F-8785-C. [72,74]. Complying with both regulations, the spiral mode falls within level 1 of requirements in most cases. For very low speeds near stall speed and at higher altitudes, the spiral mode level may decrease to the level of 2 or 3. At all points, however, the results remain within the acceptable range of the criteria.

5. Conclusions

In conclusion, the proposed multidisciplinary possibilistic approach successfully provides the capability to size the empennage package of a multi-engine propeller-driven aircraft by altering the tail geometry in the early design stages where several uncertainties are present due to low fidelity estimations. In this work, four disciplines were considered including aerodynamics, stability and control, propulsion and weight and balance. Considering various design criteria, including longitudinal/lateral/directional trim, stability characteristics, and other requirements, the optimized tail configuration provided the design with the desired characteristics to ensure that airworthiness requirements were met and predefined constraints were complied with while offering an optimal form. Results showed that the aircraft empennage sizing using both scenarios of mass minimization and the drag optimization in trimmed flight produced very similar configurations. The proposed approach may be used to improve the preliminary design process of multi-engine propeller driven light aircraft.