Instead, the length and number of the hybridised mission phases affected the shape and presence of a Pareto front. The results clearly showed that, while the longest mission phases contributed the most to the reduction of fuel consumption, if most mission phases were hybridised, the NOx emissions would go down together with an increment in the overall propulsive efficiency. Moreover, the maximum take-off mass limit affects how much trade-off is present between the objectives. These phenomena can be explained in the limited amount of available electric energy due to the limited amount of battery mass. If this energy is distributed over more mission phases, on average, less gas turbine power is used, and therefore, less NOx is emitted as it runs at a lower combustion temperature. However, the DoH would be lower on the longer mission phases, requiring most fuel burned. Indeed, if the electric energy is used only on the longest mission phase, cruise, the mission fuel consumption is minimised.

Lastly, the maximum take-off mass determines how much battery mass is distributed over the mission phases. While this might suggest increasing this limit indefinitely, the increment in the battery mass and the higher amount of electric power to manage would greatly impact the aircraft's structural and systems mass, increasing the aircraft's cost and complexity and ultimately reducing the fuel efficiency benefits.

This aspect was not included in this study but will be included in future work. Furthermore, from this study, the thermal management aspect of the hybrid-electric propulsion system was ignored, which is one of the current challenges in the design and production of HEA. Future studies will introduce the sizing of the thermal management system (TMS) from the electric system lost energy and its effects on aircraft mass and available power.

Finally, a critical aspect in the study and development of HE and Electric propulsion is the uncertainty of future technology. In this study, a fixed value of battery energy density was selected for a 2050 scenario. However, future studies will introduce an uncertainty into the presented methodology for studying different technological scenarios and their impact on the design parameters and requirements at the same time.

**Author Contributions:** Conceptualization, T.K. and A.S.; methodology, A.S., T.K. and H.B.E.; software, A.S. and H.B.E.; validation, A.S. and H.B.E.; formal analysis, A.S. and T.K.; investigation, A.S. and H.B.E.; resources, T.K.; data curation, A.S.; writing—original draft preparation, A.S. and H.B.E.; writing—review and editing, A.S., H.B.E., B.Z. and T.K.; visualization, A.S. and T.K.; supervision, T.K. and P.L.; project administration, T.K.; funding acquisition, T.K. and P.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This project has received funding from the European Union's Horizon 2020 Research and Innovation programme under Grant Agreement No 875551.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Results not available, raw data are cited in text.

**Acknowledgments:** We thank all the researchers of Cranfield University working in the FUTPRINT50 project. We would also like to thank the reviewers for their constructive comments which helped to improve the quality of the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

#### **Appendix A**

**Figure A1.** Objective space of cruise-only hybridised strategies.

**Figure A2.** Objective space of Climb and cruise hybridisation strategies.

**Figure A3.** Objective space of Uniform strategies with different hybridised mission phases.

**Figure A4.** Objective space of Linear strategies with different hybridised mission phases.

**Table A1.** Detailed results of the optimal points where Climb and Cruise are hybridised. (Section 5.2).

**Figure A5.** Optimal energy management strategies of points in Figure 12a.

**Figure A6.** Test 4 optimal energy management strategies.

**Figure A7.** Optimal energy management strategies for different take-off mass constraint levels.

#### **References**

## *Article* **Wing Structural Model for Overall Aircraft Design of Distributed Electric Propulsion General Aviation and Regional Aircraft**

**Raquel Alonso Castilla †, Florent Lutz \*,†, Joël Jézégou and Emmanuel Bénard**

ISAE-SUPAERO, University of Toulouse, 31000 Toulouse, France;

