Multidisciplinary structural optimization of novel high-aspect ratio composite aircraft wings

Novel high-aspect ratio airframe designs pave the way for a more sustainable aviation future. Such configurations enhance the aerodynamic efficiency of an aircraft through induced drag reduction mechanisms. Further performance gains, mainly in terms of structural mass, are accomplished via composite materials airframes. Nevertheless, undesired phenomena such as geometric nonlinearities and aeroelastic couplings due to elevated flexibility may often rise, rendering the design and optimization of such airframes extremely intricate and prohibitive in terms of computational cost. Low-fidelity tools, often preferred on the early design stages, accelerate the design process, albeit suffering from reduced accuracy and ability to capture higher-order phenomena. Contrastingly, high-fidelity computational methods incur excessive computational cost and are therefore utilized at the later, detailed design stages. There arises, therefore, the need for a combination of the various fidelities involved in a cost-effective manner, in order to drive the design towards optimal configurations without significant performance losses. In our approach, variable fidelity analyses are initially conducted in order to shed light on their effect on the structural response of a high-aspect ratio composite materials reference wing. An optimization framework combining low and high-fidelity tools in a sequential manner is then proposed, aiming at attaining a minimum mass configuration subject to multidisciplinary design constraints. As demonstrated, reasonable mass reduction was obtained for a future aircraft wing configuration.