<EMAIL_ADDRESS>(R.A.C<EMAIL_ADDRESS>(J.J<EMAIL_ADDRESS>(E.B.)

**\*** Correspondence<EMAIL_ADDRESS>
† These authors contributed equally to this work.

**Abstract:** In the context of reducing the environmental footprint of tomorrow's aviation, Distributed Electric Propulsion (DEP) has become an increasingly interesting concept. With the strong coupling between disciplines that this technology brings forth, multiple benefits are expected for the overall aircraft design. These interests have been observed not only in the aerodynamic properties of the aircraft but also in the structural design. However, current statistical models used in conceptual design have shown limitations regarding the benefits and challenges coming from these new design trends. As for other methods, they are either not adapted for use in a conceptual design phase or do not cover CS-23 category aircraft. This paper details a semi-analytical methodology compliant with the performance-based certification criteria presented by the European Union Aviation Safety Agency (EASA) to predict the structural mass breakdown of a wing. This makes the method applicable to any aircraft regulated by EASA CS-23. Results have been validated with the conventional twin-engine aircraft Beechcraft 76, the innovative NASA X-57 Maxwell concept using DEP, and the commuter aircraft Beechcraft 1900.

**Keywords:** aircraft design; distributed propulsion; general aviation; regional aircraft; structural mass; wing design; aviation; MDAO; FAST-OAD; innovative aircraft architecture

#### **1. Introduction**

The goals of reducing the environmental footprint and decreasing the noise levels determine the trends of future aviation. In this sense, distributed electric propulsion (DEP) aircraft represent a potential key approach for meeting the overall objectives without compromising costs and safety requirements. The concept of blowing air at high speed over the wing generates a relative speed of the air with respect to the aircraft. As a result, increasing the number of propellers demands less thrust to reach the minimum speed during take-off, and it increases the lift coefficient and requires a smaller wing surface [1]. From a structural point of view, the incorporation of DEP along the span leads to the alleviation of the bending moment and the reduction of the wing structural mass [2]. Therefore, DEP affects the mass breakdown of the wing, and its assessment is one of the main challenges to be faced in the weight analysis.

The design of a new aircraft, either innovative or conventional, generally starts by a conceptual design (see Figure 1) in which the preliminary assessment of a variety of possible configurations is studied. This phase, if conducted properly, allows us to rule out unfeasible layouts and identify key features [3]. This stage is but one of the iterations of the aircraft design process and is likely to be conducted several times for each tested configuration. Additionally, many new aircraft concepts are taking shape, which creates the need for a way to quickly and reliably assess their likeliness to fly in conceptual and preliminary design steps. FAST-OAD-GA is an overall aircraft design (OAD) code that

**Citation:** Alonso Castilla, R.; Lutz, F.; Jézégou, J.; Bénard, E. Wing Structural Model for Overall Aircraft Design of Distributed Electric Propulsion General Aviation and Regional Aircraft. *Aerospace* **2022**, *9*, 5. https:// doi.org/10.3390/aerospace9010005

Academic Editors: Spiros Pantelakis, Andreas Strohmayer and Liberata Guadagno

Received: 26 November 2021 Accepted: 16 December 2021 Published: 22 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

aims to fulfill the need for general aviation (GA) aircraft. As part of the undertaking of adapting FAST-OAD-GA for novel configuration, the necessity for an improved wing mass estimation method was identified.

**Figure 1.** Aircraft design process for a typical GA aircraft. Reprinted from Ref [3].

A literature review has thus been conducted to assess the relevancy of existing wing mass estimation methods for conceptual and preliminary aircraft sizing. Three types of methods have been identified based on their fidelity level: empirical, physics-based, and finite-element methods.

Empirical models such as the one presented by Raymer [4] and Nicolai [5] are particularly well suited for the conceptual design stage and are commonly used. Their computation time is virtually instantaneous, and the level of detail required to apply the formulae matches what is usually found at this stage of aircraft design. Such methods are, however, strongly dependent on the database used to establish them, and with the lack of information on innovative designs, they would fail to capture the aforementioned effects [6]. More refined methods such as the one presented in NASA's FLOPS [7] are capable of encapsulating the design features of unconventional aircraft configurations. This is the case for high-aspect-ratio wings, strut-braced wings, or aeroelastic tailoring. Nevertheless, the aerodynamic loads and the engine relief effect are only used in part of the computation of the wing weight, and they must be specified by the user.

On the opposite end of the complexity spectrum, one finds finite element methods. Some of them, such as the one presented in [8], have been shown to be able to estimate the wing mass with satisfactory results. However, such methods present computational times that are not acceptable for conceptual design phase. Furthermore, the required inputs and precision level of those models are usually not known at this stage of the design. The final limiting step is the complexity of such codes or software that makes their integration efficient in an overall aircraft design software complex.

Finally, existing analytical and physics-based methods have been reviewed. These combine the ease of integration with flexibility and relatively low computational time. Nonetheless, they are often proprietary software or embedded in an already existing OAD tool. Other approaches such as the one developed in [9] or [10] present the appropriate level of fidelity and are described with enough precision that they can be implemented to estimate a basic mass breakdown of the wing. However, these methodologies do not consider the loads cases and structural requirements prescribed by EASA CS-23, and they are more oriented towards commercial transport aircraft. Nonetheless, these have served as the basis to build the methodology presented in this paper.

This work focuses on developing and validating a semi-analytical approach to predict the structural mass breakdown of a wing for a usage in a preliminary aircraft design tool. The choice of a semi-analytical model was made to combine the robustness of the method, the ease of the integration, and the low computational time. This model shall be applicable to aeroplanes in EASA CS-23 category and be adapted to both conventional and distributed propulsion airplanes. The methodology developed by [10] was adapted to consider both the relevant CS-23 loads certification requirements and the effects of DEP. After a validation of the developed method on three aircraft covered by the CS-23 requirements, parametric studies were conducted to investigate different design choices. This methodology was then implemented in FAST-OAD-GA, which brought its own set of challenges, chief among which was the concern for the computational time and opportunities such as benefiting from the software framework.

#### **2. Materials and Methods**

#### *2.1. FAST-OAD-GA (Future Aircraft Sizing Tool—Overall Aircraft Design—General Aviation)*

FAST-OAD-GA is an open-source plugin of FAST-OAD, an open-source aircraft sizing and optimization tool developed by ISAE-SUPAERO and ONERA based on the OpenM-DAO framework [11]. The links to the code are available in the Supplementary Material section. A FAST-OAD run consists in the assembly of several OpenMDAO components with inputs and outputs that represents the application of an aircraft sizing methodology. Moreover, FAST-OAD allows us to define specific problems by using a configuration file. This file provides an easy way to switch between models or to remove/add some modules. Coupled with OpenMDAO's multidisciplinary design analysis capabilities, it allows for a preliminary estimation of the viability of possible aircraft configurations. Currently, FAST-OAD-GA is bundled with statistical models for mass estimation targeted at conventional airplane architectures. However, in the extension of this software to cover innovative aircraft designs, this was targeted as a limiting factor, which prompted the creation of the presented methodology.

**Figure 2.** Representation of the typical FAST-OAD-GA process as an N2 diagram.

A representation of the typical FAST-OAD-GA sizing process is presented in Figure 2 in the format of an N2 diagram. Each of the colored lines on the bottom-left side of the diagonal represents a loop in the process, while the one on the upper-right are forward computation. Henceforth, the outermost loop, i.e., the loop on the wing area, will be referred to as an iteration. Figure 2 shows that each of the modules (gray blocks in the main diagonal) can be used several times in each iteration. More specifically, the weight module that computes the overall empty weight (OEW) is inside 3 loops in each iteration, meaning it is run at least three times. The OEW is currently computed by summation of the weights of every component of the aircraft broken down into categories adapted from [12]. The wing corresponds to one of the main components that contributes to the OEW. The work presented in this paper consists in replacing the current way of computing the wing weight by the semi-analytical methodology mentioned in the introduction to improve the precision of the sizing process while compromising the computational time. Figure 2 also illustrates that this model will be used in a complex environment with a methodology at various levels of fidelity, making it unnecessary to have a model with a level of detail in the inputs greater than what the rest of the process can provide, which was one of the reasons why a finite element method was not retained.

Furthermore, Figure 2 shows that many modules are run before the computation of the mass breakdown. This means that since the outputs of some components can be the inputs for others, it is easy to reuse the results of previous models without requiring their computational time. An example of a module in this situation is the computation of the aerodynamic performance of the aircraft, which contains the computation of the lift distribution along the span.

In order to be able to use a real lift distribution as part of the computation of the aerodynamics loads, FAST-OAD-GA relies on OpenVSP through an interface that allows an automated processing of the computation and results. It is an open-source software released under NASA Open Source Agreement as a parametric geometry tool. One of OpenVSP's tools is VSPAERO which allows aerodynamic analysis based on either a vortex lattice or a panel method. OpenVSP is well suited in early aircraft design stage because of its low computational time and ability to run various types of analyses. Another advantage of using this software is that it allows the simulation of the wing-propeller interaction due to the induced velocity through the use of an actuator disk model [13], which is a key feature to capture the effect of DEP in the aerodynamic loads.

#### *2.2. Wing Mass Model*

For the purpose of initial design, the wing is assumed to be trapezoidal, with N engines distributed along the wingspan and divided into a primary and secondary structure. The leading and trailing edge devices such as ailerons, flaps, and slats form the secondary structure of the wing. The primary structure is modeled as a wing-box, represented in Figure 3 with a miscellaneous weight to account for the joints, cut-outs, and connections. The wing-box structure is then broken into different components. Each of these components is assumed to withstand only one type of stress:

**Figure 3.** Wing-box model. Reprinted from Ref [3].

This set of assumptions allows a simplified analytical estimation of their individual mass based on É. Roux's model [10], which is only valid for isotropic materials. The use of composites is accounted for in the miscellaneous weight through the FCOMP factor [7], ranging from 0 for no composites to 1 for its maximum utilization.

The choice was made to ultimately not consider buckling for the sizing of the ribs. Indeed, this method was elaborated to be used as part of a conceptual design tool and does not substitute a full load analysis, which means we had to rule out some failure scenarios such as buckling.

In addition to the previous model, Sui An [14] shows in her studies the impact of flutter on the wing mass when using distributed propulsion. She proves that the flutter constraint increases the optimum mass by 6–14%. Therefore, the current method will approximate such effects in the applicable cases through a mean coefficient of 1.10.

#### *2.3. Methodology Description*

The following methodology is applicable to airplanes regulated by EASA CS-23 certification specification. Requirements from CS-23 at amendment 4 [15] are considered throughout this paper as an acceptable means of compliance to amendment 5, the latest applicable issue written in a performance-based language. The goal is to compute for each part of the wing the quantity of material necessary to support the stresses induced by the load cases prescribed in the certification requirements. Figure 4 describes the iterative process developed for the evaluation of the mass breakdown of the wing, which begins with an initial estimation of the total wing weight for the computation of the structural loads and then follows the steps presented hereafter.

**Figure 4.** Flow diagram of the method.

1. Load factors: CS 23.321 and CS 23.333 [15] state that the aircraft must comply with the load requirements for any combination of load factor and airspeed within the limits of the flight envelope. This envelope results from the application of the maneuver and gust criteria defined for every weight and flight altitude computed according to CS 23.335, CS 23.337, and CS 23.341. To simplify the process, two constraining weight cases for the wing structure at the typical cruise altitude were selected: the aircraft

at MTOW and with minimum fuel in the wings [15]. The choice of considering only the cruise altitude and cruise equivalent airspeed *Vc* was made since these conditions lead to the highest gust load factor. In addition to that, according to CS 23.305, the wing should also demonstrate compliance to the strength requirements at ultimate loads. These are calculated by applying a safety factor to the maximum/minimum load factors obtained for the aforementioned conditions.

2. Lift distribution: The computation of the lift distribution was separated into two parts: the clean lift, which is the lift without accounting for engine wake and the slipstream lift, meaning the additional lift generated by the engine wake. The clean lift distribution is obtained from OpenVSP and was assumed to be of constant shape along the wingspan. It is then scaled to match the lift coefficient required to sustain the extreme load factors. In the computation of the clean lift, the tail down force necessary to maintain the aircraft level was also considered, leading to a 5% increase assumed to balance the weight in cruise (L = W) [12]. The slipstream lift distribution is computed using OpenVSP through a simple actuator-disk engine model to capture the effects of DEP on the aerodynamic loads. For the computation of the slipstream effect it was considered that two types of propellers could be installed on the wing: the High-Lift Propeller (HLP) [16] to provide extra lift at low speed and the Wingtip Propeller (WTP) to provide thrust during cruise and reduce induced drag. In cruise, the HLPs are folded, but at the climb cruise transition, they are still activated and working close to their maximum operating speed. At this time of flight, the potential risk of whirl flutter require setting the torque at the minimum level with a consequent reduction of power [17]. These few seconds of blown wing under cruise conditions are taken into account in the wing sizing process under a steady approximation to ensure conservative results. This lift distribution is then multiplied by the chord at the corresponding station along the span to obtain the aerodynamic loads at each point of the wing. The results of the computation of the lift load for the X-57 are shown in Figure 5.

**Figure 5.** Lift load distribution along the dimensionless wing span.

the density of the material *ρ<sup>m</sup>* multiplied by the volume of the part. The following computations will be done for the two previously mentioned aircraft masses: MTOW and minimum fuel in the wing. Consequently, two values of the mass will be obtained for each element, with the maximum being retained. Hereunder, a non-dimensional reference coordinate system will be used, where the *Y* axis (*Y* = *y*/(*b*/2)) is defined positive to the right side of the aircraft, with *b* the span.

• Web mass in shear *Mweb*: the spar web is a vertical panel designed to react, under assumptions of the presented methodology, to all the shear force of the wing *T*(*y*). CS 23.305 stress criteria is applied, which states that the stress on the web must be lower than the maximum shear stress of the material |*τweb*| ≤ *τmax*. So, the mass of the web for the whole wing is calculated integrating its section area *Sweb* along its elastic axis *η* which sweep angle is *φe*:

$$M\_{\rm uvb} = 2\rho\_m \frac{b/2}{\cos \phi\_\ell} \int\_0^1 S\_{\rm uvb} dY \qquad ; \qquad S\_{\rm uvb} = \frac{|T\_{\rm max}|}{\tau\_{\rm max}} \tag{1}$$

• Spar caps mass in bending *Mcaps*: the lower and upper spar caps of the main spar react to the bending moment *Mf*(*y*). Similarly to the web sizing, the CS 23.305 stress criteria imposes that the stress on the spar caps cannot surpass the maximum normal stress of the material, *σf l* <sup>≤</sup> <sup>|</sup>*σmax*|, neither in traction (*σmaxt* > 0) nor in compression (*σmaxc* < 0). In general, the resultant shear force of the wing is in the upwards direction. Hence, the upper spar cap will work under compression, while the lower spar cap will be subjected to traction. The consequence is that the upper spar cap will face higher risk of failure under buckling than the lower spar cap [10]. With those considerations, the mass of both spar caps can be calculated with Equation (2), where *h*(*Y*) is the distance between the center of gravity of the spar caps retrieved from [10] and *SE*, *SI* are the cross-section areas of the upper and lower spar caps, respectively.

$$M\_{\rm caps} = 2\rho\_m \frac{b/2}{\cos \phi\_c} \left[ \int\_0^1 S\_E dY + \int\_0^1 S\_I dY \right] \tag{2}$$

$$\begin{cases} S\_I = \frac{1}{h(Y)} \max\left(\frac{Mf(y)^+}{\sigma \max\_l}, \frac{Mf(y)^-}{\sigma \max\_c}\right) \\\\ S\_E = \frac{1}{h(Y)} \max\left(\frac{-Mf(y)^+}{\sigma \max\_c}, \frac{-Mf(y)^-}{\sigma \max\_l}\right) \end{cases} \tag{3}$$

• Skin mass in torsion *Mskin*: The skin is a thin panel of a thickness assumed constant that provides the necessary strength to bear the torsion moment. This assumption makes sense from a manufacturing point of view, where varying the thickness for wings of this size represents additional manufacturing complexity that is not desirable for general aviation manufacturers. Besides having a varying skin thickness did not significantly improve the results in [10], hence our choice to keep the wing skin thickness constant. The regulation imposes two sizing criteria: CS 23.455 that limits the deformation of the wing in rolling conditions and one that corresponds to the maximum admissible stress at ultimate loads stated in the CS 23.305. Hence, the skin mass will be determined by the most demanding requirement. In previous versions of the model [18], both criteria were computed for various aircraft. The results revealed that the deformation dimensioning also meets the stress requirements. Consequently, only the CS 23.455 rule was implemented in the current edition of the code. Then, as presented in Equation (4), the skin mass is obtained by integrating the product of the skin thickness, *e* and the perimeter of the wing-box, *l*. The objective of this criterion is to guarantee the efficiency of the ailerons control by limiting the roll speed at the design airspeeds *pA*, *pC*, *pD*. When ailerons are deflected, this norm

restricts the wing deformation that changes the angle of attack Δ*α*. This Δ*α* can be approximated by the variation of the twist angle due to the torsion moment produced by the ailerons. The skin weight will then be computed based on the necessary thickness to react that torsion moment. For further detail on the calculation of *e*, *pA*, *pC*, *pD* the user is encouraged to see [10] (pp. 84–104).

$$M\_{skin} = 2\rho\_m \int\_{\frac{R\_{fus}}{\text{cm}\_{\text{f}}}}^{\frac{k/2}{\text{cm}\_{\text{f}}}} c(\eta) l d\eta \qquad ; \qquad \begin{cases} p\_C \ge p\_A \\ 3p\_D \ge p\_A \end{cases} \tag{4}$$

The mass of the ribs *Mribs* and miscellaneous *Mmiscellaneous* items are obtained empirically according to the expressions proposed by É. Roux [10] and NASA FLOPS [7], respectively. As for the secondary structure, Torenbeek [2] states that it contributes to around a 25% of the total mass of the wing. Therefore, its weight will be empirically obtained with this proportion. Finally, the total mass of the wing will be the sum of all its components:

$$M\_{\text{primary}} = M\_{\text{web}} + M\_{\text{caps}} + M\_{\text{skin}} + M\_{\text{ribs}} + M\_{\text{misccllanocons}}\tag{5}$$

$$M\_{\text{win}} = M\_{\text{primary}} + M\_{\text{secondary}}$$

5. Iteration: the method finishes by comparing the computed mass of the wing with the one initially estimated. If these values are different, the tool iterates on the solution with the computed value until the error is lower than 0.1%. This ensures the obtained wing weight is close enough to the structural weight used to compute relief effects.

#### *2.4. Methodology Validation*

Validation of the proposed methodology and related models is performed through an assessment for three aircraft use-cases: the Beechcraft 76 as a conventional twin-engine airplane, the innovative NASA Maxwell X-57 as a DEP concept, and the Beechcraft 1900 as a commuter aircraft. In this way, the robustness and flexibility of the method can be proved through a realistic range of weights and sizes of aircraft. Input data used for those use-cases are detailed in Appendix A. The structural safety factor for the Maxwell X-57 was set to 1.8 according to the airworthiness approach set for the Mod III version of the aircraft [19].

#### *2.5. Parametric Analysis*

Two parametric studies were conducted to show the importance of considering the impact of DEP and engine placement on the wing sizing. For the analysis described in Sections 2.5.1 and 2.5.2, the aircraft we will use as a reference is the Maxwell X-57. For the study in Section 2.5.3, we will use the Beechcraft 76. As previously introduced, the weight of the engines and the wing-propeller interaction modify the structural loads with subsequent effects on wing mass. For a better comprehension, each of these effects is presented and analysed independently.

#### 2.5.1. Effect of the Engines Weight Distribution

The investigation presented in this section aims to analyse exclusively the engines mass effect on the wing sizing. The complete study of the effects of a change in the propulsion system is an Overall Aircraft Design problematic that is not addressed in this paper. The results of this study, however, can be seen as the first iteration of the global OAD process. For such purpose, the quantity of high-lift propellers per semi-wing has been varied from 1 to 6, and the following assumptions have been retained:

• The high-lift engines are evenly spaced along the wingspan and it is assumed that the distance between the last HLP and the WTP was greater than the distance between two HLP. The position used for this study is detailed in Table 1.

**Table 1.** HLP position considered in the parametric study on the effect of HLP weight (in % of the semi-span) per half-wing.

#### 2.5.2. Effect of Slipstream Consideration

The second parametric study shows the sensitivity of the model to the slipstream effect in cruise conditions. This analysis is made to evaluate the impact of the choice made in this paper to consider that HLP were still active in the cruise conditions. The original configuration of the X-57 Maxwell with 14 engines is used. In order to achieve this objective, the mass of the wing has been obtained with and without adding the slipstream lift.

#### 2.5.3. Effect of Engine Placement

The third parametric study aims at analysing the influence of the position of the engine on the overall wing mass for a conventional design. The baseline configuration for this study is the Beechcraft 76. The position of the engine varies between the root and the tip of the wing. For a correct interpretation of the conclusions of this study, it is necessary to mention that, for the computation where the engine is located closer to the fuselage centerline than the reference position, results are physically irrelevant. The baseline engine location is indeed already a lower limit based on the propeller clearance, which can be specified as a design constraint. The results will nevertheless be included as an indicative basis.

#### **3. Results**

#### *3.1. Methodology Validation*

For the purpose of validation of the methodology, the outputs of the use-cases have been evaluated against data available for Beechcraft 76, Beechcraft 1900 and NASA X-57 Maxwell with DEP. For the latter, the weight of the wing was directly available in the literature [20]. For the Beechcraft 76 and the Beechcraft 1900, as no real values are publicly available, the comparison is made with a traditional statistical method [4,5,12]. Said statistical methods are applicable to normal category aircraft and commuter aircraft and were established based on database containing aircraft which technological level is close to that of the Beechcraft 76 and 1900. They are written in the Appendix B.

A summary of the results of the validation study is presented in Table 2.

For the Beechcraft 76 use-case, the deviation of the wing mass obtained with the presented semi-analytical model compared to the reference value is less than 1%. The actual value of the wing mass is unknown; however, this result shows that the newly established methodology gives the same results as a classical empirical formula. This prompted the authors to compare the result with an empirical methodology [5]. In this case, the difference computed was of around 6%, which confirms that no significant precision is lost when using the semi-analytical approach.

**Table 2.** Outputs of the presented methodology applied on the Beechcraft 76, the X-57 Maxwell, and the Beechcraft 1900.

The same conclusion can be drawn from the validation of the method against the commuter Beechcraft 1900. The empirical formula used to compare the results obtained for this aircraft is selected from [12], which suggests the use of commercial transport aircraft models as the most appropriate choice to describe turbo-propeller aircraft with a maximum take-off mass over 5670 kg, as is the case for this aircraft.

For the NASA X-57 Maxwell use-case, the empirical formula from the literature gives results with an error of around 20%. This can partly be explained by the fact that the correlation was established on aircraft whose wing loading was fairly low because of the flaps technology used. Given that DEP allows one to reach a lift coefficient of over 4.0 in certain configurations [21], this leads to small wings but with heavy structure. By contrast, the semi-analytical methodology provides results that are much closer to the real wing mass. The reason behind this is that it actually takes into account the structure required to withstand the loads compliant with CS-23 requirements.

In addition to the total wing mass, the presented methodology also provides a further breakdown of the wing structural components shown in Table 3.

From this table, a noticeable remark can be formulated about the change in the proportional importance of the spar caps on the total wing weight. Indeed, for the X-57 Maxwell they take a much larger share. This can be explained by the increase in the wing aspect ratio (AR), which is equal to 15 for X-57 Maxwell, compared to 8 for the Beechcraft 76 and 10.8 for the Beechcraft 1900. The AR is one of the main factors in the computation of the bending moment [3], and since the bending moment is the defining factor for the spar caps computation, this explains such an increase.

**Table 3.** Detailed mass breakdown (absolute and share of total wing weight) of the Beechcraft 76, the X-57 Maxwell, and the Beechcraft 1900.

#### *3.2. Parametric Studies*

3.2.1. Effect of the Engines Weight Distribution

The goal of this study is to see if the presented methodology was capable of replicating the relief in the bending moment that comes from the distribution of engines along the wing. The results are presented in Table 4 and plotted in Figure 6, where the dashed lined is the reference weight for the X-57 Maxwell.

**Table 4.** Mass breakdown [kg] of the X-57 Maxwell wing for different number of engines on the wing.

The results show an alleviation effect on the total wing mass of around 3.5 kg for an increasing number of engines. It can be observed that the webs, spar caps, and secondary structure masses are the only values that change. As explained before, the ribs and miscellaneous weights are computed through empirical methods. These techniques do not take into account the loads or the number of engines, which justifies that no change is observed. As for the skin, its mass only depends on the characteristic speed of the aircraft. Therefore, considering this study's assumptions of constant MTOW and no slipstream effects, the skin weight does not vary with the number of engines.

Based on these observations and the hypothesis made regarding the secondary structure mass, which will account for 25% of the total wing mass, the weight variation is due to a change in the primary structure. This is consistent, since the primary mass is mainly obtained as the weight of the material required to withstand the loads. Since the weight of the engines relieves the aerodynamic forces to react to the structure of the wing, a difference in the engine weight and position will impact the structural weight.

**Figure 6.** Effect of the number of engines on the overall wing mass.

For a deeper analysis into the reason for this lightening of the primary structure mass, Figure 7 plots the contribution of each structural component to the primary structure mass. This bar chart was plotted based on the results obtained for the two limit cases with only one HLP and with six HLP.

**Figure 7.** Effect of the number of engines on the primary structure mass.

As seen in Table 4 and Figure 7, the change in the primary structure mass comes from spar caps and web. In particular, the shear force seen by the web and the bending moment borne by the spar caps are lower. Thus, the structures will be lighter. Similar effects were observed in the literature [22], albeit on bigger aircraft, which reinforces the idea that this model can capture the structural benefits that distributed propulsion can bring forth.

#### 3.2.2. Effect of Slipstream Consideration

This study investigates the impacts of the slipstream effects on the wing sizing. The results are presented in Table 5.

**Table 5.** Mass breakdown (kg) of the X-57 Maxwell wing with and without slipstream effects considered.

The results reveal a sensitivity of almost 5 kg ( 3%), resulting in a heavier wing when the induced velocity of the HLP is considered. The air blown by the engines on the wing increases the aerodynamic loads that the wing structure reacts to. Consequently, reinforced and heavier web and spar caps are required.

On the other hand, the current analysis highlights an important point. Depending on the power considered for the HLPs, the slipstream can heavily constrain the wing sizing. In the current approach, a minimum thrust is set, since the study is performed immediately at the end of the climbing phase at *Vc* speed. However, in more detailed design stages, the take-off phase will have to be analyzed with full engine power. In this case, the slipstream effect will significantly exceed the actual 3%, challenging the choice of using *Vc* as the sizing point. The point of design was chosen based on the maximum load factor obtained within the maneuver and gust envelopes. However, for commuter-category aircraft, the

certification specification requires us to also consider the point of maximum gust from the gust envelope [15].

The V-n diagram we obtained in our study on the Beechcraft 1900 is shown in Figure 8, in which point A corresponds to the point we considered and point B corresponds to the maneuvering point. At this point, even if the load factor is no greater than at cruise, the speed is lower and thus the slipstream effects are greater. As a result, this might lead to a less homogeneous aerodynamic loads and a need for stronger primary structure.

#### *3.3. Effect of Engine Placement*

This study focuses on the effect of changing the location of a sole engine along the wing span for a twin-engine aircraft. The reference case is the Beechcraft 76. The results are presented in Figure 9.

**Figure 9.** Effect of the location of the engine on the wing and caps mass.

The results show that the further the engine is from the fuselage, the lighter the spar caps are and thus the lighter the wing mass. This echoes results presented in Section 3.2.1 for the study of distributing the propulsion over the whole wing. The amplitude of the changes, however, is much more significant, around 16% of the baseline weight. These results lead to the same conclusion developed in [22]: moving the thrust generation on the outer part can have a positive impact on the wings' structural weight. Nevertheless, the last points of the figure may be too optimistic, as the model currently does not take into account any structural reinforcement that would be needed to withstand such a mass at the wingtip. It also omits the case of hard landing, which might constrain the stiffness of the wing.

Placing the engine outboard can have significant positive or negative impacts on the overall aircraft design. On one hand, studies on the aerodynamic effect of placing an engine at the tip have shown that clear reduction of the induced drag have been obtained [21,23]. On the other hand, the closer the engines are to the wingtips, the higher the roll inertia of the aircraft will be. This will result in either a decreased roll rate or in larger aileron, which means more complexity and weight. It can also affect the sizing of the vertical tail plane in case of an engine failure during low speed conditions. A rough estimation of the vertical area needed based on the one- engine-inoperative certification requirements was made using FAST-OAD-GA. The results show that from the baseline configuration to the case where the engine is at the tip, the required area is more than doubled. This illustrates that while wing weight savings may be obtained, the overall empty weight might increase, highlighting the importance of coupling this method with an overall aircraft design process.

#### *3.4. Computational Time*

One of the main objectives of the newly developed methodology was that, once implemented in an overall aircraft sizing tool, its computational time was to remain acceptable for conceptual and preliminary design. Numerous experiments were conducted on each use-cases configuration, and the time to converge was registered. The average results presented in Table 6 were obtained.

**Table 6.** Comparison of the computational time between the newly developed methodology and statistical formula.

As expected, the newly developed methodology takes more times than the statistical formula suggests. However, the time is fairly low compared to some other modules of FAST-OAD-GA (mainly propeller performances and aerodynamic computations). With the default sizing process, the computation of the OEW (where the developed method is used) is entered three times in each iteration. Knowing that it takes 10 to 15 iterations to converge, the new model only adds around 7 s to the overall process. Compared to the 5–10 min needed for the whole FAST-OAD-GA process to converge, the time added by the new methodology is deemed acceptable. It is also worth noting that there is a strong dependency of the computational time on the initial estimated wing mass. During all those studies, the initial point was taken as the result provided by a statistical formula [4]. This partly explains why the process was much faster on a traditional aircraft for which both results are fairly close. By using the real wing mass of the X-57 Maxwell as the initial value, the time necessary to converge on the wing mass was decreased to about 120 ms. There is, however, a difficulty in repeatedly predicting an initial wing mass that is close to the results. This implies that for unusual wing configurations, the process might take longer to converge than for a conventional design.

#### **4. Discussion**

The newly developed methodology has proved to be accurate for a conventional twin-engine aircraft, an innovative DEP concept and a commuter aircraft, and it has provided valuable estimations at preliminary aircraft design stage. However, results must be considered carefully. Indeed, even if FAST-OAD-GA allows those values to be changed, a certain number of the assumptions retained in this paper for the structure of the wing (e.g., height of the spar caps, wingbox depth) are based on commonly used values, which may be a source for errors and would require attention from aircraft designers. Some of the computations could also be further improved, namely the ribs' mass, which is currently

obtained with an assumed constant spacing. This could be improved by considering, for instance, that the ribs are sized to prevent the wing from buckling [10,24].

For the skin weight, early results in conventional design seemed to show that the deformation criterion was more constraining than the stress one. However, the use of leading edge propellers placed further from the leading edge may reverse this common tendency. This position has been identified as a major design parameter when considering the sizing of a wing with distributed propulsion [16]. It also plays a major role in the computation of the torsion moment and has an impact on the aerostructural properties of the wing by moving the section's center of gravity forward [25].

On the topic of aero-structural interaction, it was considered that for every aircraft using distributed propulsion, an additional 10% should be added to the wing mass to prevent flutter. This percentage was the mean value of several experiments on different DEP configurations [14]. Further improvements could be added to the presented methodology by running tailored aero-structural analysis to better estimate this coefficient. This, however, would come with an increased computational time and would be an optional module in FAST-OAD-GA.

#### **5. Conclusions**

In this paper, and in the context of the continuous improvement of FAST-OAD-GA, a methodology has been developed to estimate the structural mass breakdown of a wing that while taking into account the effect of distributed propulsion. This was done by adapting a semi-analytical methodology to take into account a realistic load distribution. This new method was specifically designed for aircraft covered by the EASA CS-23 certification specification and was adapted for use in the conceptual and preliminary design stages. Good agreement with actual available aircraft data or previous methodology was shown, and documented effects could be replicated, which justifies its use inside the mass estimation step of the overall aircraft design process. The novelty of this method comes not only from the adaptation of the models to conform to CS-23 and to enable the consideration of DEP, but also from the work that has been done regarding its ease of integration inside an overall aircraft sizing tool. Thus, the present work enables a preliminary consideration of the structural synergies brought forth by DEP. In addition to that and thanks to FAST-OAD-GA's API, this model can be used in any Python-based OAD tools.

**Supplementary Materials:** FAST-OAD is available online at https://github.com/fast-aircraftdesign/FAST-OAD, accessed on 26 November 2021; FAST-OAD-GA is available online at https: //github.com/supaero-aircraft-design/FAST-GA, accessed on 26 November 2021.

**Author Contributions:** Conceptualization E.B., J.J. and R.A.C.; methodology, R.A.C. and F.L.; software, R.A.C. and F.L.; validation, R.A.C. and F.L.; investigation, R.A.C. and F.L.; data curation, F.L.; writing—original draft preparation, R.A.C., F.L., E.B. and J.J.; writing—review and editing, F.L., R.A.C. and J.J.; visualization, F.L.; supervision, J.J. and E.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** The presented work is supported by Daher through the research chair ISAAR "Innovative Solutions for Aviation Architecture & Regulation".

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available in the FAST-OAD-GA repository in Supplementary Materials and in the article.

**Acknowledgments:** The authors acknowledge the support of Daher through the research chair ISAAR "Innovative Solutions for Aviation Architecture & Regulation".

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

#### **Appendix A. Statistical Formula Used for Methodology Validation**

This formula taken from [4] estimates the wing mass based on the wing area *Sw* in ft2, the fuel weight inside the wing *Wf w* in lb, the aspect ratio *A*, the sweep angle Λ at 25% MAC in rad, the dynamic pressure at cruise *q* in lb/ft2, the taper ratio *λ*, the thickness-to-chord ratio t/c, the ultimate load factor *Nz*, and the flight design gross weight W*dg* in lb.

$$\mathcal{W}\_{\text{wing}} = 0.036 S\_w^{0.758} \mathcal{W}\_{fw}^{0.0035} \left(\frac{A}{\cos^2 \Lambda}\right)^{0.6} q^{0.006} \lambda^{0.04} \left(\frac{100t/c}{\cos \Lambda}\right)^{-0.3} (\mathcal{N}\_{\tilde{z}} \mathcal{W}\_{d\tilde{g}})^{0.49} \tag{A1}$$

The formula from [5] computes the wing mass based on the ultimate load factor *Nz*, the flight design gross weight *Wdg* in lb, the aspect ratio *A*, the sweep angle Λ at 25% MAC in rad, the wing area *Sw* in ft2, the taper ratio *λ*, the thickness-to-chord ratio t/c, and the maximum level airspeed at sea level *VH* in KEAS.

$$\mathcal{W}\_{\text{wing}} = 96.948 \left[ \left( \frac{N\_z \mathcal{W}\_{d\_\mathcal{S}}}{10^5} \right)^{0.65} \left( \frac{A}{\cos^2 \Lambda} \right)^{0.57} \left( \frac{S\_W}{100} \right)^{0.61} \left( \frac{1+\lambda}{2t/c} \right)^{0.36} \sqrt{1 + \frac{V\_H}{500}} \right]^{0.993} \tag{A2}$$

The formula from [12] used for the Beechcraft 1900 uses the maximum zero fuel weight *WMZF* in lb, the span b in ft, the sweep angle Λ<sup>50</sup> at 50% MAC in rad, the ultimate load factor *Nz*, the wing area *Sw* in ft2, and the maximum thickness at wing root *tr* in ft.

$$\mathcal{W}\_{\text{wing}} = 0.0017 \mathcal{W}\_{\text{MZF}} \left( \frac{b}{\cos \Lambda\_{\text{50}}} \right)^{0.75} \left[ 1 + \sqrt{\frac{6.3 \cos \Lambda\_{\text{50}}}{b}} \right] \mathcal{N}\_{z}^{0.55} \left( \frac{b \mathcal{S}\_W}{t\_r \mathcal{W}\_{\text{MZF}} \cos \Lambda\_{\text{50}}} \right)^{0.30} \text{ (A.3)}$$

#### **Appendix B. Reference Technical Data for the X-57 Maxwell and the Beechcraft 76**

**Table A1.** Aircraft weights.

**Table A2.** Aircraft geometry.

**Table A3.** Aircraft aerodynamics.

**Table A4.** Aircraft design airspeeds (EAS).

#### **Table A5.** Propulsion.

#### **Table A6.** Mission.

#### **References**

## *Article* **Multi-Fidelity Optimization of a Composite Airliner Wing Subject to Structural and Aeroelastic Constraints**

**Angelos Kafkas 1,\*,†, Spyridon Kilimtzidis 2,†, Athanasios Kotzakolios 2,†, Vassilis Kostopoulos 2,† and George Lampeas 1,†**

**Abstract:** Efficient optimization is a prerequisite to realize the full potential of an aeronautical structure. The success of an optimization framework is predominately influenced by the ability to capture all relevant physics. Furthermore, high computational efficiency allows a greater number of runs during the design optimization process to support decision-making. The efficiency can be improved by the selection of highly optimized algorithms and by reducing the dimensionality of the optimization problem by formulating it using a finite number of significant parameters. A plethora of variable-fidelity tools, dictated by each design stage, are commonly used, ranging from costly high-fidelity to low-cost, low-fidelity methods. Unfortunately, despite rapid solution times, an optimization framework utilizing low-fidelity tools does not necessarily capture the physical problem accurately. At the same time, high-fidelity solution methods incur a very high computational cost. Aiming to bridge the gap and combine the best of both worlds, a multi-fidelity optimization framework was constructed in this research paper. In our approach, the low-fidelity modules and especially the equivalent-plate methodology structural representation, capable of drastically reducing the associated computational time, form the backbone of the optimization framework and a MIDACO optimizer is tasked with providing an initial optimized design. The higher fidelity modules are then employed to explore possible further gains in performance. The developed framework was applied to a benchmark airliner wing. As demonstrated, reasonable mass reduction was obtained for a current state of the art configuration.

**Keywords:** wing design optimization; structural optimization; composite materials; computational fluid dynamics; multi-fidelity optimization

### **1. Introduction**

Ever since the early years of aviation, the design, analysis, and optimization of aircraft structures has been receiving ever-increasing attention from the scientific community and has, therefore, been subjected to extensive research studies. As a result, a wide variety of multi-fidelity computational tools, with application to a plethora of analysis disciplines, has emerged. During an aircraft design process, consisting mainly of the conceptual, preliminary and detailed design stages [1], the various computational tools are designated to a corresponding design stage mainly based on the computational time, the modeling complexity and the ability to capture the physical phenomena present. Early on, low-fidelity models, associated with fast turnaround times, yet incapable of modeling higher order phenomena, are employed, along with empirical knowledge, to steer the design towards optimality. As the design knowledge on the current candidate configuration matures, higher fidelity tools are employed, aiming to replicate the relevant phenomena with greater accuracy, albeit at an elevated and often prohibitive computational cost. Throughout

**Citation:** Kafkas, A.; Kilimtzidis, S.; Kotzakolios, A.; Kostopoulos, V.; Lampeas, G. Multi-Fidelity Optimization of a Composite Airliner Wing Subject to Structural and Aeroelastic Constraints. *Aerospace* **2021**, *8*, 398. https://doi.org/ 10.3390/aerospace8120398

Academic Editors: Spiros Pantelakis, Andreas Strohmayer and Liberata Guadagno

Received: 17 October 2021 Accepted: 12 December 2021 Published: 15 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the stages, all relevant tools are coupled with optimization algorithms in order to gain further knowledge regarding the design space, as well as to obtain the corresponding optimized solution.

A key enabler to the mass reduction trend of aeronautical structures is the rise of composite materials in the aeronautics field. Lighter yet stiffer configurations have emerged with improved static and dynamic aeroelastic response through tailoring of their properties [2]. One of the main challenges, nevertheless, remains the optimization of such structures, since the introduction of composite materials induces further computational and algorithmic complexity. On one hand, ply angles and thicknesses can be expressed as design variables, resulting into a more complex and irregular design space, but with facilitated algorithmic applicability. On the other hand, the use of lamination parameters for optimization purposes has been extensively studied over the past few years. Miki and Sugiyama [3] was among the first to generate optimum designs of laminated composite plates for required in-plane and maximum bending stiffness, buckling strength and natural frequency utilizing lamination parameters. Fukunaga et al. [4] explored the effect of bend-twist coupling on the fundamental frequency of symmetric laminated plates, indicating that this type of coupling reduces the fundamental frequencies. The optimal laminate configuration that maximizes the fundamental frequencies was also obtained. In his work, Liu et al. [5] maximized the buckling load of composite materials panels using flexural parameters and compared its results with a stacking sequence optimization design generated via a genetic algorithm, indicating a close correlation between the two methods. In later works, efforts were directed towards closure of the up-to-then incomplete feasible design space of the lamination parameters [6], and the aeroelastic tailoring of regular and variable stiffness composite materials wings [7,8], as well as the stiffness optimization subject to aeroelastic constraints [9]. Macquart et al. [10], as well as Bordogna et al. [11], extended the capabilities of the state-of-the art lamination parameters optimization algorithms, introducing blending constraints in order to guarantee a certain degree of ply continuity inside a variable stiffness composite wing.

It becomes evident that various disciplines must be involved in the design optimization process of a modern airliner composite wing, the main contributors being the structural and aerodynamic analysis, since they significantly influence the performance and safety of wing structures. Additional disciplines may be included on a case by case basis. Regarding the structural representation, Finite Element Analysis (FEA) models are commonly used at the preliminary and detailed design stages, with reduced order, equivalent beam or plate models aiming to provide greater insight into the candidate configurations, as well as reduce the design space at an early stage of the design process. The equivalent plate methodology (EPM), often preferred due to its inherent facilitated coupling capabilities with concurrent aerodynamic and optimization codes, initially developed by Giles in a series of research papers [12–14], was proven to be a robust and accurate enough numerical tool for the structural analysis of aircraft wings, as the comparison with a 3D FEA wing model showed. Later on, Livne et al. [15] included the EPM in a multidisciplinary structural, aerodynamics and control analysis and optimization framework of a composite aircraft wing, and went on to also include transverse shear effects in the EPM formulation [16]. The large deformation non-linear static and dynamic behavior of wings via the EPM formulation was also studied by Livne and Navarro [17], with good correlation up to relatively large levels of loading being obtained. In their work, Kapania and Liu [18] introduced Legendre polynomials in the EPM in order to avoid the numerical ill-conditioning problem often accompanied by the usage of simple polynomials as trial functions. The accuracy and efficiency of the method was demonstrated via a series of static, as well as free, vibration loading scenarios. Through a series of works, Krishnamurthy and Eldred [19], Krishnamurthy and Tsai [20], Krishnamurthy [21] developed an optimization framework for obtaining an equivalent plate by minimizing an error function accounting for the displacements and frequencies between the two numerical models under consideration. The methodology was initially tested in plates with and without

discrete damage and was later used to evaluate the static and dynamic structural response of a baseline aircraft wing. In each case, an optimum thickness and concentrated mass distribution for the equivalent plate model was obtained. A scaling down methodology was also developed and put into test via a series of numerical tests. The developed methodology indicated that the static response of the wing structure was accurately predicted, with the equivalent plate model being capable of reproducing the first five frequencies of the wing structure within five percent. Over the last few years, Na and Shin [22], along with Henson and Wang [23], further enhanced the current EPM capabilities to account for control surfaces, as well as for tow steered composite wing skins, respectively.

On the other hand, implementation of dedicated computational analysis tools for the relevant disciplines into optimization frameworks is of paramount importance. Pioneering research was conducted by Triplett [24] and Love and Bohlman [25], where one of the earliest multidisciplinary design and optimization tool, aeroelastic Tailoring and Structural Optimization (TSO), was developed. The Rayleigh-Ritz EPM for the structural model, coupled with the Doublet Lattice Method were combined to optimize the thickness distribution and laminate orientations of a fighter wing subject to strength and flutter velocity constraints. Haftka [26] developed an automated procedure for the design and optimization of composite wings subject to strength and flutter constraints. A major conclusion drawn from this work was that the flutter speed may not be a continuous function of the structural stiffness; hence, in optimization under flutter constraint, the constraint should be formulated in terms of other, continuous parameters relative to the flutter phenomenon. Haftka [27] also coupled a lifting-line theory-based aerodynamics model with a wingbox FEA model to perform aeroelastic analysis and optimization of a wing under stress and drag constraints. Weight versus drag trade-off studies for aluminium and composite wings were conducted. It was found that the composite wings were associated with lower mass and drag but their increased flexibility resulted in a larger variation in drag in the Pareto front. One of the earliest aerostructural optimization analysis was performed by Grossman et al. [28], where, utilizing lifting-line aerodynamic models, along with beam equations for the structural representation, the optimum shape, and structural configuration of a sailplane wing was obtained. In a similar but more complex fashion, by means of numerical tools, Grossman et al. [29] optimized a transport aircraft wing, while in parallel developing methodologies for the calculation of the sensitivity derivatives, as well as a sequential approximate optimization module.

While optimization procedures using low-fidelity components for the representation of the structural and aerodynamic problem have to date offered major gains in the overall performance of wing structures, they have inherent limitations. The limitations result from the main assumptions of such methods; for aerodynamics, non-linearities are typically neglected and simplified wetted geometries are used. On the structural side, the assumption of linear behavior, together with simplified FEA models, sees widespread use in optimization. Increasing the fidelity of the aerodynamic and structural solvers can influence the optimization results [30] and represents a step beyond the current state-of-the-art [31,32]. For the coupled aeroelastic problem, linear approaches, both in terms of structure and aerodynamics represented in the frequency domain, are commonly used. Studying the dynamic aerostructural response with such formulations should be limited to linear or weakly non-linear cases due to the underlying assumptions. The dynamic response of wings that operate in the transonic region, as well as cases of boundary layer-shockwave interaction and flutter from flow separation (common in rotor blades), requires a high-fidelity representation of the aerodynamics to be accurately captured [33–35]. Even when the aeroelastic behavior is correctly captured, the increased detail in the unsteady aerodynamic loads can influence the optimum result. In Reference [36], different unsteady aerodynamic methods were used in a multi-fidelity aerostructural optimization of a helicopter rotor, leading to variation in the final optimum blade geometry. Even for the static aeroelastic case, differences due to varying fidelity of the aerodynamics solvers have indeed been pinpointed [37]. Furthermore, the presence of structural non-linearities can additionally have a serious

effect on the aeroelastic behavior of a wing structure [38]. Such non-linearities are expected to become more common in the future due to the tendency to move to ever more efficient wings of very high-aspect ratio combined with highly flexible composite materials.

It becomes evident that the obtained optimized solution corresponds to a specific level of fidelity. Depending on the case, the assumptions in modeling each discipline can be justified. However, they could also represent an oversimplification of the actual physics of the problem and yield a sub-optimal design. The optimized solution can be in doubt due to the level of fidelity with which the various disciplines are approximated regardless of how tight the optimization constraints are set and how rigorous the optimization algorithm becomes. Thus, the fidelity of each component of the optimization framework must be carefully evaluated by taking into account the outstanding properties of each problem or by confirming the validity of the solution using higher fidelity analysis.

In the following brief overview, the emphasis will be mainly upon the impact of increasing the fidelity of the aerodynamic formulation. This is a conscious choice because solution and modeling procedures for even detailed 3D Shell-3D Beam Finite Element representations of the structure have become quite efficient. However, some examples covering the impact of high-fidelity structural modeling are also mentioned below.

Gains from the inclusion of high-fidelity analysis component to the optimization of future unconventional aircraft wing concepts have been reported in literature. In References [39,40], high-fidelity Computational Fluid Dynamics (CFD) was used to provide realistic predictions of engine-pylon-strut interference and viscous drag that was fed to the optimization procedure. Significant reduction in structural mass for a strut braced wing based upon the Common Research Model (CRM) wing planform was obtained by Variyar et al. [41] by using a complete multi-fidelity optimization framework that encompassed both high-fidelity aerodynamics and structural dynamics. For a similar aircraft configuration, the accurate derivation of aerodynamic loads is stressed as a crucial parameter by Variyar et al. [41] and a multi-fidelity approach even in conceptual stage deemed advantageous. Qian and Alonso [42] have included a high-fidelity structural model in an aerostructural optimization procedure to discover gains in range and a reduction in structural mass. Smith et al. [43] mention that including high-fidelity CFD and non-linear structural effects in the aeroelastic modeling of high aspect ratio wings can lead to important differences compared to linear methods. In this particular case, linear methods were found to lead to an overly conservative result for tip displacement and aerodynamic loads.

Unsurprisingly, most recently published studies regarding high-fidelity optimization frameworks deal with innovative wing configurations that are prone to aeroelastic instabilities (high-aspect ratio wings, transonic aircraft) and may also contain structural non-linearities. Quantifying the impact of higher-fidelity optimization to current state-ofthe-art wings has not been as prominent in the literature. Liem et al. [44] have shown that significant benefits from high-fidelity aerostructural optimization can also be realized for current state-of-the-art wings.

The decision on the required analysis capabilities that dictates the adequate level of fidelity has important ramifications from a computational-cost perspective. Further numerical analysis with high-fidelity tools is typically performed to assess the validity of the low-fidelity results and support the decision on the adequate level of fidelity. This approach has been used in industry and research centers [45–47]. Multiple strategies to bridge the gap in computational cost and capabilities between high and low-fidelity potential flow solutions represented by the 3D Euler Equations and simple Vortex Lattice and Doublet Lattice Method (VLM and DLM) formulations, respectively have been presented in the literature. A practical approach is to employ medium-fidelity aerodynamics. This level of fidelity is represented by advanced 3D Panel Methods. These can allow for a more detailed geometrical modeling of aeronautical structures, including non-planar aerodynamic shapes. Furthermore, by employing advanced wake modeling techniques, reasonable approximations for re-circulation effects and induced drag can be achieved at a lower cost than high-fidelity solutions. Kennedy and Martins [48] have employed a 3D

Panel Method in a parallel aerostructural optimization framework and have additionally checked the validity of the obtained aerodynamic solutions against high-fidelity ones. Mieloszyk and Goetzendorf-Grabowski [49] have used a 3D Panel method in conjunction with corrections for viscous effects to evaluate the flight dynamics stability constraints in the multidisciplinary optimization of a joined wing box configuration. Recently, a finite element beam model was coupled to 3D Panel Aerodynamics and employed in aeroelastic optimization by Conlan-Smith and Schousboe Andreasen [50].

A natural path to reduce the computational cost associated with increased fidelity is to search for an optimized solution in smaller regions of the parameter space (thus requiring a smaller number of optimization iterations). This approach requires a priori some knowledge of the region in parameter space where the optimum is likely to be detected. In our multi-fidelity optimization framework, the low-fidelity part is entrusted with obtaining an initial optimized design with a high numerical certainty. Then, the high-fidelity modules explore further gains in performance by modeling the physics that the low-fidelity modules cannot. Overall, only a fraction of the iterations using highfidelity are performed (in the range of 5–10%). Such an approach prevents a drastic increase in computational cost. The procedure of gradually moving to increasing level of fidelity but with a smaller number of solutions is typical of a multi-fidelity optimization framework [51,52].

Still, a very high computational cost due to each high- fidelity solution is incurred, particularly so if a higher number of runs is needed. This cost is further exacerbated in the aeroelastic solutions. The cost of solving the unsteady aerodynamics with a high-fidelity method (typically a form of the Navier–Stokes equations with spatial approximation provided by the finite volume or element methods) tends to be higher than the structural solution. A combined computational cost of solving the structure and unsteady aerodynamics plus the interpolation of data between the two systems is then incurred. To mitigate this cost significantly, particularly for multiple runs around the same reference flow conditions, Order Reduction techniques can be used. Such techniques strive to approximate the physical behavior of the system on a greatly reduced (but physically meaningful) set of degrees of freedom. Most reduced order models are trained using samples from full solutions. The cost of obtaining those samples represents the main cost considering their use, since, afterward, they provide extremely rapid solution times. Several types of reduced order models (ROMs) have been used to treat aeroelasticity, including those based on (a) Mathematical decompositions, such as proper orthogonal decomposition (POD) [53,54], (b) Series approximations (for example, Volterra or Wiener) [55], and (c) Machine Learning techniques [56]. In our work, we integrate a ROM to the optimization framework that uses a Volterra series approximation for the unsteady aerodynamics coupled with mode superposition on the structural side. It serves as an intermediate step between low-fidelity aeroelastic solutions and the full coupled high-fidelity final run.

Although multiple analysis techniques of varying fidelity have been successfully used in optimization, the different computational cost that they incur combined with their setup complexity leads to a different combination of methods being adopted in each design stage of an aircraft structure. Specifically in industry, the current state-of-the-art is to include higher fidelity modules as design maturity increases (preliminary and detailed design stages), while exploiting the low computational cost of low-fidelity methods in the conceptual stage [47,57]. In the conceptual stage, efficient low-fidelity methods (panel method aerodynamics and equivalent beam or plate structural models) are commonly preferred due to the requirement to evaluate the feasibility of a vast number of design configurations.

Advances in computational capabilities combined with robust and optimized highfidelity analysis tools has resulted in a trend to incorporate such methods already from the preliminary stage. In order to bridge the gap between capabilities and computational cost between high and low-fidelity methods, two of the most attractive methodologies consist of employing correction factors to the lower-fidelity tools [37,58], as well as constructing

surrogate models. Furthermore, the construction of multi-fidelity optimization frameworks represents a highly active research topic. Such frameworks strive to combine the advantages of both worlds while mitigating the associated computational cost [47].

The primary research goal of the present study is to evaluate the gains in the structural optimization of a current state-of-the-art composite airliner wing by implementing a multi-fidelity optimization framework. Structural and aeroelastic constraints are used to formulate the optimization problem. The developed framework and the insight gained in the structural design and optimization of the composite wing model was used in the GRETEL project [59]. As a case study, a modified, planar (untwisted) version of the Common Research Model (CRM) [60] wing has been treated. The optimization framework is based on the Equivalent Plate Model (EPM) and the Mixed Integer Ant Colony Optimization (MIDACO) algorithm. At the first stage, an optimized configuration for the wing under low-fidelity aerodynamic loading is obtained. Higher fidelity aerodynamic solutions are then generated, and possible changes in the optimized solution are explored. Finally, the static aeroelastic response of the wing under consideration is also investigated. An overview of the performed actions of the developed optimization framework for the current case study is presented in Appendix B, highlighting input-output relations.

Although multiple studies in the literature assess the gains of using multi-fidelity optimization frameworks for innovative wing configurations, the gains for the currentstate-of-the-art wings are less prominently covered. The present study aims to contribute to closing this gap. We demonstrate that, through the use of this optimization tool, reasonable gains in structural mass of the test case wing were realized.

#### **2. Materials and Methods**

The developed multi-fidelity structural optimization framework is comprised of several distinct modules of varying fidelity and computational cost. The flowchart for the optimization framework is presented in Figure A4. These are combined in a manner that enables efficient turnaround times, as well as an optimized solution, with high levels of certainty accounting for the most relevant physics. In the next paragraphs, each of the consisting modules of the optimization framework will be mentioned, and their contribution, advantages, and disadvantages assessed. In addition, we will then present their combination in the developed optimization framework. The path from initial design to optimized solution referring to the potential changes at each separate stage will be described.

#### *2.1. Low-Fidelity Modules*

The present module, consisting mainly of low-fidelity computational tools and physics approximations, initializes the design process and, hence, is responsible for narrowing down the design space, while simultaneously steering the design towards optimum solutions, as well as accelerating the design procedure. The following sub-modules are included and described in detailed fashion:

#### 2.1.1. EPM Structural Model

Starting from the cornerstone of this work, the EPM enables a great reduction in the dimensionality of the problem since under the assumption of a relatively small thicknessto-chord ratio one can assume that an aircraft wing pertains a plate-like behavior. The problem is then reduced from a 3D representation to a 2D equivalent plate with smeared properties, reducing the associated DOFs. Plate kinematics, and specifically the First Order Shear Deformation Theory (FSDT) [61], can be utilized to describe and calculate the stress and strain state of each wing component. In particular, the wing skins, spar, and rib webs are treated as laminated plates, while, on the other hand, and as per common practice, spar and rib caps are assumed to resist axial loads only and are, thus, considered

as one-dimensional rods. To obtain the governing equation of the system the principle of minimum total potential energy is applied [18]. Under the FEA method, interpolation functions are used to describe the variation of the displacement field within an element, as well as to obtain the stiffness and mass matrices, [*K*] and [*M*], respectively. Another important assumption of the EPM is that the superposition principle holds; hence, the individual stiffness and mass matrices of each wing component can be summed in order to obtain the respective wing global matrices:

$$\begin{aligned} [\mathbf{K}]\_{\text{win}\mathfrak{g}} &= [\mathbf{K}]\_{\text{skin}} + [\mathbf{K}]\_{spar\_{\text{cap}}^{\text{wbb}}} + [\mathbf{K}]\_{rib\_{\text{cap}}^{\text{wbb}}} \\ [\mathbf{M}]\_{\text{win}\mathfrak{g}} &= [\mathbf{M}]\_{skin} + [\mathbf{M}]\_{spar\_{\text{cap}}^{\text{wbb}}} + [\mathbf{M}]\_{rib\_{\text{cap}}^{\text{wbb}}} \end{aligned} \tag{1}$$

Let *F* be a generic function that describes the stiffness and mass distribution along the threedimensional space for each wing component. The contribution of each EPM element to the global matrices is derived via the following integrals. Gaussian numerical integration is employed for the solution of the integrals, along with a coordinate transformation between the physical coordinates (*x*, *y*) and the generalized coordinates (*ξ*,*η*). The limits of integration in the through-the-thickness (z axis) coordinate correspond to the maximum extent of the local part of the actual wing structure.

$$\int\!\!\int\!\!\int\_{V} F(x,y,z)dV = \int\_{-1}^{1}\!\!\int\_{-1}^{1} \left(\sum\_{n=1}^{N\_{\tilde{x}}} \int\_{z\_{in}}^{z\_{j\tilde{n}}} \!\!\{x(\xi,\eta),y(\xi,\eta),z\} \ |f| \,dz\right)d\xi \,d\eta,\tag{2}$$

where *Nz* is the number of integration points in the z-direction, *zin* and *zjn* the integration limits of the n-th integration point, and |*J*| the determinant of the Jacobian of the coordinate transformation. Each wing component pertains different geometric characteristics; hence, specific treatment is required for the integration procedure:

1. Skins:

$$\int\!\!\int\!\!\int\_{V} \mathbf{F}(\mathbf{x}, y, z)dV = \int\_{-1}^{1} \int\_{-1}^{1} \left(\int\_{z\_{L} - \frac{1}{2}t\_{L}}^{z\_{L} + \frac{1}{2}t\_{L}} \mathbf{F} \cdot |f| \, dz + \int\_{z\_{ll} - \frac{1}{2}t\_{L}}^{z\_{L} + \frac{1}{2}t\_{L}} \mathbf{F} \cdot |f| \, dz\right) \, d\xi \, d\eta\_{\prime} \tag{3}$$

where *tL*,*<sup>U</sup>* = *t*<sup>0</sup> 1 + *tan*<sup>2</sup>*aL*,*U*, the thickness of the upper and lower skins, respectively, and *aL*,*<sup>U</sup>* the local airfoil angle of the lower and upper skin, as illustrated in Figure 1a. For a laminated plate, the numerical integration in the z-direction is performed for each lamina individually, along with its corresponding constitutive matrix. The rest of the components contribute in an analogous manner, albeit the integration limits are modified according to their relevant geometric characteristics, as presented in Figure 1b.

2. Spar Webs:

$$\begin{split} \int \left| \int\_{V} \mathbf{F}(\mathbf{x}, \boldsymbol{y}, \boldsymbol{z}) dV \right| &= \int\_{-1}^{1} \int\_{\frac{\pi}{6} (n) - \frac{t\_{1}}{c}}^{\frac{\pi}{6} (n) + \frac{t\_{1}}{c}} \int\_{\frac{\pi}{2} (\boldsymbol{z} - \frac{1}{2} t\_{L} + h\_{1})}^{\frac{\pi}{6} (n - \frac{1}{2} t\_{L} - h\_{1})} F \left\{ \mathbf{x}(\boldsymbol{\xi}, \boldsymbol{\eta}), \boldsymbol{y}(\boldsymbol{\xi}, \boldsymbol{\eta}), \boldsymbol{z} \right\} \mid \boldsymbol{j} \, |dz d\boldsymbol{\xi} d\boldsymbol{\eta} = \\ \int\_{-1}^{1} \int\_{-1}^{1} \frac{t\_{1}}{c} \int\_{\frac{\pi}{2} (\boldsymbol{z} - \frac{1}{2} t\_{L} + h\_{1})}^{\frac{\pi}{2} (n - \frac{1}{2} t\_{L} - h\_{1})} F \left\{ \mathbf{x} \left[ \frac{t\_{1}}{c} \boldsymbol{\xi} + \boldsymbol{\xi}\_{s}(\boldsymbol{n}), \boldsymbol{n} \right], \boldsymbol{y} \left[ \frac{t\_{1}}{c} \boldsymbol{\xi} + \boldsymbol{\xi}\_{s}(\boldsymbol{n}), \boldsymbol{n} \right], \boldsymbol{z} \right\} \mid \boldsymbol{j} \, |dz d\boldsymbol{\xi} d\boldsymbol{\eta} \end{split} (4)$$

3. Spar Caps:

$$\begin{split} \iiint\_{V} \mathbf{F}(\mathbf{x}, y, z)dV &= \int\_{-1}^{1} \int\_{\tilde{z}\_{s}(n) - \frac{h}{c}}^{\tilde{z}\_{s}(n) + \frac{h}{c}} \left( \int\_{z\_{L} + \frac{1}{2}t\_{L}}^{z\_{L} + \frac{1}{2}t\_{L} + h\_{1}} F\{x(\tilde{z}, \eta), y(\tilde{z}, \eta), z\} \mid |I| dz \\ &+ \int\_{z\_{L} - \frac{1}{2}t\_{L} - h\_{1}}^{z\_{L} - \frac{1}{2}t\_{L}} F\{x(\tilde{z}, \eta), y(\tilde{z}, \eta), z\} \mid |I| dz \right) d\tilde{z} d\eta = \\ &\int\_{-1}^{1} \int\_{-1}^{1} \frac{l\_{1}}{c} \left( \int\_{z\_{L} + \frac{1}{2}t\_{L}}^{z\_{L} + \frac{1}{2}t\_{L} + h\_{1}} F\left\{x\left[\frac{l\_{1}}{c}\xi + \xi\_{s}(n), n\right], y\left[\frac{l\_{1}}{c}\xi + \xi\_{s}(n), n\right], z\right\} \right) |I| dz \\ &+ \int\_{z\_{L} - \frac{1}{2}t\_{L} - h\_{1}}^{z\_{L} - \frac{1}{2}t\_{L} - h\_{1}} F\left\{x\left[\frac{l\_{1}}{c}\xi + \xi\_{s}(n), n\right], y\left[\frac{l\_{1}}{c}\xi + \xi\_{s}(n), n\right], z\right\} |I| dz \right) d\tilde{z} d\eta \end{split} (5)$$

4. Rib Webs:

$$\begin{split} \int \int\_{V} \int\_{V} \mathbf{F}(\mathbf{x}, y, z) \, dV &= \int\_{-1}^{1} \int\_{\eta\_{r}(\xi) - \frac{t\_{2}}{s}}^{\eta\_{r}(\xi) + \frac{t\_{2}}{s}} \int\_{z\_{L} + \frac{1}{2}t\_{L} + h\_{2}}^{z\_{Ll} - \frac{1}{2}t\_{L} - h\_{2}} \, F\{\mathbf{x}(\xi, \eta), y(\xi, \eta), z\} \, |f| dz d\xi d\eta = \\ \int\_{-1}^{1} \int\_{-1}^{1} \frac{t\_{2}}{s} \int\_{z\_{L} + \frac{1}{2}t\_{L} + h\_{2}}^{z\_{Ll} - \frac{1}{2}t\_{L} - h\_{2}} \, F\left\{\mathbf{x}\left[\overline{\mathbf{g}}, \frac{t\_{2}}{s}\eta + \eta\_{r}(\xi)\right], y\left[\overline{\mathbf{g}}, \frac{t\_{2}}{s}\eta + n\_{r}(\xi)\right], z\right\} |f| dz d\xi d\eta \end{split} \tag{6}$$

5. Rib Caps:

$$\begin{split} \left| \iiint\_{V} \mathbf{F}(\mathbf{x}, y, z) dV = \int\_{-1}^{1} \int\_{\eta \cdot (\xi) - \frac{l\_{2}}{s}}^{\eta\_{\tau} \cdot (\xi) + \frac{l\_{2}}{s}} \left( \int\_{z\_{L} + \frac{1}{2}t\_{L}}^{z\_{L} + \frac{1}{2}t\_{L} + h\_{2}} F\{\mathbf{x}(\xi, \eta), y(\xi, \eta), z\} \right) |f| dz \\ &+ \int\_{z\_{L} - \frac{1}{2}t\_{L}}^{z\_{L} - \frac{1}{2}t\_{L}} F\{\mathbf{x}(\xi, \eta), y(\xi, \eta), z\} |f| dz \bigg) d\xi d\eta = \\ \int\_{-1}^{1} \int\_{-1}^{1} \frac{l\_{2}}{s} \left( \int\_{z\_{L} + \frac{1}{2}t\_{L}}^{z\_{L} + \frac{1}{2}t\_{L}} F\left\{ \mathbf{x}\left[ \underline{\xi}, \frac{l\_{2}}{s}\eta + \eta\_{r}(\xi) \right], y\left[ \underline{\xi}, \frac{l\_{2}}{s}\eta + \eta\_{r}(\xi) \right], z \right\} |f| dz \\ &+ \int\_{z\_{L} - \frac{1}{2}t\_{L} - h\_{2}}^{z\_{L} - \frac{1}{2}t\_{L} - h\_{2}} F\left\{ \mathbf{x}\left[ \underline{\xi}, \frac{l\_{2}}{s}\eta + \eta\_{r}(\xi) \right], y\left[ \underline{\xi}, \frac{l\_{2}}{s}\eta + \eta\_{r}(\xi) \right], z \right\} |f| dz \right) d\xi d\eta. \end{split} (7)$$

where *h*1,*h*<sup>2</sup> the thickness of a spar/rib cap, *l*1,*l*<sup>2</sup> the width of a spar/rib cap, and *t*1,*t*<sup>2</sup> the thickness of a spar/rib web, respectively.

**Figure 1.** Wing cross-section.

2.1.2. Low-Fidelity Aerodynamics

One of the most attractive merits of the EPM is the facilitated coupling with contemporary aerodynamic panel methods, and specifically the Vortex Lattice Method (VLM) [62], since both methods are based on flat-plate representations of their respective discipline. The VLM, a well-established numerical method used in the early design stages of aircraft design, assumes an incompressible, inviscid and irrotational flow field, along with horseshoe vortex elements, in order to model an aerodynamic surface, thus neglecting the presence of thickness and viscosity in the solution. For an aerodynamic mesh consisting of N panels, at each collocation point i, the contribution of all the horseshoe vortices is added, thus generating an Aerodynamic Influence Coefficient (AIC) matrix. Introducing, also, the

Neumann boundary condition of zero normal velocity across the surface, we obtain the following set of linear equations, which are solved for the unknown vortices strength, **[Γ***j***]**:

$$[AIC\_{i\bar{j}}] \* [\Gamma\_{\bar{j}}] = [b\_{\bar{j}}],\tag{8}$$