Background
Ever since its inception at the dawn of the twentieth century, the airline industry has been dramatically expanding, with the latest predictions (Boeing 2018; Airbus 2018) foreseeing a 4.3% annual traffic growth over the next 20 years. This forthcoming rise in fuel consumption raises environmental concerns, mainly due to the associated increase of CO 2 and NO x concentration in the atmosphere that constitutes a severe public threat. On that end, the European Commission envisages through the definition of ACARE 2050 (2011), among others, a 75% reduction in CO 2 emissions for future commercial transport aircrafts. A key enabler for this effort is the introduction of novel, more efficient airframe designs with improved aerodynamic as well as structural efficiency. In the field of aerodynamics, the aerodynamic efficiency of an aircraft is closely related to the aspect ratio of the wings, with higher aspect ratios reducing the induced drag and increasing the lift-to-drag ratio, which in turn improves fuel efficiency. Therefore, and as clearly depicted in Fig. 1, an increasing trend in the aspect ratio of commercial aircraft wings over the past decades is observed. On the other hand, the advent of advanced composite materials in the aeronautics scene, replacing conventional metallic alloys, has enabled the development of lighter yet stiffer airframe configurations with aeroelastic tailoring and active load alleviation capabilities. Nevertheless, the aforementioned aircraft design novelties are accompanied by certain shortcomings, the most prominent being the elevated flexibility of the wing that on one hand induces a material and/or geometric nonlinear behavior and on the other hand a closer coupling between the structure and the surrounding fluid, aggravating static and dynamic aeroelastic phenomena. Furthermore, the use of composite materials perplexes the design optimisation procedure by vastly increasing the design space via the introduction of new variables such as ply thickness and ply orientation. As a result, the development of novel design methodologies and computational tools capable of predicting the whole spectrum of phenomena inherent to this new generation of structures as well as producing optimum configurations is of paramount importance.
Unfortunately, contemporary high-fidelity tools are commonly associated with elevated and often prohibitive computational time, hence their use is limited to the later, more detailed, design stages. Since in these stages the design freedom is bounded and any design alteration could induce severe cost penalties, as illustrated in Fig. 2, low and medium fidelity tools, capable of representing a portion of the phenomena that might rise, need to be exploited for the conceptual and preliminary design stages. In this manner, design space exploration is reinforced by alleviating the computational burden, steering the design process towards optimality.
Early efforts, confined by modest computational resources, focused primarily on low-fidelity representations of the structure as well as of the aerodynamic disciplines. Pioneering research was conducted by Triplett (1980) and Love and Bohlman (1989), where TSO, one of the earliest multidisciplinary design and optimization tools, was developed. An equivalent plate model of the structure coupled with the doublet lattice method (DLM) for aerodynamics were combined to optimize the thickness distribution and laminate orientations of a fighter wing subject to strength and flutter velocity constraints. Haftka (1973Haftka ( , 1977 was also among the first to combine computational tools of similar fidelity along with optimization algorithms to perform aeroelastic analysis and optimize aircraft wings under stress and drag constraints. The superiority of a composite wing design against the up-to-date traditional aluminum one was demonstrated in terms of mass and induced drag. The discontinuity of the flutter constraint was also pinpointed. In a similar fashion, Grossman et al. (1988) conducted an aerostructural optimization study of a sailplane wing using lifting-line aerodynamics along with beam equations and went on to optimize a transport aircraft wing, while developing methodologies for calculating sensitivity derivatives as well as a sequential approximate optimization module (Grossman et al. 1990).
With the advent of the computational prowess, higher fidelity computational tools and mostly computational fluid dynamics (CFD) methods started to emerge. Traditional panel methods were gradually getting replaced by nonlinear aerodynamic solvers and on the other hand beam or plate structural representation techniques were substituted by global finite element method (FEM) models. Aeroelastic effects were also introduced in the formal coupled analyses, due to their criticallity in the aerodynamic and structural design. Maute et al. (2001) presented a methodology for optimizing the aerodynamic and structural parameters of a steady-state aeroelastic system, using 3D Euler finite volume method (FVM) and FEM models. A plethora of novel solution methodologies was developed and applied, aiming at the single and multi-objective optimization of various wing test cases. A highly efficient and robust Schur-Newton-Krylov solution method to the coupled, steady-state aeroelastic problem was given in Barcelos et al. (2006), with the methodology being extended to also include Reynolds averaged Navier-Stokes (RANS) CFD (Barcelos and Maute 2008). Moreover, coupled aerostructural optimization studies started to draw the attention of the scientific community, with Reuther et al. (1999) presenting an integrated, highfidelity aerostructural optimization framework of complete aircraft configurations. Several optimization studies were   (Mavris and DeLaurentis 2000) conducted and optimal configurations of complete configuration flight and wind-tunnel models were obtained. Later on, Martins et al. (2004) proposed a coupled sensitivity analysis of aerostructural systems, demonstrating the accuracy and efficiency of the coupled-adjoint method by optimizing the aerostructural design of a supersonic business jet. Nearly concurrently, exploratory efforts on unconventional configurations commenced, with emphasis given on highaspect ratio wings due to their increased aerodynamic performance. Smith et al. (2001) underlined the necessity of high-fidelity CFD in the aeroelastic modeling of high-aspect ratio wings, since significant differences emerged in comparison to linearized methods. Particularly, linear methods tend to underestimate the tip bending and twist, leading to overly conservative divergence and flutter speeds predictions. Furthermore, in their work Patil et al. (2001) demonstrated a dramatic reduction in flutter speed due to the nonlinear coupling among edgewise bending and torsion as well as in the flight dynamic characteristics of an aircraft due to the wing flexibility, vastly affecting the trim solution. The differences between a linear and a nonlinear beam model of a high-aspect ratio wing coupled with RANS CFD at the transonic regime were also presented in Garcia (2005). The static aeroelastic response of unswept and swept configurations was investigated, indicating a reversal in wing twist due to nonlinear torsion-bending coupling effects for the former and a reduction in the amount of washout, in comparison with the linear solution, for the latter. Recently, Coggin et al. (2014) studied the aeroelastic response of a trussbraced wing configuration. The NASTRAN solver was modified in such manner that the nonlinear aeroelastic response of the wing was accounted for, with the resulting pre-stressed stiffness matrix being used to calculate the prestressed mode shapes which in turn were fed to the flutter analysis module. In Bartels et al. (2015), the flutter and limit cycle oscillation (LCO) of a similar configuration was investigated using low-fidelity panel methods as well as CFD simulations. Both studies indicate a discrepancy in the flutter speed calculation and LCO onset between the linear and nonlinear approaches, stressing the necessity for the inclusion of nonlinear structure and aerodynamic effects for the accurate prediction of the aeroelastic behavior of the trussbraced wing. In parallel, the development of high-fidelity aerostructural optimization tools, employing hundreds of variables and constraints was also of great interest to the scientific community. An approach towards the resolution of such vast optimization problems was conducted in Kenway and Martins (2014), where a multipoint high-fidelity aerostructural optimization of the common research model (CRM) was presented. The optimization problem was divided into two problems, aiming at the takeoff gross weight and fuel burn minimization under five cruise and two maneuver conditions respectively. The proposed framework managed to effectively balance aerodynamic and structural performance, despite the vastness of the design space and the constraints posed at each discipline. In particular, a 4.2% reduction in takeoff gross weight with a corresponding 6.6% fuel burn reduction were accomplished for the first of the two optimization problems, whereas for the second one a 11.2% fuel burn reduction with no significant change in the takeoff gross weight was observed. A similar study for the single and multiple point aerodynamic optimization of the CRM wing was also conducted in Lyu et al. (2015). The geometry of the wing was paremetrized using the free-form deformation (FFD) technique (Kenway et al. 2010), resulting into 720 variables. The RANS equations constituted the aerodynamics model, with the resulting aerodynamic shape optimization problem being handled via the MACH framework. The minimization of the drag coefficient subject to lift, pitching moment, and geometric constraints formed the objective function and constraints of the optimization problem respectively. For the single point optimization, a drag reduction of 8.5% was achieved, while a more robust design was obtained through a multipoint optimization. Various randomly generated starting points were also examined in order to assess the multimodality of the design space, with the resulting geometries pertaining similar geometrical characteristics. The effect of varying the number of shape design variables was also investigated. Various aeroelastic modeling techniques of high-aspect ratio wings are demonstrated in Howcroft et al. (2016) and Castellani et al. (2017). Specifically, the linear models tend to overestimate the wing deflection, bending moment and lift force. The large influence of the drag force, often disregarded in panel aerodynamics, on the increase of the resulting twisting moment is also highlighted. The research efforts of Calderon et al. (2018Calderon et al. ( , 2019 focused on the ramifications of including the geometrically nonlinear behavior of the structure on the structural optimization. Based on a coupled geometrically nonlinear VLM and beam models, a series of studies were conducted, obtaining optimum panels and stiffeners configurations under strength, skin buckling and skin-stringer failure constraints. For a baseline wing with an aspect-ratio of 18, a lighter configuration of 5% emerged for the geometrically nonlinear model, with the root bending and torsional moment being reduced by 8% and 50% respectively. A modest 1.3% improvement in the Breguet range was also noted. In the latter study, the effect of the aspect-ratio using linear and nonlinear models on the structural sizing was also investigated. Wings optimised using geometrically nonlinear analysis emerged lighter with similar aerodynamic efficiency, by means of lift-to-drag ratio, and therefore greater Breguet ranges than those optimized using a linear approach. Of particular interest is the fact that due to the overestimation of loads in the linear analysis the optimum aspect ratio increases from 18 to 19, when accounting for geometric nonlinearities. The inclusion of geometrically nonlinear flutter constraints was accomplished in Lupp and Cesnik (2019), where a blended wing body (BWB) aircraft was optimized under a bending curvature (static aeroelastic) constraint, a linearized about the jig shape as well as a geometrically nonlinear flutter constraint. The addition of the flutter constraint in the aerostructural optimization showed a reduced wing span and aspect ratio when compared to the bending curvature constraint. Introducing the geometrically nonlinear flutter constraint yielded a more conservative design with a higher fuel burn. The authors conclude that the necessity for a geometrically nonlinear flutter constraint is certainly problem and configuration dependent. Straying off the vastly used beam models and raising the computational fidelity in terms of structural representation, Verri et al. (2020) coupled a full-order geometrically nonlinear wing box FEM model of an Embraer regional jet with a wing aspect ratio of 12 to high-fidelity CFD. A 2.5 g pull-up maneuver limit load was considered and the differences between the linear and nonlinear structural behavior were underlined. In particular, and from the aerodynamics point of view, a modified pressure coefficient distribution acting on the wing due to the tip displacement difference of 16% was obtained. Furthermore, as also observed in the majority of the aforementioned studies, the lift force was shifted towards the wing tip direction due to the outboard rotation of the wing. The resulting internal shear loads remained unaltered in the root, but higher in the outboard portion by 14% for the nonlinear case. Regarding the bending moment, a slightly higher at the root but much elevated in the outboard portion was observed. A complete integrated aerostructural design and optimization effort was performed in Gray (2021), where two modified versions of the original CRM wing model were compared and optimized. Initially, and for a similar loading scenario, the nonlinear analysis resulted in greater bending stresses in the upper and lower skins which in turn increased the Von Mises and buckling failure criteria by around 10% for the highaspect ratio wing. The presence of the Brazier loads (1927) was also demonstrated, resulting in substantially greater compressive axial stresses in the ribs. Several structural and aeroelastic design optimisation studies were subsequently executed, with the increase in bending stresses leading to a 6% and 4% increase in mass for the optimised high and moderate aspect ratio wingbox respectively. A computational cost analysis was also included, indicating an order of magnitude increase when geometric nonlinearities are considered.
The majority of the aforementioned studies steered its attention mainly towards the identification of the influence of the geometrically nonlinear phenomena in the aerodynamic and structural disciplines, with any effort towards the implementation of composite materials in the analysis and optimization frameworks being nearly absent. On that front, a common way of optimizing composite structures involves the declaration of ply angles and thicknesses as design variables, resulting in a more complex and irregular design space, but with facilitated algorithmic applicability. On the other hand, optimization via lamination parameters, extensively studied over the past years, is presented as an alternative. Miki and Sugiyama (1993) was among the first to generate optimum designs of laminated composite plates via lamination parameters for maximum in-plane and bending stiffness, buckling strength and natural frequency. Fukunaga et al. (1994) explored the effect of bending-twisting coupling on the fundamental frequency of symmetric laminated plates cast in the lamination parameters space, indicating that this type of coupling reduces the fundamental frequencies. An optimal laminate configuration for maximum fundamental frequencies was also generated. The buckling load of a composite panel was maximized using flexural parameters in Liu et al. (2004). The resulting stacking sequence was compared to one generated via a genetic algorithm, indicating a close correlation between the two methods. In more recent works, efforts were directed towards the implementation of lamination parameters optimization frameworks for the aeroelastic tailoring of regular and variable stiffness composite materials wings (Thuwis et al. 2009;IJsselmuiden et al. 2010), as well as the stiffness optimization of composite wings subject to aeroelastic constraints (Dillinger et al. 2013). Blending constraints, guaranteeing a certain degree of ply continuity inside a variable stiffness and thickness composite wing were studied in Macquart et al. (2018) as well as Bordogna et al. (2020), extending the capabilities of the state-of-the art lamination parameters optimization algorithms.

Motivation
Despite the existence of a wide variety of low-fidelity aeroelastic optimization frameworks, studies towards high-fidelity optimization frameworks including geometric nonlinearities for composite materials high-aspect ratio wings are rarely documented in the literature. This knowledge gap is aggravated when examining sequential optimization frameworks, combining low and high-fidelity computational tools. Aiming to address this issue, this research study bridges the various fidelities of the numerical tools involved in order to provide an efficient sequential structural optimization framework of a high-aspect ratio composite materials aircraft wing subject to local panel buckling, strength and flutter constraints, combining low and high-fidelity computational tools. The geometric nonlinear behavior of the structure as well as the follower forces nature of the aerodynamic forces is accounted for at each stage of the optimization framework, since as demonstrated in the literature as well as in later sections they alter the deformation and stress field, which eventually translates into possible structural mass reduction.
Reference analyses are initially conducted in order to assess the effect of fidelity of the associated computational tools on the structural response of the reference wing, guiding the development of the proposed optimization approach. Low-fidelity tools form the backbone of this optimization framework and are initially tasked with guiding the structural design towards promising regions of the design space. Higher fidelity methods are then employed, aiming at exploring possible further gains in performance. The present approach intends to contribute to sequential structural optimization frameworks for composite materials future wing configurations, which are less prominently covered in the literature. We demonstrate that, through the use of this optimization approach, reasonable reduction in structural mass of the test case wing can be realized.

Methodology
In this section, the various variable fidelity computational tools used for the aerodynamics and structural disciplines are rigorously described. The methodology starts with the particulars of the high-fidelity CFD analysis and proceeds with the low-fidelity 3D panel and DLM method. Aspects of the nonlinear FEM analysis are also manifested, which will eventually lead to the core of this research work, the proposed sequential structural optimization framework.

Reference wing geometry
Within the context of this work, the undeflected CRM, namely uCRM − 13.5, has been chosen as the reference wing model. Constituting a modified, high aspect ratio derivative of the original CRM wing (Vassberg et al. 2008), the uCRM − 13.5 model (Brooks et al. 2018) serves as a benchmark configuration for CFD and aerostructural optimization studies of realistic, contemporary as well as future aircraft configurations operating at the transonic regime. The relevant geometric data are summarized in the following Table 1.
Regarding the internal configuration, two spars located at 10% and 60% of the local chord are present along with 54 evenly distributed ribs and upper and lower skin stiffeners. Spar and rib caps are also included withing the framework of this study, in contrast to the internal configuration presented in Brooks et al. (2018). These primary load-carrying components are commonly found in a commercial airliner's wingbox, hence their effect and influence in the overall stiffness and strength of the structure were deemed necessary to be investigated, as demonstrated in later sections. An exploded view of the internal structure is illustrated in the following Fig. 3. For the sake of clarity, the spar and rib caps as well as the stiffeners are depicted via spanwise straight lines.
For the subsequent analyses a 2.5 g pull-up maneuver limit load, assumed to be exerted at the conditions indicated in Table 2, is considered, with the angle of attack (AoA) being modified accordingly at each level of fidelity in order to generate the appropriate aerodynamic loading (with MTOW being the maximum take-off weight of the aircraft, set to 268e 3 kg Gray 2021).

High-Fidelity CFD
For the high-fidelity CFD analysis the compressible RANS equations (Crovato et al. 2020;Economon et al. 2016) are discretized and solved via the FVM in a C-grid shaped    Fig. 4 along with the relevant dimensions and boundary conditions. Turbulence in the form of RANS methods is also modeled and introduced via the Spalart-Allmaras one-equation turbulence model (1992) along with wall functions providing near wall treatment. A first cell wall distance, y + , necessary for capturing the evolution of the boundary layer using RANS turbulence models, in the region of 50-70 was targeted due to the size of the wing, constituting the use of lower values prohibitive in terms of computational cost. Given the target y + value, as well as the Reynolds number of the analysis, which in our case is equal to 35.5e 6 , the first layer height is calculated based on equations derived from analyses of the boundary layer development in flat plates, as described in White (2010). Moving on, and to account for temperature induced changes in viscosity, Sutherland's law, based on the kinetic theory of ideal gases and an idealized intermolecular-force potential, is used. The spatial terms of the Navier-Stokes equations are discretized on the computational finite volume mesh with appropriate discretization schemes. In particular, the spatial convection terms are discretized using second order upwind schemes, mitigating numerical oscillations while maintaining good accuracy, while spatial diffusion and source terms are discretized using central differences. Moreover, a Green-Gauss method has been chosen for the evaluation of gradients. The resulting structured hexahedral CFD mesh is shown in Fig. 5. For verification purposes, a mesh convergence study for the major aerodynamic lift, drag and pitching moment coefficients, C L , C D and C M respectively, has also been conducted and is presented in the subsequent Table 3. The resulting aerodynamic loads are transferred to the structural FEM mesh via interpolation techniques, as discussed in Kilimtzidis (2022) at a greater extent.

3D panel method
On the medium to low-fidelity aerodynamic analysis front, an in-house 3D panel method framework has been developed in MATLAB. This module accepts a discretized geometry of an aircraft wing along with the reference flight conditions and calculates several aerodynamic related quantities. In particular, and under the assumption of an inviscid and irrotational flow, there exists a potential function that satisfies the Laplacian equation (Katz and Plotkin 2001): Applying Green's third identity and the impermeability boundary condition on the wing's surface, while introducing sources and doublets as elementary solutions, denoted as and p respectively, Eq. 1 is recast in the following form: With n being the surface normal vector, r the distance vector and ∞ the free-stream potential function. For a discretized surface consisting of N wing and N w wake panels, integration where A 1k and B 1k doublet and source influence coefficient for the kth panel respectively. Setting the source terms on the right-hand-side of Eq. 3 equal to k = n k ⋅ U ∞ , a set of linear equations that can be solved for the unknown doublet distribution is obtained. Making use of spatial interpolation schemes one can calculate the induced velocities and via the Bernoulli equation the pressure coefficients on each panel. To account for possible compressibility effects in the flow, the Prandtl-Glauert correction factor has been also implemented, as in Katz and Plotkin (2001). Additionally, the wake of the wing has been considered to be fixed and extending in the chordwise direction, since wake shape calculations are beyond the scope of this research study. In terms of load transfer to the structural FEM mesh and for computational efficiency purposes, the resulting aerodynamic panels and structural meshes are coincident, allowing for a direct load transfer between the two disciplines. The results of the 3D panel method developed within the framework of this study were also verified and validated against the ONERA M6 wing case study (Schmitt and Charpin 1979) and are presented in Kilimtzidis (2022). The resulting aerodynamic mesh is displayed in the following Fig. 6.

DLM aerodynamics
Three-dimensional panel methods are particularly advantageous in problems where thick airfoil configurations exist. When the airfoil thickness is adequately small (less than 12%), further simplifications can be applied, reducing the dimensionality of the aerodynamic surface to planar representations. Among of the existing panel methods, the DLM constitutes the cornerstone for aerodynamic calculations at this level of fidelity. A comprehensive work on the mathematical derivation of the method is provided in Blair (1992). Aerodynamic surfaces are typically divided into small trapezoidal panels and a constant pressure distribution is assumed. The DLM is capable of modeling unsteady flows in the frequency domain and is based on the linearized small disturbance potential flow equation, hence neglecting large perturbations and shockwaves: where M ∞ and U ∞ the free-stream Mach number and velocity respectively and xx , yy , zz , xt , tt second derivatives of the potential function. Assuming small amplitude, harmonic motion, solution to Eq. 4 is given by the so-called acceleration potential and by introduction of doublets across the discretized planar lifting surface. The doublet line of constant value is placed at the 1/4th of the chord of each panel, while the normalwash collocation point, where the impermeability boundary condition is applied, is placed at the 3/4th of the chord. The resulting pressure coefficient for a panel i, denoted as p i , is expressed in terms of a complex matrix A ij (Albano and Rodden 1969), including the contribution of the generated downwash from panel j to panel i, and the downwash value w i at the current panel of interest: Complex valued matrix A ij , widely termed as the aerodynamic influence coefficient (AIC) matrix, is now subject to evaluation given the discretized aerodynamic surface, the operating Mach number as well as the reduced frequency k.
The DLM aerodynamic surfaces of the uCRM − 13.5 wing are modeled in NASTRAN, with the resulting aerodynamic mesh being presented in Fig. 7 along with the calculated collocation points. Despite the attractive characteristics of the DLM aerodynamic modeling, certain corrections need to be applied in order to replicate the actual airfoil geometry to a greater extent. Inclusion of the realistic camber line and twist angle of the corresponding airfoil, allowing for a proper calculation of the aerodynamic coefficients of the wing at various angles of attack is often sought. This is accomplished by populating the W2GJ entry vector in NASTRAN with the values of the slope of the camber line of each airfoil of the wing under consideration calculated at the collocation points of the aerodynamic panels. For the calculation of the relevant correction factors for the uCRM − 13.5 wing, a two-stage procedure, similar to the one presented in Demirer (2021) has been implemented herein: The resulting W2GJ distribution for the current aerodynamic mesh is shown in the subsequent Fig. 7. Regarding the coupling between the aerodynamic and structural mesh, infinite plate splines have been used in this study to transfer the aerodynamic loads to the structure. For the splining procedure, nodes close to the ribs and spars as well as to the DLM mesh have been chosen, allowing for a proper load distribution and avoiding possible local loading scenarios. For the sake of completeness, a verification and validation analysis of the DLM aerodynamic model as well of the flutter analysis has been conducted (Kilimtzidis et al. 2018).

FEM model development
The well-established and reliable FEM analysis has been vastly used among researchers and the aerospace industry to model structures and to predict their behavior under the application of loads. The accuracy of the method for wing structures has been validated against experimental data in a plethora of literature studies (Ritter et al. 2021;Dessena et al. 2022;Keimer et al. 2022 where t the thickness,A 11 , A 12 , A 2 and A 66 the corresponding terms of the extensional stiffness matrix of a laminate. Regarding the boundary conditions, the wing is assumed to be clamped at its root, thus fixing all relative nodal degrees of freedom (DoF). For the sake of completeness, external Table 4 Composite materials properties mass as well as gravitational loads have been accounted for. In particular, a fuel load of 56.000 kg (Brooks et al. 2018) has been considered as traction loads and introduced in the model at each wingbox bay via RBE2 elements acting on the lower skin. Furthermore, the engine mass, equal to 2000 kg, is also modeled and connected to the wingbox via a combination of concentrated mass and RBE2 and RBE3 elements. The resulting FEM model of the uCRM − 13.5 wing is shown in the subsequent Fig. 8. The NASTRAN solvers SOL 101 and SOL 400 (Lee 1992) were chosen to carry out the linear and nonlinear static analysis respectively. In order to avoid any possible divergence issues mainly during the optimization procedure, a full Newton-Raphson stiffness matrix update technique has been preferred for the nonlinear analyses.

Global-local panel buckling analysis
Due to the criticality of the buckling phenomenon and its impact on the design of aircraft wings, a global-local FEM modeling technique was implemented for the panels of the upper and lower skins within the framework of this study. A panel is considered as the intersection between the front and rear spar as well as two adjacent ribs, as indicated in Fig. 9. As far as the global-local modeling technique is concerned, the resulting nodal displacements at each of the four edges of each panel of the global FEM act as input boundary conditions to the local level FEM analysis. Since dissimilar meshes are considered, interpolation schemes are employed at each edge of each panel in order to calculate the local displacement field at intermediate positions. Furthermore, and at this level, blade stiffeners are considered and modeled explicitly via CQUAD4 elements. Buckling analyses are executed for each of the 106 panels present in the wing via the SOL 105 solver implemented in NASTRAN, with the first eigenvalue, namely the critical buckling load expressed as a percentage of the originally applied load, being the quantity of interest (NASTRAN 2021).

Reference analyses
Prior to the optimization study, reference analyses are conducted in order to ascertain the differences among each level of fidelity on the structural response of the uCRM − 13.5 wing. Various levels of fidelity have been considered from an aerodynamics and structural point of view, with the relevant analyses conducted being listed in Table 5.
In all types of analyses, the AoA is modified accordingly in order to attain the required loading conditions, as described in Table 2 and listed in Table 6.
Significant changes in AoA are required for the various aerodynamic analysis tools in order achieve a similar lift load. For the 3D panel method a smaller angle is required since pressure recirculation as well as flow separation near the trailing edge phenomena are not captured in such   DLM elastic Linear methods, resulting in higher pressure difference at the outer regions of the wing. It should be point out, nevertheless, that a flat wake has been assumed throughout the analysis. Wake shape calculation techniques are reported in Katz and Plotkin (2001) and influence the calculation of the aerodynamic coefficients. For the DLM, the camber correction method presented earlier is not sufficient to obtain the required lift load at the high-fidelity AoA, the main culprit being thickness effects, which are not captured in this method, that allow for a greater pressure difference between the upper and lower surfaces of each wing section. The induced displacement field in terms of deflection and tip torsion angle, the static strength as well as the buckling critical loads for the upper and lower skins constitute the quantities of interest for the analyses. On the static strength front, a first-ply failure (FPF) via the Tsai-Wu criterion (Jones 2018;Tsai and Hahn 2018) is used for the strength prediction in terms of the failure index (FI) of the composite skins, spar and rib webs. For the rest of the parts, a direct comparison between the maximum stress and the corresponding material strength value is made in order to obtain the FI for the wing component under consideration. As far as the buckling critical load is concerned, buckling analyses are executed for each panel of the upper and lower skin with the first buckling eigenvalue (BE) being extracted. The minimum eigenvalue should then be greater than unity in order to avoid any buckling failure. 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 (Kreisselmeier and Steinhauser 1979;Poon and Martins 2006) for the static strength as well as buckling constraints. For a problem consisting of N c quantities of interest, g, having a maximum value of g max , the KS functions are of the following form: and are formed separately for the calculated FI of each of the components involved in the static strength as well as for each upper panel buckling critical load 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 (Lambe et al. 2016;Lambe and Martins 2015) and is shown to provide accuracy of the optimal solution.
Substantial differences in the results can be observed depending on the level of fidelity of the aerodynamics as well as of the structural representation of the uCRM − 13.5 wing configuration. For the former, and despite acquiring similar loading conditions, the 3D panel method seems to be overestimating the maximum deflection and tip torsion induced at the wing. Pressure recirculation as well as flow separation near the trailing edge effects are not captured in such methods as stated earlier, resulting in higher pressure difference and eventually loads at the outer regions of the wing, which in turn are responsible for this increase. Higher bending moments induce a higher stress field and lower the critical buckling eigenvalues, as seen in Table 7. On the contrary, the DLM seems to be more accurately predicting the displacement and stress field when compared to the highfidelity CFD linear analysis, as indicated in Fig. 10a and b. Of particular interest is the aeroelastic behavior of the uCRM − 13.5 wing, allowing for aerodynamic load redistribution and overall lower displacement and stress fields due to the bending-torsion coupling phenomenon. From a structural point of view, inclusion of nonlinearities affect the solution and result in lower displacements and stresses as expected mainly due to the geometric stiffening of the structure. This effect is on the contrary less pronounced for the torsion angle, as depicted in Fig. 10c. Other nonlinear phenomena also rise, with one of the most prominent being the tip shortening effect, as illustrated in Fig. 10b. Moreover, the presence of Brazier loads drastically alters the stress field in the rib caps, as demonstrated in Table 7, resulting in highly elevated stresses by nearly 50% in comparison with the linear models.
The distribution of the maximum FI per ply for the linear and nonlinear structural analyses under the CFD loading, along with the relative differences is presented in Fig. 11a. The resulting distributions are quite similar, with the majority of the differences lying near the Yehudi break of the wing. In particular, the reduction of the effective span of the wing as well as the rotation of the outboard sections of the wing reduce the moment arm and the bending moment, which cascade in lower stress values near this spanwise region. A similar situation for the 3D panels method is portrayed in Fig. 11b. In terms of the effect of the aerodynamic loading on the distribution of the FI, the 3D panel method overpredicts the induced root bending moment and as a result the developed stress field, as indicated in Fig. 11c.
Regarding the computational time, a wall-clock time comparison for the linear and nonlinear solvers has been conducted for the reference analyses presented. In particular, the wall-clock time was 5.68 min and 17.8 min for the linear and nonlinear solver respectively, indicating the superiority of the linear solver in terms of computational efficiency. In terms of convergence, nonlinear analyses at specific points of the design space indicated divergence issues for the spar caps CBEAM elements near the yehudi break of the wing at certain points of the design space, even for a full Newton-Raphson stiffness matrix update technique. As a result, the lower bounds of those elements were increased in order to avoid divergence issues during the optimization process.

Optimization framework
The particulars of the optimization framework are discussed in detailed fashion in this section, starting from the definition of the variables and proceeding with the formulation of the optimization problem.

Variables definition
The ply count of the 0 • , 90 • and 45 • ∕ − 45 • plies of each laminate and for each wing component constitute the variables of the optimization problem. Setting the number of 45 • and −45 • plies equal implies generation of symmetric and balanced lay-ups, thus achieving conformity to composite materials design guidelines (Kassapoglou 2013). Within the framework of this optimization study, the wing is also divided into eight spanwise and evenly spaced zones, as illustrated in Fig. 12, allowing for a wider design space as well as increased structural design freedom. A summary of the variables present along with the respective lower and upper bounds as well as their type is provided in Table 8. It should be pointed out that the two materials presented in Table 4 are also used as variables and can be assigned during the optimization process to the rib webs and caps, spar webs and caps, upper and lower skin as well as the stringers. This particular parametrization allows for a rapid cost estimation for each part of the wingbox and provides the designed with greater freedom with respect to the material of each component.

Objective function and constraints
In our study the minimization of the mass of the wing represents the objective function. Design constraints in terms of static strength, buckling and flutter velocity, formulate the optimization problem. Since no explicit static stiffness requirements do exist in the literature for this generation of aircraft wings, such constraints are not imposed. Furthermore, the static strength and critical buckling load of the upper and lower skin panels for each candidate design are also included and expressed in terms of the KS functions, as described in earlier sections. Completing the set of constraints, the dynamic aeroelastic instability by means of flutter velocity of the candidate design solutions is also investigated. In particular, the p−k method implemented in NASTRAN SOL 145 is used to identify any possible divergence and flutter instability that might be present. For the purpose of the analysis, a set of 20 reduced frequencies are also used to calculate and interpolate in user-defined velocities the reduced-frequency-dependent aerodynamic loads. Possible flutter instabilities are investigated for the first ten structural modes via the velocity-damping ( V−g ) and velocity-frequency ( V−f ) plots. In particular, the trends of each mode are monitored regarding the damping values, with positive damping indicating a possible flutter instability. The   corresponding flutter velocity is then calculated via linear interpolation between the previous and subsequent velocities and damping values. In case of no flutter point, the flutter velocity is set to a value outside of the velocity range of the analysis. Classical infinite plate splines, similar to the ones described earlier, have been used to transfer the aerodynamic loads to the structure. The objective function as well as the constraints are summarized in the following Table 9.

Optimization algorithm parameters
Within the framework of our study, the MIDACO solver (Schlüter et al. 2013) has been chosen to carry out the optimization problem. MIDACO adopts a combination of an extended ant colony optimization algorithm (ACO; ) along with the Oracle Penalty Method (Schlüter and Gerdts 2009), an advanced method developed for metaheuristic search algorithms for constraint handling of the solution process. Despite the abundance of optimization algorithms in the literature, the nature of the optimization problem, featuring black-box objective function and constraints along with discrete variables, as described in Sect. 4.2, limit the applicability of many of the optimization algorithms present, such as gradient-based ones. Despite their reduced efficiency with increased number of variables, gradient-free algorithms (e.g. simulated annealing, genetic algorithms) tend to make less assumptions about the modality and smoothness of the design space, thus present increased robustness characteristics (Martins and Ning 2021). Among these algorithms, the MIDACO algorithm has proven its efficiency and accuracy in a plethora of mathematical benchmark optimization problems as well as in engineering applications (Schlueter 2014;Kontogiannis and Savill 2020). Regarding the solution procedure and in order to increase the effectiveness of the optimization problem, a sequential approach with multiple runs has been adopted. On that end, the predefined number of executions is divided into runs pertaining different algorithmic parameters as well as level of fidelity. Particularly, initial runs are executed via the 3D panel method and mainly focus on extensive design space exploration, bounding the structural design towards promising regions of attraction. From a structural point of view, and based on the results presented earlier, the inclusion of nonlinearities even at this design stage strongly affects the structural behavior of the reference wing and are therefore deemed necessary throughout the design stages. Early runs are also accompanied by a relaxed constraint satisfaction tolerance. As the solution advances, higher-fidelity CFD solutions are employed, while the search becomes increasingly local. This is accomplished by modifying 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 multi-kernel Gaussian probability density functions (PDF). For a generic variable k with upper and lower bounds x u and x l respectively, the FOCUS parameter applies an upper bound for the standard deviation of a Gaussian PDF given by 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, with this procedure being repeated for a predefined number of iterations, satisfying user-defined stopping criteria. The number of iterations was defined based on the computational resources available, while the associated parameters adopted for the current optimization problem were selected based on a similar study (Kilimtzidis et al. 2023). The cascading runs approach along with certain algorithmic values is also suggested in MIDACO (2018). Overall, the optimization algorithm parameters are summarized in Table 10: A flowchart of the proposed optimization framework is also provided in Fig. 13.

Results
The convergence of the objective function for the sequential optimization approach is presented in Fig. 14. Convergence to a minimum mass, along with non-violating design constraints, has been achieved. In particular, significant mass reduction is obtained by employing high-fidelity CFD solutions, since the corresponding loading condition results into lower stress levels, which eventually translate into reduced    Table 11. The resulting mass distribution for each wing component, along with the corresponding material, is provided in Fig. 15. The PW Fabric material is assigned to the majority of the parts, resulting in lower wing mass due to its lower density value. On the other hand, and despite the associated higher density, the Hexcel IM7 UD is assigned to the lower skin due its superior performance in tension in comparison with the PW Fabric. On the constraints front, the resulting optimal design is clearly stiffness driven, since the aggregate minimum buckling eigenvalue for the upper and lower skin closely approximate their limit value. As illustrated in Fig. 16, the critical panel buckling load decreases in the spanwise direction, indicating that the outboard panels are more prone to buckling. On the other hand, strength constraints are less tight, indicating possible further mass reduction. However, additional strength safety factors are commonly introduced to account for various stress risers, such as manufacturing defects, holes, impact phenomena, etc. The developed stress field, by means of the maximum FI as well as the deflected wing are depicted in Figs. 17 and 18.
The thickness distribution among the wing parts is provided in Table 12 and also visualized in Fig. 19. As a general trend, a spanwise thickness decrease can be observed for the majority of the wing components, specifically occurring at the outboard section of the wing, since the bending moment is decreased towards the wing tip. The region around the yehudi break, (Zones 3 and 4) appears to be the highest stressed one in contrast with the root of wing, and is therefore is accommodated with elevated thickness values. On a component-wise analysis, the spar caps and skins obtain the highest thickness values, increasing the secondary moment of inertia of the cross-sections, thus allowing for greater stiffness of the wing. On the contrary, the front and rear spar webs are assigned with lower thickness values, since they are closer to the neutral axis of bending of the wing resulting in lower stresses. Shear forces and induced torsional moments, however, highlight the need for elevated thicknesses even for these components. Of particular interest  are the thickness values obtained by the rib webs and caps, which are mainly driven by the Brazier loading. The former obtain the majority of the aforementioned loading and therefore an increase in their thickness is required, while the latter appear to obtain low thickness values, transferring the developed loads to the webs. As far as the stringers is concerned, relatively high thickness values throughout the span of the wing are deemed necessary in order to sustain the critical buckling load.
The percentages of the 0 • , 90 • and 45 • ∕ − 45 • plies of each laminate and zone are also presented in Fig. 20. The majority of the components, and specifically the front and plies mainly due to their initial percentage in the baseline lay-up. 0 • plies are present in the skins and caps to accommodate for the presence of normal stresses in these components, while on the other hand 90 • plies, enhance the strength of the components at the normal to the fiber direction.
Regarding the flutter analysis, the V−g and V−f plots are presented in the following Fig. 21. Clearly, no indication of instabilities at each of the first ten structural modes is observed for the optimal design and the velocities range of the flutter analysis.
The comparison of the obtained optimal mass of the wing with similar results present in the literature is also discussed. In particular, upon removal of the external as well as the leading and trailing edge skin masses, which are not included in the literature, the optimized mass of the wing is 29.247 kg, while a mass of 30.032 kg is obtained in Brooks et al. (2018). In Brooks et al. (2020), the mass of the optimized composite materials wing ranges between 20.000 and 40.000 kg, indicating good agreement between the optimized mass values. A similar study was conducted in Brooks et al. (2019), with the optimized mass being 19.796 kg. The results of the aforementioned studies differs to the one presented herein in the sense that combined aerostructural tools and measures such as fuel burn were considered for the optimization cases. The geometry of the wing was also parametrized in terms of twist, which produced an optimized lift loading distribution, decreasing the bending moment by decreasing the lift produced in the outer sections of the wing. The results presented in the optimization study conducted of Gray (2021) in terms of structural mass can be compared more accurately to the ones of the present study, since no aerostructural performance optimization was conducted and the nonlinear response of the structure was also accounted for. The resulting 2.5 g mass was around 14.000 kg, with the major difference lying in the mass of the stringers. In our case, buckling was deemed critical in the current configuration. As a result, and in order to support the developed buckling loads, the thickness of the stringers was directed towards greater thickness values, as also demonstrated in Fig. 19. This phenomenon along with the predefined width and spacing of the stiffeners led to their elevated total mass. Another important issue is the fact that the approach towards buckling analysis is different in the two studies, with the study in Gray (2021) following a smeared modeling approach along with a global buckling analysis. In our case, nevertheless, the stiffeners have been modeled explicitly and buckling was also examined through a more detailed, global-local modeling approach, which could have eventually shed light into failure mechanisms that cannot be predicted by a combination of a smeared modeling and global buckling analysis approach. Nevertheless, different spacing and width values should be explored in future frameworks in order to assess their effect at a greater extent.

Conclusion
A sequential optimization framework combining variable fidelity computational tools for future high-aspect ratio composite aircraft wings, and specifically the undeflected, uCRM − 13.5 wing, based on the MIDACO algorithm has been presented herein. The wing structure was represented using a 3D shell and beam elements FEM model, while aerodynamics were treated via RANS, 3D panel as well as DLM solutions, depending on level of fidelity. The latter provided also insight on the static aeroelastic behavior of the reference wing. Comparison between the different fidelities was initially performed with the aim of highlighting their effect on the structural response of the test case wing. For the aerodynamics, significant changes to the AoA to attain the target critical aerodynamic loading scenario compared to the high-fidelity CFD were observed. Furthermore, the developed 3D panel method resulted into higher vertical deflection and torsional angles in comparison with the high-fidelity CFD, while the corrected DLM provided a closer approximation. From a structural point of Fig. 21 Optimal design-V−g and V−f plots view, the inclusion of follower forces and geometric nonlinearities greatly altered the displacements field, inducing lower vertical deflections yet higher tip shortening. The aforementioned models were combined and cast in a sequential optimization framework aiming at obtaining a minimum mass configuration subject to strength, stiffness in terms of panel buckling as well as flutter constraints. The framework was initially guided by low-fidelity aerodynamics into promising basins of the design space, with high-fidelity CFD runs providing refinement of the optimal solution. The resulting wing structure was critical to panel buckling, with strength criteria being more relaxed. Zones near the yehudi break, and specifically skins and spar caps, obtain higher thickness values in order to provide the sufficient stiffness and strength for the structure.
In terms of future work, the present framework could be benefited by the addition of aerodynamic shape optimization procedures, which due to the lack of computational tools and limited pre-processing capabilities were not accounted for in this study. Geometry parametrization techniques, such as the FFD are capable of altering geometric variables of the wing, such as the twist, sweep, span etc. during the optimization stages, thus enabling the optimization of combined measures of aerostructural performance. When cast into a sequential optimization approach, additional effort needs to be directed towards analyzing the effect of low-fidelity tools on drag predictions, since they can significantly affect the overall performance. Turning our attention towards composite materials, the inclusion of a library of materials and lay-ups in order to support decision making in manufacturing and cost analysis at a greater extent as well as a wider variety of global-local analysis techniques (spars buckling, presence of holes and rivets, etc.) could further expand the capabilities of the present framework. In order to demonstrate the robustness of the optimization algorithm, several starting points could be provided in subsequent studies, comparing the optimal solutions obtained. The algorithmic parameters, as presented in Table 10, should also be further evaluated in order to shed more light on their effect in the sequential optimization framework. The ratio of the low to highfidelity runs as well as the effect and potential benefit from early on high-fidelity CFD samples in the optimization framework should be also investigated.