where **[***bj***]** the product of the free stream velocity and the panel surface normal vectors.

Tornado VLM [63], a VLM implementation in MATLAB, has been used in our case for the generation of the critical aerodynamic loading case, corresponding to a 2.5G pull-up maneuver at sea-level conditions. Prior studies [64,65] have identified the 2.5G condition at a Mach number of 0.64 as the critical load case for the studied wing, as summarized in Table 1. Equation (8) can then be evaluated, and the aerodynamic pressure is distributed to the nodes of each element of the EPM model, since similar meshes for the structural and aerodynamic analysis have been considered. Nevertheless, more elaborate 2D interpolation schemes for load transfer purposes between the two meshes can be easily adapted, as demonstrated in subsequent modules.

**Table 1.** Critical aerodynamic loading summary as a function of the Maximum Take-Off Weight (MTOW).

#### 2.1.3. Optimization Framework

The present low-fidelity module is completed by integrating the EPM, along with the aerodynamic loading procedure, into an efficient optimization scheme for the sizing of the structural components of the wing. In a typical structural optimization scheme, the minimization of the wing structural mass is set as the objective function with stiffness, static strength, modal and dynamic aeroelastic constraints completing, and rationalizing the optimization problem. Concerning the optimization variables, the number of the 0◦, (45◦, −45◦), 90◦ plies of a baseline layup, corresponding to a component of the wing, and expressed as multiples of an integer number of the material under consideration, along with the depth of the spar and rib caps, constitute the design variables vector [66]. As an example, the composite layup parametrization technique for a wing component is presented in Equation (9).

$$|0\_{(n\_i)}/45\_{(n\_{i+1})}/-45\_{(n\_{i+1})}/90\_{(n\_{i+2})}|\_s \tag{9}$$

where *n* an integer-valued variable allowed to range between its upper and lower bounds, and *i* an internal variable counter, ranging from 1 to the number of variables of the optimization problem. This formulation was deemed advantageous due to a straightforward and less complex implementation than advanced lamination parameters techniques, such as References [3,4,11]. The results, expressed as the required integer number of plies to achieve an optimized percentage of different ply orientations up to the required thickness is deemed satisfactory for the purpose of this study.

Conformity to composite materials design guidelines [67], by means of generating symmetric and balanced lay-ups, is achieved by setting the number of 45◦ and −45◦ plies equal. The upper and lower bounds, along with the nature of the variables under consideration, are also presented in Table 2. The spar/rib caps depth variables are allowed to vary between their designated bounds with a 0.5 mm step. To allow for a wider design space and an increased structural design freedom, the wing is divided into 6 spanwise evenly spaced zones, as illustrated in Figure 2.

**Figure 2.** Optimization zones.

**Table 2.** Optimization variables type, lower and upper bounds.

Regarding the constraints of the optimization problem, and although not clearly stated as design requirements per regulation, static stiffness constraints are accounted for in the majority of optimization studies of aircraft wings present. For the present study, the maximum deflection at the tip of the wing is constrained, along with a maximum twist angle induced at the tip of the wing [68]. The corresponding values are chosen to represent typical maximum deformation values reached for current state-of-the-art airliner wings [69]. Additionally, aeroelastic analysis was carried out to ensure safety against divergence and flutter. To avoid possible local structural designs during the optimization process and to ensure the required amount of modes for an accurate dynamic aeroelastic analysis, the first eigenfrequency is also constrained. The static strength of the structure is also examined via the following algorithm:

To avoid local maximum stress driven designs, while simultaneously keeping the number of constraints for the optimization problem to a minimum, the constraint aggregation technique of Kreisselmeier–Steinhauser (KS) has been employed [70,71]. For an

optimization problem consisting of *Nc* constraints, *g*, having a maximum value of *gmax*, the KS functions are of the following form:

$$\text{KS}\left(\mathcal{g}\_{j}\right) = \mathcal{g}\_{\text{max}} + \frac{1}{\rho} \ln[\sum\_{j}^{N\_{\mathcal{E}}} \exp\left(\rho(\mathcal{g}\_{j} - \mathcal{g}\_{\text{max}})\right)] \le 1 \tag{10}$$

and are formed separately for the calculated FI of each of the components involved in the static strength evaluation procedure. The aggregation parameter *ρ* is set to 100 for all KS functions. This value has been reported in structural optimization studies of the CRM wing [72,73] and is shown to provide accuracy of the optimal solution. The aeroelastic stability by means of flutter speed is also investigated for each candidate design. A MATLAB code developed by NASA, EZ-ASE [74], is utilized for generating the unsteady aerodynamic loads calculation and conducting the flutter analysis in the frequency domain via the p-k method. EZ-ASE features flat-plate aerodynamics and particularly the VLM and DLM, thus allowing for straightforward coupling with the existing structural mesh of the EPM model; however, since, for higher reduced frequencies, a denser aerodynamic mesh is needed for accuracy purposes, a bilinear interpolation scheme is used within. The optimization objective function and constraints are summarized in the following, Table 3.

**Table 3.** Optimization problem setup.

The flowchart of the resulting low-fidelity optimization framework is presented in the following, Figure 3. The MIDACO solver [75], adopting a combination of an extended Ant Colony optimization algorithm (ACO) [76], along with the Oracle Penalty Method [77], an advanced method developed for metaheuristic search algorithms for constraint handling of the solution process has been chosen for carrying out the optimization problem. A sequential approach with multiple runs has been adopted, implying that the predefined number of runs is divided into multiple runs pertaining different algorithmic parameters. Particularly, initial runs are mainly focused on extensive design space exploration and are accompanied by a relaxed constraint satisfaction tolerance. As the solution advances, the search becomes increasingly local by tweaking accordingly the internal FOCUS parameter that forces the MIDACO solver to focus mostly on the current best solution. In particular, the ACO algorithm implemented in MIDACO generates samples of iterates based on multikernel Gaussian probability density functions (PDF). For a generic variable k with upper and lower bounds *xu* and *xl*, respectively, the FOCUS parameter applies an upper bound for the standard deviation of a Gaussian PDF given by *xu*(*k*)−*xl*(*k*) *FOCUS* and *max*( *xu*(*k*)−*xl*(*k*) *FOCUS* , <sup>√</sup> <sup>1</sup> *FOCUS* ) for continuous variables and integer variables, respectively. As a result, smaller values of the FOCUS parameter is recommended for the initial runs, with larger ones used for refinement purposes. In parallel, the constraint satisfaction tolerance is tightened. At each succeeding run, the previous best solution obtained serves as the starting point for the current run. As stated earlier, this procedure is repeated for a predefined number of runs, satisfying user-defined stopping criteria. The parameters of the optimization approach are summarized at the following Table 4.

**Table 4.** Optimization algorithm parameters.

**Figure 3.** Low-Fidelity Module Optimization Flowchart (Figure A4).

#### *2.2. High-Fidelity Modules*

The high-fidelity analysis modules encompass aerodynamics of the undeformed body and static and dynamic aeroelasticity. The aerodynamic component is treated using the OpenFOAM framework while the high-fidelity structural system is formulated in ANSYS APDL Mechanical. The aerodynamic solution is provided by the Navier–Stokes equations discretized using the Finite Volume Method (FVM), including turbulence modeling using the Spalart-Allmaras model. Either the EPM approach described above or a 3D-shell 3D-beam FEA model is used to represent the structure. Each of those components are discussed in detail below.

2.2.1. Steady (Undeformed Body) Aerodynamic Loads Module

The main components that comprise this module are:

For the CFD solution, the 3D Navier–Stokes equations in compressible form are discretized with the finite volumes method and are solved with an implicit integration scheme. Turbulence modeling in the form of RANS methods can be included. There is voluminous bibliography on CFD with Finite Volumes Methods, to which the reader is referred for a detailed description of the respective numerical methods [78,79]. Below, only

a brief overview of the CFD formulation for the studied wing geometry will be provided, with emphasis on important choices that influence the solution procedure.

The compressible Navier–Stokes equations in integral form were solved. The spatial terms of the Navier–Stokes equations are discretized on the computational finite volume mesh with appropriate discretization schemes. The mesh itself is a hybrid C-grid consisting mainly of unstructured tetrahedron volumes but using inflated prism layers near the wing's surface to achieve the first cell wall distances (Y+) necessary for capturing the boundary layer using RANS turbulence models. In particular, a Y+ of 50–70 was targeted. The CFD grid is shown in Figure 4.

**Figure 4.** Finite Volume C-Grid Mesh for the coarse density (Table A4).

The unstructured internal volume mesh is advantageous for mesh deformation purposes, as will be emphasized in the description the aeroelastic module at the next stage of analysis. The Pressure Implicit Splitting of Operators (PISO) method is used to treat the problem. More details can be found in Demirdzic et al. [80].

The spatial convection terms are discretized using second order upwind schemes. These schemes offer a good balance between accuracy and mitigation of numerical oscillations. The spatial diffusion and source terms are discretized using central differences.

First order implicit time integration is used for the numerical solution. The form of the system of equations is represented in Equation (11):

$$
\boldsymbol{\phi}^{n} = \boldsymbol{\phi}^{n-1} + \boldsymbol{\Delta t} \cdot \boldsymbol{f}(\boldsymbol{\phi}^{n}),
\tag{11}
$$

where *φ* denotes any flow quantity that takes part in the Navier–Stokes equations.

The implicit formulation allows to use a wide range of timesteps with Courant numbers higher than unity. This is important for using the same mesh in the aeroelastic system as the flexibility in selecting the timestep makes the synchronization of data transfers between the structural and aerodynamic systems easier.

Turbulence was successfully treated with the Spalart-Allmaras model with wall functions providing near wall treatment. This turbulence model is integrated in most modern

CFD solvers. A first cell Y+ of 30–100 is required for this formulation. Due to the size of the wing, turbulence models that require lower values of Y+ incur a very high computational cost. The mesh convergence study at a load corresponding to a 2.5G maneuver is provided in Table A4.

For the structural component, the finite element method using either the EPM (lowfidelity) or a 3D-Shell (Shell 181 of ANSYS MEchanical APDL), 3D-Beam (Beam 188 of ANSYS MEchanical APDL) medium fidelity representation of the wing structure is used. The latter uses 4-node, 3D layered Shell Elements to discretize the wing skins, spar and rib webs, and 3D beam elements to discretize the spar and rib caps. Linear shape functions were used. The structural mesh is presented below in Figure 5.

**Figure 5.** Finite Element Structural Mesh.

In the current module, this component was only used to receive the steady aerodynamic loads and provide a steady state solution of the stress/strain field for optimization purposes. The interpolation of the aerodynamic loads from the CFD mesh to the structural mesh was realized using the GGI methodology [81]. This technique has seen widespread use in the transfer of loads between sliding mesh interfaces in CFD. A significant advantage compared to distance-based interpolation after projection of nodes upon a common boundary is that it takes into account the overlap between the elements. As a result, the area loads, such as pressure, are transferred consistently. This is highlighted in Figure 6 and Equation (12).

**Figure 6.** Transfer of pressure loads upon a common boundary with GGI.

$$P\_{\text{target}\,facc} = \sum\_{i=1}^{n} (\frac{A\_{\text{target}\,facc} \cap A\_{\text{source}\,facc}, i}{A\_{\text{target}\,facc}} P\_{\text{i}}). \tag{12}$$

A second interpolation layer is used if the loads must be transferred to the EPM. Then, after transfering them to the 3D Finite Element Mesh with the GGI approach, each node is mapped to the nearest ones of the EPM using a distance-based approach. This is less accurate, but, in any case, if maximum accuracy in interpolation is required the 3D-Shell, 3D-Beam structural formulation is preferred, and this layer not used. The complete flowchart of the High-Fidelity Steady Aerodynamic Loads Module is given in Figure 7.

**Figure 7.** Steady Airloads High-Fidelity Module (Figure A4).

2.2.2. Formulation of the Aeroelastic Problem and Order Reduction

The methodology to consistently interpolate surface pressure from the CFD to the structural mesh with the GGI technique was outlined in the previous subsection. To enable a coupled aeroelastic solution, the structural displacement must additionally be passed into the wing surface boundary of the CFD mesh and the internal mesh deformed, accordingly. For passing the structural displacements to the aerodynamic mesh, distancebased interpolation is sufficient. This contrasts with surface loads (aerodynamic pressure), for which the overlap of element face areas between the structural and aerodynamic mesh is taken into account by the GGI technique. An overview of the final staggered coupling scheme is given in Figure 8.

**Figure 8.** Staggered coupling scheme.

During deformation of the CFD mesh, the deterioration of element quality must be avoided. To that end, a Laplacian smoothing function with variable diffusivity based on distance from moving boundary is used. A detailed overview of the implementation of this mesh deformation technique, along with application examples, can be found in References [82,83]. Remeshing is only used if several cells violate the quality limits. The usage of unstructured tetrahedron internal meshes is advantageous because it is easier to maintain mesh quality upon deformation compared to structured meshes. A triangularshaped element can retain quality upon movement of each nodes much more robustly than a rectangular one. The same applies to 3D for tetrahedrons.

Order reduction was achieved by coupling a series approximation of the unsteady aerodynamics with a modal space description of the structure. The unsteady aerodynamics component of the aeroelastic reduced order model (ROM) is based on the Volterra series theory. Volterra series ROMs for aerodynamics have been pioneered by Silva [84] and have seen implementation in complex aeroelastic problems [85].

In our work, a more basic form of a Volterra series ROM is used for unsteady aerodynamics. Linear kernels up to 50–100 timesteps in length were used to model the unsteady aerodynamic response. Although the first 2–3 non-linear kernels of the series were built, their influence was minor. This is a result of the combination of flow regime (0.64 Mach at atmospheric conditions) and the behavior of test case wing that does not generally exhibit non-linear aerodynamic phenomena at those conditions.

A Volterra series consists of linear and/or non-linear convolution integrals. The final assembled series represents an expansion of the linear convolution theorem and associates the time history of input states to the one of the outputs. For discrete systems, the convolutions are formulated as multiplications of scaled and shifted characteristic responses (or kernels of the series) with the time history of the input. For a multi-input multi-output (MIMO) system, an output state is then dependent on the input states of the system and the kernels of the series as briefly presented in Equation (13):

$$\begin{aligned} y\_{out,i} &= y\_{0,i} + \sum\_{j=1}^{mode\,number} \left[ \sum\_{k=1}^{n} h\_{ji}(n-k) \cdot x\_j(k) \right] + \\ &\sum\_{j=1}^{mode\,number} \left[ \sum\_{k=1}^{n} \sum\_{l=1}^{n} H\_{ji}(n-k, n-l) \cdot x\_j(k) \cdot x\_j(l) \right], \end{aligned} \tag{13}$$

where:

*y*0,*<sup>i</sup>* denotes the initial state of the *i*th modal coordinate, *yout*,*<sup>i</sup>* is the output state of the *i*th modal coordinate, *xj* is the *j*th modal force,

*h*, *H* are the linear and non-linear kernels, respectively, and

*n* denotes the length of kernels in timesteps.

The Volterra series kernels were obtained by interpolating impulsive motions in each generalized coordinate to the aerodynamic system and obtaining the associated aerodynamic load response. For more information considering the construction of Volterra series models for aeroelasticity, the reader is referred to the corresponding references given in the start of this subsection.

Order reduction of the structural system was achieved by casting it in modal space and then formulating the resulting equations of motion in state-space form. Their form is provided in (14) and (15):

$$\begin{array}{l}\dot{[\dot{x}]} = A\_{ss}[\dot{x}] + B\_{ss}[F] \\\dot{[y]} = \mathbb{C}\_{\mathfrak{sl}}[\dot{x}] + D\_{\mathfrak{sl}}[F] \end{array} \tag{14}$$

$$A\_{ss} = \begin{bmatrix} [0]\_{n\times n} & [I]\_{n\times n} \\ [-\omega\_n^2]\_{n\times n} & [-2\zeta\_n\omega\_n]\_{n\times n} \end{bmatrix} \quad B\_{ss} = \begin{bmatrix} [0]\_{n\times n} \\ [I]\_{n\times n} \end{bmatrix} \tag{15}$$

$$\mathbf{C}\_{ss} = \begin{bmatrix} [I]\_{2n\times 2n} \end{bmatrix} \qquad \qquad D\_{ss} = \begin{bmatrix} [0]\_{2n\times n} \end{bmatrix} . \tag{15}$$

where *ωi*,*ζ<sup>i</sup>* are the natural frequencies and modal damping ratios.

The assembled aeroelastic ROM, together with the associated input parameters that can be used, is presented in Figure 9 below.

The output of the ROM consists of the time history of the state variables (modal forces and modal displacements). The aeroelastic oscillation frequency at each mode is obtained by pinpointing the position of maximum for the absolute value of the Fast Fourier Transform (FFT). A window of greater length than the period of oscillation is used to scan the signal for each mode. The window moves a period of oscillation corresponding to the aeroelastic oscillation frequency at each time if it is found to be non-zero; otherwise, this is a period of oscillation corresponding to the specific structural eigenfrequency. The

maximum of the absolute value of the windows is recorded. Termination is triggered when one of the following is true:

**Figure 9.** Aeroelastic reduced order model.

The third criterion is only used for dynamic aeroelastic analyses.

The time history of the state variables is also printed to file to allow for manual inspection. The complete module for aeroelasticity using order reduction as implemented in the interpolation framework is shown in Figure 10.

**Figure 10.** Implementation of aeroelasticity using order reduction to the optimization framework (Figure A4).

The final validation consists of a full order staggered (or loosely) coupled solution between high-fidelity aerodynamics and structural solvers. This full order module is run only once to validate the final optimized design candidate.

#### **3. Results and Discussion**

#### *3.1. Benchmark Solution*

A benchmark structural optimization case study is initially solved in order to enhance our confidence in the MIDACO solver. This problem refers to the optimization of the dimensions of the cross sections of a stepped cantilever beam under the application of a point load at the free end [86]. The minimization of the beam volume subject to various engineering design constraints formulate the mathematical optimization problem. In particular, the bending stress in each part of the cantilever should not exceed the maximum allowable stress, *σmax*. The deflection at the free end is also constrained by the maximum deflection, *δmax*. For each cross section, the height of the beam should also be twenty times less the corresponding width. The stepped cantilever beam, along with the constraints and the input data, are illustrated in the following, Figure 11.

**Figure 11.** Stepped cantilever beam benchmark problem.

Regarding the nature of the variables of the problem, the height and width (*Hi* and *Bi* for the *i*-th cross-section) of the first section are integers; for the second and third section, these are chosen from a discrete set, while, for the fourth and fifth, they are continuous, resulting into a mixed-integer design problem. Several optimization techniques, including the non-linear branch and bound method, simulated annealing or genetic algorithms, and approximation methods, were investigated in Thanedar and Vanderplaats [86]. The results of the optimization study, along with those obtained by the MIDACO solver, are reported in Table 5.

Overall, the optimum height and width of each cross section, as well as the overall beam volume values, are similar to the ones obtained by the MIDACO solver. The convergence history of the objective function and of the constraints satisfaction, expressed as the L-1 norm, are illustrated in the following, Figures 12 and 13.

**Table 5.** Benchmark problem optimization results.

**Figure 13.** Stepped cantilever beam benchmark problem.

A study on the effect of tightening the constraints satisfaction tolerance value is included in Appendix C.

#### *3.2. Main Case Study*

A modified version of the original CRM wing serves as a baseline model for the analyses described earlier, with the relevant geometric parameters being presented at Table 6. The baseline CRM wing is modified by neglecting the twist distribution, the resulting geometry being planar. This modification was necessary due to the limitations of the VLM code at our disposal. Additional effort would be required for the implementation of twist angle in the low-fidelity aerodynamics software, and, since the present work constitutes a structural optimization study, it was decided to proceed with a planar representation of the original CRM wing. Furthermore, the lack of accounting for the twist distribution constitutes a limitation of the current realization of the EPM.

**Table 6.** CRM wing geometric data.

As far as the internal configuration of the CRM wing is concerned, a three-spar configuration, with each spar laying at 10%, 30%, and 70% of the local chord, respectively, along with thirty-eight evenly spaced and aligned with the airflow ribs, similar to the one presented in Jutte et al. [87], is considered. Spar caps, along with full-depth rib caps of rectangular cross section, are also present. A view of the Outer Mold Line (OML), as well as of the internal structure of the CRM wing, can be seen in Figure 14. Concerning the material model and pertaining to all numerical models, all relevant components are assumed to be manufactured of composite materials with linear elastic behavior, and specifically the Hexcel 8552/IM7 prepreg [88] with unidirectional carbon fibers, whose properties are summarized in Table 7.

**Figure 14.** CRM wing OML and internal structure.

**Table 7.** Hexcel 8552/IM7 UD CFRP properties.

#### *3.3. Low-Fidelity Module*

The dimensionality reduction capability of the EPM in comparison with the corresponding 3D FEA model is initially investigated. As demonstrated in the following, Table 8, the EPM achieves a significant reduction of DOFs when compared with the 3D FEA model of the CRM wing, which eventually translates into faster computational time by a factor of 4.

**Table 8.** Numerical models size comparison.

The MIDACO optimizer was executed for the optimization problem presented in Table 3 with the parameters demonstrated in Table 4. In brief, the KS stress constraints represent the aggregate value of the Failure Indices of each wing component considered, with the limit value of 1 indicating failure under the maximum stress criterion. A strategy of tightening the constraints tolerance value was employed with the final value at each stage serving as an improved initial prediction for the next, as presented in Table 4. The convergence of the objective function, as well as of the constraint satisfaction, expressed via the L-1 norm, being demonstrated in the following, Figures 15 and 16.

**Figure 15.** Wing mass convergence.

**Figure 16.** Constraint satisfaction, expressed via the L-1 norm.

Convergence to a minimum mass, along with non-violating design constraints, has been achieved. The corresponding mass, as well as design constraints, at the optimized solution are presented in Table 9. The values of the constraints, and especially the static stiffness and strength related ones, approach their limit values, indicating a well approximated true optimum solution.

**Table 9.** Low-Fidelity module optimization results summary.

As far as the thickness distribution for the relevant wing components is concerned, a general spanwise decrease in thickness for the relevant components is demonstrated in Figures 17 and 18. Of particular interest is the fact that the aforementioned decrease appears to start at the outboard section of the wing, indicating a highly stressed region near the yehudi break, namely the coalescing region between the inboard and outboard sections of the wing. On the front of the individual thickness of each component, the thickness values associated with the wing skins are higher in comparison to the other components mainly due to the fact that they contribute more to the secondary moment of inertia of the wing. The thickness values of the spars are close to the corresponding skin thicknesses, albeit lower, being closer to the neutral axis of the wing section. Nevertheless, they take torsional loads; thus, relatively high thickness values are necessary. On the other hand, the wing ribs are associated with relatively low thickness values, as expected, since they are not the main load carrying components of an aircraft wing; they, however, aid at the relief of a portion of the developing shear stresses. The spanwise ply percentages for each component of the wing are illustrated in Figure 19.

**Figure 19.** CRM wing components spanwise layups percentages.

#### *3.4. High-Fidelity Module*

High-fidelity aerodynamics using the RANS formulation described in Section 2.2.1 were incorporated in structural optimization. The optimization was performed for a specified lift load (2.5G maneuver condition) rather than a constant angle of attack.

A single optimization pass that includes multiple high-fidelity solutions to pinpoint the new optimized solution is conducted. The main goal is to pinpoint a good balance between loss in lift and reduced structural mass. The change to the required load for a 2.5G maneuver due to the change in wing structural mass was accounted for. For the termination of the high-fidelity aeroelastic optimization, an attained lift within 5% of the rigid wing was demanded. A tighter termination criterion for the attainable lift could have been implemented; nevertheless, due to limited computational resources, it was deemed satisfactory for the purpose of this study.

Significant changes in the required angle of attack to attain the targeted lift load were obtained by the inclusion of high-fidelity aerodynamics and are summarized in Table 10.

**Table 10.** Change in required angle of attack to attain a 2.5G Lift Load.

The initial optimization pass uses static structural solutions for the wing with the loads provided by high-fidelity rigid body aerodynamics. Specifically, the loads were interpolated onto the EPM mesh, providing a more realistic loadset than VLM for the static solution. Following convergence to a new optimized design, a final optimization pass that treated the coupled problem by using a high-fidelity aeroelastic formulation of the wing was performed. In this stage, the loads are transferred directly onto the wetted surface of the 3D-wing model. Both stages fail to satisfy the minimum lift that is required for the 2.5G maneuver scenario at the angle of attack of the base geometry. The inclusion of viscous effects reduces the maximum lift compared to the VLM method, the main culprit being pressure recirculation effects near the wingtip and trailing edge.

Using high-fidelity rigid body aerodynamics, an increase of the angle of attack to 3.3 degrees was deemed necessary to attain the same *CL*. The static deformation of the new optimized wing geometry after the inclusion of the change in angle of attack is presented in Figure 20, for both the EPM method and a full 3D model of the structure.

**Figure 20.** Deflection to static aerodynamic RANS load. EPM (**left**); 3D Shell Model (**right**).

The reasonable agreement between the two approaches points to the success of the EPM structural representation and a proper analogy between the loads on the 3D model and the ones on the plate. The structural mass remained largely unchanged after the inclusion of the new loads. Possible mass gains were deemed too small (in the region of 1–2%) to conduct further optimization iterations at this stage. From an aerodynamic perspective, however, the change in angle of attack is not negligible.

The high-fidelity aeroelastic ROM pinpointed both a further loss in lift compared to the high-fidelity static structural solution and a lower aeroelastic tip deflection. The wing candidate is characterized by strong bending-torsion coupling at structural level, and, as a result, the aeroelastic twist tends to cause negative twist angle near the wingtip upon aeroelastic deformation. The loss in lift is mainly a result of this bend-torsion coupling of the wing combined with the planar (untwisted) geometry at the undeformed state; upon aeroelastic deformation, negative twist angle is attained near the wingtip resulting in loss of lift. Near the root the aeroelastic twist angle changes are smaller; thus, a portion of the lift continues to be generated. In order to compensate for the loss of lift, a higher angle of attack for the entire wing is selected. An additional factor that leads to a reduction in lift (although minor in this specific case study due to the requirement for a reasonable bending deflection) is the curvature of the outer part of the wing due to aeroelastic bending deflection. The generated load at each section of the wing, being normal to its surface, becomes slightly tilted with respect to the vertical axis upon aeroelastic deflection. A loss of a part of the total lift ensues due to the lift load becoming tilted with respect to the vertical axis. This reduction due to bending is minor compared to the one due to change in twist angle for the current case. However, it could be significant for wings that exhibit a high magnitude of tip deflection (for example, in high aspect ratio configurations). Consequently, a further increase in AoA compared to the rigid case is required to attain the target lift load, as presented in Table 10.

Additionally, the difference between negative twist angles at the tip and much lower aeroelastic twist values towards the root also provides an opportunity for a reduction in structural mass. At the new increased angle of attack, the negative twist at the tip is increased due to the higher aeroelastic deflection magnitude and a restoring bend-torsion coupling tendency. At the same time, a much smaller change in aeroelastic twist results near the root. Thus, an even greater portion of the of the lift is generated near the root, and the bending moment is reduced significantly as a result. This reduction in bending moment at a similar lift load causes a smaller maximum deflection for the coupled problem, as presented in Figure 21, in which the static (right) and static aeroelastic (left) equilibrium deformation for the final design candidate of each stage are compared. The reduction of bending moment due to aeroelastic twist cascades to a reduction in structural mass as captured in Table 11.

**Table 11.** Change in structural mass.

\* high-fidelity aerodynamics.

In a non-planar baseline backswept wing with an optimized twist distribution, lower gains are expected. Since negative twist angles are already used near the wingtip together with a higher incidence angle at the root a lower bending moment even for the static case would be obtained compared to the planar wing results. It needs to be, therefore, stressed that the non-optimized aerodynamic shape of the considered wing (due to lack of proper twist distribution) exacerbates the reductions in structural mass that were attained in the present study.

A structural mass reduction in the order of 9% was realized using a static aeroelastic description of the wing in optimization as opposed to the static high-fidelity one. A loss of approximately 4% resulted. The final mass is compared against the static case in Table 11.

The wing was additionally checked against divergence and flutter using high-fidelity aeroelasticity. Divergence was pinpointed as the dominant aeroelastic instability, but the dynamic pressure of occurrence is more than three times higher than the operational one for a constant Mach number of 0.64. The modified wing is, therefore, Flutter- and Divergence-free at the considered point of the flight envelope, confirming the results of the low-fidelity aeroelastic analysis.

**Figure 21.** Comparison of static deflection distribution versus aeroelastic solution (*z*-axis scale not actual). Note that a higher angle of attack was required for the aeroelastic solution to attain the same coefficient of lift. Only deformation is shown in the figure; angle of attack is not plotted.

#### **4. Conclusions**

A multi-fidelity optimization framework for current state-of-the-art composite aircraft wings, based on the Mixed Integer Distributed Ant Colony Optimization (MIDACO) has been presented. The main case study consisted of a modified CRM wing with a planar geometry. The structure was represented using the Equivalent Plate Method (EPM) or a 3D shell and rod model, while aerodynamics were treated with a Vortex Lattice Method (VLM) or Reynolds Averaged Navier–Stokes (RANS) solutions, depending on level of fidelity. The aeroelastic solutions were provided by Reduced Order Models (ROMs). Upon completion of the optimization procedure, optimized layups for the various wing components were obtained, as illustrated in Figures 18 and 19. As a general trend, a spanwise decrease in thickness was observed, as expected. A region near the yehudi break was accompanied with an increase in the relevant thicknesses, indicating that this region is in fact a highly stressed one due to the slope change in geometry of the wing.

Comparison between the different fidelity modules was performed with the aim of highlighting the merits and contribution of each to the optimization process. The EPM enabled a great reduction of the associated DOFs of the problem, as discussed in Table 8 and in Section 3.3, while a good comparison with the high-fidelity structural model was demonstrated in Figure 20. Significant changes to the angle of attack to attain the target lift load resulted from the inclusion of high-fidelity aerodynamics and aeroelasticity. The differences have been attributed to negative aeroelastic twist angles near the wingtip and resulted in a lift degradation in the order of 4% and reduced bending moment which cascaded to a structural mass reduction in the order of 9%. It should be stressed that the reduction in mass is exacerbated by the planar baseline geometry of the wing since, for an optimized twist distribution, a smaller bending moment is exhibited, even in the undeformed state, as discussed in Section 3.4. The inability of the current framework to account for a realistic twist angle distribution at all levels of fidelity is acknowledged as an important shortcoming.

Overall, the present methodology presents a successful blend of current state-of-the-art multi-fidelity structural, aerodynamic, and fluid-structure interaction analysis for structural optimization. It is shown that employing such frameworks can result in performance gains even for current state-of-the-art-wings. Further possible enhancements to the current optimization framework could include the modification of the VLM and EPM code to account for twist angle distributions, implementation of buckling constraints in the highfidelity structural model, and incorporation of aeroelastic tailoring capabilities. Another major addition would be the development aerodynamic shape optimization procedures.

**Author Contributions:** Conceptualization, A.K. (Angelos Kafkas) and S.K.; methodology, A.K. (Angelos Kafkas) and S.K.; software, A.K. (Angelos Kafkas) and S.K.; validation, A.K. (Angelos Kafkas) and S.K.; formal analysis, S.K., A.K. (Angelos Kafkas), A.K. (Athanasios Kotzakolios), V.K. and G.L.; investigation, S.K. and A.K. (Angelos Kafkas); resources, S.K., A.K. (Angelos Kafkas), A.K. (Athanasios Kotzakolios), G.L. and V.K.; data curation, S.K. and A.K. (Angelos Kafkas); writing original draft preparation, S.K., A.K. (Angelos Kafkas), A.K. (Athanasios Kotzakolios), V.K. and G.L.; writing—review and editing, S.K., A.K. (Angelos Kafkas), A.K. (Athanasios Kotzakolios), G.L. and V.K.; visualization, S.K., A.K. (Angelos Kafkas), A.K. (Athanasios Kotzakolios), V.K. and G.L.; supervision, A.K. (Athanasios Kotzakolios), V.K. and G.L.; project administration, A.K. (Athanasios Kotzakolios), V.K. and G.L.; funding acquisition, S.K., A.K. (Athanasios Kotzakolios), V.K. and G.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was financially supported by the State Scholarships Foundation of Greece (IKY), which granted a scholarship to Spyridon Kilimtzidis. The scholarship is co-financed by Greece and the European Union (European Social Fund—ESF) through the Operational Programme << Human Resources Development, Education and Lifelong Learning >> in the context of the project "Strengthening Human Resources Research Potential via Doctorate Research" (MIS-5000432), implemented by the State Scholarships Foundation (IKY).

The authors also gratefully acknowledge the financial support of EU, H2020 CS2 project under the acronym GRETEL (Grant agreement ID: 737671).

**Institutional Review Board Statement:** Not applicable. **Informed Consent Statement:** Not applicable. **Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

#### **Appendix A. EPM Validation Test Case**

In this section, the accuracy and efficiency of the EPM is thoroughly investigated via a test case concerning a built-up wing composed of composite materials skins, spar/ribs webs and caps. A typical internal geometry, consisting of two spars placed at 20% and 70% of each chord and 12 evenly spaced and aligned with the airflow ribs is considered. In addition, spar and rib caps of rectangular cross-section are also present. The relative dimensions of the various components were selected such that no local modes, which the EPM is incapable of capturing, will be present to the results as much as possible. Additionally, two numerical models are developed and compared, namely a 3D FEA model in NASTRAN and its corresponding EPM model. Regarding the latter, the FEA is employed for the solution of the numerical integrals present in the Materials and Methods section, with the solution procedure being carried out entirely in MATLAB. Initially, and based on the external geometry of the corresponding wing, the geometry of the equivalent plate is generated, upon which a FEA mesh is built. To avoid shear locking phenomena and reduce the spurious mechanisms present, isoparametric nine-noded Lagrange plate elements with selective integration for the membrane, bending and shear terms have been chosen for the analysis. Eigenfrequencies and eigenmodes extracted from a free vibration analysis constitute the quantities of interest for the comparison between the two methods, since both the stiffness, as well as the mass matrix, are involved in this type of analysis. Concerning the 3D FEA mesh, the wing upper and lower skins, along with the spar and rib webs, are modeled via shell elements, (CQUAD4), while the spar and rib caps are assumed to resist axial loads only and are, thus, modeled via rod elements (CROD). For the EPM, QUAD-9 elements have been chosen to model the equivalent plate. The particulars of the two FEA meshes along are presented in Table A1, while the resulting FEA mesh for both of the methods is presented in the following Figure A1a,b. Each component pertains a [45/0/ −45/90]s quasi-isotropic, balanced and symmetric layup, with similar material properties to the one presented in the Results section. Regarding the boundary conditions, the wing is considered clamped at its root for both of the numerical models, fixing all 6 DOFs.

**Table A1.** Numerical models comparison.

The efficiency of the EPM is further emphasized for this test case, since the addition of more components (spars/ribs webs and caps) is accompanied by an increase in the DOF's for the 3D FEA, which is not the case for the EPM, where the associated DOFs remain unaltered. The natural frequencies, mode shapes (as illustrated in Figures A2 and A3) and the mass of the wing for the two methods, presented in Table A2, are also examined.

**Figure A1.** FEA & EPM mesh.

**Table A2.** NASTRAN-EPM natural frequencies and mass.

**Figure A2.** Eigenmodes comparison—Modes 1–3.

**Figure A3.** Eigenmodes comparison—Modes 4–6.

A close agreement between the shapes of the eigenmodes is observed in general, with the EPM being capable of predicting the trends of the 3D-FEA model.

**Figure A4.** Flowchart of the optimization framework.

**Figure A5.** Action list of the developed optimization framework for the CRM case.

#### **Appendix C. Benchmark Solution Tolerance Value Study**

**Table A3.** Effect of tolerance on the optimized solution of the benchmark case study.