Towards Structural and Aeroelastic Similarity in Scaled Wing Models: Development of an Aeroelastic Optimization Framework

: The pursuit of more efficient transport has led engineers to develop a wide variety of aircraft configurations with the aim of reducing fuel consumption and emissions. However, these innovative designs introduce significant aeroelastic couplings that can potentially lead to structural failure. Consequently, aeroelastic analysis and optimization have become an integral part of modern aircraft design. In addition, aeroelastic testing of scaled models is a critical phase in aircraft development, requiring the accurate prediction of aeroelastic behavior during scaled model construction to reduce costs and mitigate the risks associated with full-scale flight testing. Achieving a high degree of similarity between the stiffness, mass distribution and flow field characteristics of scaled models and their full-scale counterparts is of paramount importance. However, achieving similarity is not always straightforward due to the variety of configurations of modern lightweight aircraft, as identical geometry cannot always be directly scaled down. This configuration diversity has a direct impact on the aeroelastic response, necessitating the use of computational aeroelasticity tools and optimization algorithms. This paper presents the development of an aeroelastic scaling framework using multidisciplinary optimization. Specifically, a parametric Finite Element Model (FEM) of the wing is created, incorporating the parameterization of both thickness and geometry, primarily using shell elements. Aerodynamic loads are calculated using the Doublet Lattice Method (DLM) employing twist and camber correction factors, and aeroelastic coupling is established using infinite plate splines. The aeroelastic model is then integrated within an Ant Colony Optimization (ACO) algorithm to achieve static and dynamic similarity between the scaled model and the reference wing. A notable contribution of this work is the incorporation of internal geometry parameterization into the framework, increasing its versatility and effectiveness.


Introduction 1.Background
Nowadays, the continuous search for climate-friendly transportation has led engineers to explore in-depth the design space and achieve more efficient configurations in order to reduce emissions since aviation has a significant impact on climate change through the release of nitrogen oxides, water vapor, sulfate and soot particles at high altitudes [1].The current objective is to reduce emissions by at least 55% by 2030 compared to 1990 in order to become climate neutral by 2050 [2].So, ensuring sustainable and efficient aviation is imperative, with the aerospace scientific community striving towards the fulfillment of these goals.A way, among others, to reduce emissions is to design aircraft with a high Aspect Ratio (AR).The Induced Drag Coefficient, CDi, is inversely proportional to AR, hence, wings designed with the highest possible AR enable induced drag reduction, which in turn cascades into less fuel consumption [3].This trend is clearly illustrated in the subsequent Figure 1, where an increase in the AR of transport airliner wings is observed in the past decade.Increasing AR is a great way to improve aerodynamic efficiency.However, from the structural discipline perspective, the design of such structures proves to be a challenging task.In fact, high AR wings are more sensitive in bending and torsion due to their long span, thus, introducing non-linearities and aeroelastic effects.Aeroelasticity is now seen as one of the key disciplines in design and as one of the "critical" processes in the aircraft development approach [5].Because of the significance of those phenomena, a plethora of methods to predict them has been developed, including the development of computational and analytical methods, wind tunnel experiments on aeroelastically scaled models and flight testing of full-scale models.Despite the prevalence of computational methods, aeroelastic scaling is still a vital part of the design process and is mainly used for the in-flight testing of new designs, in order to reduce the development costs, time and risks of manufacturing and testing the actual full-scale aircraft [6].
In order to achieve an aeroelastically scaled model, similarity in the three main properties of the aircraft is typically sought: Flow similarity (e.g., external shape and flow conditions matching) Depending on the particular test, some of those properties must be reproduced.In the case of a divergence or a control effectiveness test, the structural stiffness, distribution and flow must be similar.For a vibration test, structural stiffness and mass distribution are important, while flow does not need to be considered.Instead, the flutter test must have all three properties, and for this reason, the flutter model is the most difficult to capture.Note, that only if a model is similar to the full scale in all those three properties, it is considered as aeroelastically scaled.The easiest way to produce a flutter model is to directly scale down the exact geometry of the aircraft's structure to achieve structural similarity.Thus, all of the components must produce a similar stiffness and mass distribution along with the external shape.However, in most cases, this is impossible to achieve because modern aircraft are already lightweight structures with very thin components and the scaled-down thickness distribution might not be realistic.As a result, the design of different configurations with manufacturable thickness distributions or different materials is explored, and thus, the problem of similarity maximization is introduced (scaling).
In terms of computational tools, and due to the rapid growth of computing power in the last 2-3 decades, traditional aeroelasticity studies have changed to methods of Computational Aided Engineering (CAE) using the coupling of Computational Fluid Dynamics (CFD) (including lower fidelity methods) solvers with Computational Solid Mechanics (CSM) solvers.The accuracy of the solvers used is based on some of the parameters below [7]: • Project phase (e.g., feasibility phase, preliminary design, detailed design, certification) The use of graded fidelity tools in aeroelastic computations has made it possible to minimize the associated costs.The introduction of low-fidelity solutions in the early stages of aircraft design enables the exploration of the design space due to the reduced computing time.They also provide valuable information for the next design stages and thus reduce project development costs and time.Increasing the fidelity of the solvers is a necessary step as the design becomes more detailed.The use of low-fidelity models in preliminary does not capture specific physical effects that may affect the aircraft's flight behavior, and therefore, introducing higher fidelity tools is essential as the design goes on.Also, highfidelity solvers are established in the later stages of aircraft design to identify and confirm wind tunnel test data, explore the greatest potential, and eliminate uncertainties [7].A review of the main aerodynamic and structural computational models used within the aerospace community is provided in [7].

Motivation
As previously stated, the aeroelastic testing of scaled models is an indispensable part of the aircraft design cycle.Due to the increasing structural flexibility of contemporary aircraft, dynamic stability and automatic control are affected.As a result, producing scaled models could possibly alleviate this problem.Traditionally, the procedure of scaling was treated by non-dimensionalizing the model's governing equations using analytical approaches.The aeroelastic scaling procedures evolved with the rapid development of CAE methods and the modern procedures and practices consist of optimizing the target model to meet the full-scale aircraft's response.The use of FEM coupled with optimization algorithms for the scaling problem including only stiffness distribution was first suggested in [8].Results were validated with holography experimental data and were very promising for the method.The aforementioned framework was also extended to take into account the mass distribution and accurately predict the flutter velocity and frequency by relatively simple structures ( [9]).Aeroelastic scaling methodologies for the design and fabrication of an aeroelastic scaled RPV of a joined wing SensorCraft were proposed and developed in [10].The lightweight and high AR design introduces non-linearities and a properly scaled model is required to explore those phenomena.The first method used one optimization routine for the direct scaling of the modal response.The second method implemented two separate optimization routines, stiffness matching, and mass matching.Both methods achieved results with adequate similarity but the second one proved to have less computing cost.A new method based on previous work that uses a non-linear aeroelastic analysis was proposed in [11] .Their method was implemented in a simultaneous loop of static and modal response optimization using a gradient-based optimization tool.Although the flutter speed was higher, they achieved a better matching of the first four frequencies to the target model than the traditional method.A framework for non-linear aeroelastic scaling with applications in a HALE aircraft was also developed in [12].Linear scaling factors were extended to consider the case of aeroelastic scaling with geometric non-linearity.An approach for aeroelastic scaling using gradient-based optimization to match the vibration and buckling modes and eigenvalues, in addition to the static aeroelastic response was proposed in [13].For the optimization, they used a multi-start strategy to explore the design space.PolyJet 3D printing technology in order to manufacture a scaled model without the undesirable skin buckling that occurred in previous models they manufactured was used in [14].Turning our attention to high-aspect ratio wings, a framework for the aeroelastic scaling by serial and parallel optimization of the static and modal response was developed in [15].Additionally, a non-linear aeroelastic scaling procedure was proposed, with the components' thicknesses being considered as the design variables.Furthermore, a methodology for designing, manufacturing and testing dynamically scaled wings was proposed in [16].Various design and manufacturing aspects were discussed and analyzed, while normal modes for both models were obtained both numerically, as well as experimentally and compared via the Modal Assurance Criterion (MAC).Furthermore, the aeroelastic scaling of a reduced scaled wing model designed with three different scaling sets of primary quantities was investigated in [17].The results in terms of aerodynamics, structural dynamics as well as the nonlinear aeroelastic response of the scaled models were compared with the one estimated for the full-size model.The prevalence of composite materials in the aeronautics scene also led the scientific community towards the exploration of aeroelastic scaling frameworks of composite materials wings [6,18].More recently, an MDO methodology for aeroelastic scaling without similar flow conditions was proposed in [19] .They used PANAIR for aerodynamics and NASTRAN for structural analysis with shape parameters such as sweep and local chords and structural parameters such as shell thicknesses, and stringers cross sections.They used COBYLA, a gradient-free optimization algorithm, and achieved an in-flight similarity error of 1%.The presented methodology was lately expanded to include flutter matching and the mode tracking technique [20,21].
Table 1 summarizes the features of some contemporary studies on aeroelastic scaling.The features are mainly based on model fidelity (or complexity), analysis complexity (e.g., linear, or non-linear) and optimization complexity (e.g., number of parameters).It is obvious that various types of work have been conducted during the past years, each of them focusing on one or more aspects of optimization.From the analysis perspective, many medium and high fidelity methods (mainly on the structural analysis) have been implemented and various types of optimization techniques have been used.Nevertheless, a major gap in the literature exists in geometry parameterization.In fact, thickness parameterization is an inextricable part of scaling but, studies with geometric parameterization of wing components have not yet been conducted.Therefore, this study is aiming towards: 1.
Development of an automated aeroelastic analysis framework that utilizes a parameterized internal structure of the wing.

2.
Incorporation of this analysis framework in an aeroelastic scaling algorithm.

3.
Investigate how beneficial this parameterization could be in aeroelastic scaling.

4.
Evaluation of the computational and development costs of the algorithm in relation to the thickness optimization of a given geometry.
As stated in [22], an optimization problem cannot achieve maximum complexity in more than two disciplines (model, analysis, or optimization complexity).So, in this research, emphasis will be placed on the complexity of optimization and the other two areas will be limited from moderate to low complexity.

Materials and Methods
In this section, the cornerstones of this research are presented individually.The theoretical background for the aeroelastic scaling analysis, based on the linear equations of aeroelasticity, is well-documented in the literature and follows the methodology presented in [11,23,24].Starting from the structural FEM model, we then proceed to the particulars of the DLM aeroelastic model and, finally, manifest the overall aeroelastic scaling framework.

Structural FEM Model
As stated earlier, the structural analysis is an inextricable part of the aeroelastic analysis.The method that will be used in this framework in order to compute the structural response of the wing is the FEM.FEM is a numerical method used to solve ODE or PDE governing a physical problem.For problems with complex geometry, it is usually impossible to apply an analytical method; hence, the solution needs to be based on weak formulation numerical methods such as the FEM.Unlike analytical methods, FEM implements a discretization method in the body and the governing equations are applied to each of the finite elements leading to a system of algebraic equations instead of differential equations.In terms of modeling, commonly thin-walled aircraft structures are modeled via a combination of shell and beam elements.The skins, spar and rib webs are modeled via linear, 4-node elements (QUAD4) and linear, 3-node elements (TRIA3) shell elements.On the other hand, spar caps and stringers are modeled via beam elements employing the Timoshenko beam theory.A typical structural mesh of the uCRM wing, developed in NASTRAN, is presented in Figure 2.

DLM Aerodynamics Model
The aerodynamic loads used during the design process are calculated using an aeroelastic quasi-stationary trim analysis.The multidisciplinary simulation utilizes the DLM [25], a potential-theoretical aerodynamic method that can be used to determine motion-induced air forces in the subsonic speed regime.
The DLM is based on the potential flow theory and the flat surface representation which suffer from the inability to capture phenomena such as:

•
Airfoil camber and thickness as opposed to the standard flat plate results obtained from DLM, • Compressibility effects including local re-compression shocks, • Strongly non-linear aerodynamic lift and drag forces and sectional moments resulting from viscous flow phenomena like separation.
In terms of modeling, the numerical model for aerodynamic loads in NASTRAN involves the generation of flat plate lifting surfaces that represent the actual wing with two sections.These surfaces are then discretized into boxes, each of which includes a control point (a point that lies on the 3/4 of the chordwise distance and at the mean of the spanwise distance) for a flow tangency reference and a point at which pressure will be calculated.
The discretization of those surfaces is utilized with 30 elements in the chordwise direction and a total of 120 elements in the spanwise direction.The model is constructed based on the Outer Mold Line (OML) of the wing.and thus, remains constant during the optimization process.The model is presented in Figure 3.In order to achieve a realistic representation of the wing and take into account the camber and twist effects, the W2GJ matrix that models the local Angle of Attack (AoA) of each panel must be populated.The local AoA consists of two components, the local twist and the inclination due to camber.The component of the twist is calculated from the provided twist distribution but the camber is unknown for the specific airfoil.The method used to calculate the mean camber line is the construction of circles that are tangent to the upper and the lower curve of the airfoil (Figure 4).The center of each circle constitutes a point of the mean camber line and the inclination of it is calculated using numerical derivation.After the calculation of the twist and inclination due to camber distributions, the products are summed and the value of the total inclination is interpolated given the coordinates of each panel's control point.A section of the wing indicating the camber and twist values is illustrated in Figure 5.For the static aerodynamic loading, the resulting loads and moments from the corrected model have been verified against high fidelity RANS CFD analyses and are presented at a greater extent in [26].The accuracy of the flutter solution computed from the DLM model has also been validated and verified against the AGARD wing case study [27], with the results of the analysis being presented in [26].

Fluid-Structure Interaction (FSI) Coupling
Since the aerodynamic forces act on the boxes of the DLM model and the forces for the static analysis should be applied to selected nodal points of the FEM model, a transformation between the two computational models is necessary.This process is called splining and several spline methods exist for the interpolation between those two models ( [28]).In terms of modeling, an infinite plate spline is used herein for the aerodynamic forces mapping and the structural displacements transfer to the DLM model.However, the topology of spline points is critical to achieve an accurate coupling.Several topologies were studied and the most accurate, based on the results were the ones with points located near the chord of each rib.This topology had the smoothest distribution of the loads because the loads are transferred through the ribs in the rest of the structure.This spline topology is illustrated in Figure 6.

Aeroelastic Scaling Framework
Starting off, we will first refer to the routines developed for the aeroelastic analysis which form the core of the optimization and scaling algorithms.The developed framework combines in-house MATLAB codes, generating and manipulating MSC Patran and NASTRAN executable files and solvers.

Static Analysis
The static aeroelastic analysis framework starts with the initialization of the design vector.This vector contains all the information about the internal geometry, thickness distribution, and material properties of each component.The design vector consists of the following parameters: Number of stringers 3.
Rear spar position • Thicknesses

1.
Six control points per component for the thickness distribution 2.
Thickness of spar cap (L-type cross section) • Material

Density of the material 2.
Poisson ratio 3.
Young's modulus After the design vector is defined, the wing's design parameters are inserted and the session file containing the PCL commands in MSC Patran is generated.Then, a NASTRAN input file that contains information about the structural model is exported.Based on the wing's design parameters, the lifting surfaces are also defined and discretized.Then, a process that calculates the local AoA takes place.This process is based on the methodology presented in Section 2.2.The W2GJ matrix is then constructed, and together with the flight conditions, is incorporated into a file that automatically generates the final executable file through MSC Patran.The combined file that contains both the structural and aerodynamic models is compiled with SOL 144 of NASTRAN.The results are then extracted and postprocessed in order to calculate the objective function and the constraints.A detailed flowchart is presented in Figure 7.

Dynamic Analysis
Continuing with the dynamic analysis, at first a modal analysis routine is conducted.More specifically, the case control section of the input file is modified and compiled with the SOL 103 of NASTRAN.The results are extracted and post-processed in order to be used for the calculation of modal constraint objectives.Constraints and objectives vary depending on the type of application.In fact, the optimization of the wing's mass does not have any explicitly defined objective related to modal response.However, in similarity optimization, one of the main objectives is the error between the two models (the definition of this error will be given in the following chapters).A detailed flowchart that summarizes the aforementioned is presented in Figure 8. Regarding the dynamic aeroelastic analysis routine, the input file of the static aeroelastic analysis constitutes the core of it.Particularly, the routine starts with the extraction of the minimum eigenfrequency in order to begin the construction of MK (Mach and reduced frequencies) NASTRAN cards.After the construction of the MK cards, the density ratio is defined externally, and both are written in the input file.The case-control section parameters are modified and the bdf file is compiled with the NASTRAN's SOL 145.The results for each MK pair are extracted and post-processed to compute the flutter speed.A detailed flowchart is provided in Figure 9.As for the optimization routine, it is executed under the serial call of the static aeroelastic, modal and dynamic aeroelastic analysis.After each of the analyses is completed, the values of the objective functions and constraints of the problem are extracted and forwarded to the optimization algorithm which is an implementation of the MIDACO [29] ACO method.A detailed flowchart is illustrated in Figure 10.

Scaling
Aeroelastic scaling, although also an optimization process, consists of a slightly different methodology without any change in the fundamental concepts.More specifically, as already mentioned, scaling is a process of optimizing the structural and flow similarity between a reference wing and its corresponding sub-scale.The reference wing is initially optimized with the structural and dynamic responses being extracted so that they can be later used in the calculation of the objective function and the constraints.Given that, the scaling factors, the objective function and the constraints are formulated.The optimization problem is also handled via the MIDACO solver.A flowchart of the proposed methodology is presented in Figure 11.The overall proposed framework is also publicly accessible via the following link: https://github.com/VagelisFilippou/Aeroelastic-Scaling(accessed on 29 December 2023).

Reference Wing Geometry
In our study, the undeflected Common Research Model (uCRM) wing constitutes the reference wing geometry.Regarding the internal configuration of the uCRM wing, a two-spar configuration is considered.Spar caps and stringers of L and Hat cross-sections, respectively, are also present.A view of the OML, as well as of the internal structure of the uCRM wing, can be seen in Figure 12.Concerning the material, all relevant components are assumed to be manufactured of aircraft-grade aluminum, whose properties are summarized in Table 2. Detailed specifications of the uCRM wing are presented in Table 3.The wing twist distribution is provided by means of CAD files in [30] and is extracted via a MATLAB script (more information about the uCRM wing can befound in: https:// commonresearchmodel.larc.nasa.gov/(accessed on 29 December 2023)).

Critical Aerodynamic Loads
The critical aerodynamic loads considered correspond to a +2.5g pull-up maneuver at sea-level conditions.Prior studies [30,31] have identified the +2.5g condition at a Mach number of 0.64 as the critical load case for the reference wing, as summarized in Table 4, with the MTOW of the aircraft being equal to 297,500 kg.In this study, the single objective is the minimization of the structural mass of the wing submitted to a series of stiffness, strength and aeroelastic instability constraints alongside a fairly wide design space that is based on the freedom of the framework developed.The objective and constraints of the problem are summarized in Table 5.
The addition of constraints completes the structural optimization problem, providing more meaningful and requirements-compliant designs.For our purposes, static stiffness, strength and aeroelastic stability constraints are considered for the optimization framework.Any form of buckling constraint critical to the structural sizing procedure has not been included in the formulation of the optimization problem.Although not clearly stated as design requirements, static stiffness constraints are accounted for in the majority of optimization studies of aircraft wings present.For our purposes, the maximum deflection at the tip of the wing corresponding to the aforementioned loading scenario is constrained as in [32]: where s is the wing semi-span.Additionally, a maximum twist angle of 6 degrees at the tip of the wing is also imposed in order to prevent undesirable aeroelastic phenomena [33]: 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 as follows: For the evaluation of the static strength of each candidate design, we used the Von Mises equivalent stress.A safety factor equal to 1.5 was used for all the components with a yield strength of 420 MPa (Table 2): The aeroelastic stability by means of flutter speed is also investigated for each candidate design.The PK method is employed by NASTRAN for the calculation of the stability parameters.The flutter analysis is performed at the relevant flight conditions described in Table 4. Regarding the constraint, the flutter speed shall exceed the dive speed with a safety factor of 1.2:

Variables Definition
As far as the design variables are concerned, they are those stated at Section 2.5.Because of manufacturing constraints and for the avoidance of oversized design evaluations, the range of thickness values is constrained with upper and lower bounds.Also, the location of the front and rear spar is constrained in order to avoid bugs due to spar inversion.The number of stringers and ribs is also constrained in order to avoid bugs caused by geometry definition.The DV and their bounds are summarized in Table 6.The optimization procedure, as we already stated, is handled via the MIDACO solver, which adopts a combination of an extended ACO algorithm with the Oracle Penalty Method, an advanced method developed for metaheuristic search algorithms, for constraint handling of the solution process.In order to achieve a more efficient and accurate result, a cascading approach with multiple runs has been adopted in combination with MIDACO's FOCUS parameter.FOCUS is an internal parameter of the algorithm that forces the MIDACO solver to focus on the current best solution found.This parameter was initially set to zero in order to reach the global optimum solution.Then, the best solution obtained is used as a starting point for the next optimization run.In this run, the algorithm focuses more on the region around the initial point by increasing the FOCUS parameter.This procedure is repeated until user-specified stopping criteria like maximum evaluations or time is reached.The parameters are presented in Table 7.By running the previously described optimization problem setup, the optimizer converged in an optimal solution after two runs (as defined in Table 7) and 350 iterations in total.The optimization results in terms of the objective function and constraints are depicted in Table 8.The resulting internal geometry is presented in Figure 13.Table 9 summarizes the results of geometric parameters and in Figure 14 the optimum thicknesses are illustrated.At first linear static analysis was performed in the condition of full tanks.The AoA was adjusted in a way to produce a +2.5g load as pinpointed earlier.The maximum deflection occurs at the tip of the wing as shown in Figure 15 with a value of 2.66 m ≤ 0.15• Semi-Span (constraint according to [32]).Deflection gradually arises in the spanwise direction as was expected by the gradual decrease in thickness.
An important result of the analysis is the deflection in the longitudinal direction (x direction), as illustrated inFigure 16.The maximum deflection in this direction is around 0.11 m, a value that is relatively small compared to the vertical deflection's maximum value.That fact indicates the inability of the DLM method to accurately capture the drag distribution due to the assumption of inviscid flow.Another product of these considerations is the convergence of the position of the spars in the region of maximum thickness of the airfoil, a fact which gives a very high stiffness in the x direction but reduces the stiffness in the z direction quite a bit.This comes from underestimating the drag distribution.Another important result of the analysis was the spanwise twist of the wing that is presented in Figure 16.The twist arises gradually from the root of the wing and increases with a constant gradient until it reaches an almost constant value towards the tip of the wing.This phenomenon arises from the sweep of the wing and the relatively aft positioning of the spars that also move the elastic axis aft.According to [33], the maximum wing torsion must not exceed the value of 6 0 .In this analysis, the Tip Torsion = 3.85 0 ≤ 6 0 .As far as the stresses are concerned, the maximum occurs at the Yehudi break as expected with the material operating at its linear region.The rear spar is highly stressed in the root of the wing and the stresses gradually decrease spanwise (except the Yehudi break).The other components follow the same pattern except for the fact that they do not stress as much at the root of the wing.The stress distribution of all the components is illustrated in Figure 17.
ersion February 20, 2024 submitted to Aerospace 16 of 3 At first linear static analysis was performed in the condition of full tanks.The AoA wa adjusted in a way to produce a 2.5G load as pinpointed earlier.The maximum deflection occurs at the tip of the wing as shown in Figure 15 with a value of 2.66m ≤ 0.15• Semi-Span (constraint according to Dababneh et al.).Deflection gradually arises in the spanwis direction as was expected by the thickness gradual decrease.An important result of the analysis is the deflection in the longitudinal direction (X direction), as it is seen in Figure 16.The maximum deflection in this direction is around 0.11m, a value that is relatively small compared to the vertical deflection's maximum value.That fact indicates the inability of the DLM method to accurately capture the drag distribution due to the assumption of inviscid flow.Another product of these consideration is the convergence of the position of the spars in the region of maximum thickness of th  the material operating at its linear region.The rear spar is highly stressed in the root of the 403 wing and the stresses gradually decrease spanwise (except the yehudi break).The other 404 components follow the same pattern except of the fact that they do not stress as much at 405 the root of the wing.The stress distribution of all the components is illustrated in Figure 17.406 As far as the aerodynamics are concerned, the pressure seems to have a reasonable 407 distribution based on the twist and camber of the wing, with lift almost pertaining an ellip-408 tical distribution.A further reduction in lift in the spanwise direction takes place because 409 of the downward motion of the leading edge.The coefficient of pressure distribution is 410 illustrated in the Figure 18.As far as the aerodynamics are concerned, the pressure seems to have a reasonable distribution based on the twist and camber of the wing, with lift almost pertaining to an elliptical distribution.A further reduction in lift in the spanwise direction takes place because of the downward motion of the leading edge.The coefficient of the pressure distribution is illustrated in Figure 18.

Dynamic Aeroelastic Analysis
The wing flutter analysis was performed with full fuel load in this section.At first, a modal analysis was performed because the state of the structure is expressed as a function of its natural modes.Also, the modal analysis estimates the first eigenfrequency, and thus, the lowest reduced frequency in order to use it for the flutter analysis as a starting point.The required inputs of density and velocity correspond to the flight conditions, and the reduced frequencies studied, starting from the reduced frequency obtained from the first eigenfrequency.The results are depicted in Figure 19.Since the flutter is defined as the point where the oscillations (naturally damped before the flutter speed) start to diverge, the flutter point is identified when the damping plotted against the airspeed (denoted as V) turns from negative to positive values.In Figure 19 and through linear interpolation we find the first mode turns from negative to positive damping values at around 530 m/s.

Scaling Parameters Definition
For the scaled model, we consider a reduced span at the level of 6m in order to meet typical wind tunnel dimensional constraints (ONERA Modane wind tunnel [34]).Let us start by defining the length ratio, λ l .As we know the full-scale airplane span and the desired span for the scaled model, the relation is straightforward: with the full-scale characteristics represented by the subscript w and the scaled models by the subscript m.The size and aerodynamic shape of the scaled model are completely defined by this length ratio.The next step will be the matching of the Froude number, so that: this relationship will give the airspeed ratio, λ V , between scaled and full-scale models, the scaled airspeed is V m = √ λ l V r = 98.83 m/s.At the scaled model altitude, and using the international standard atmosphere model, this gives a Mach number of M = 0.29.Now, having both the length ratio and the velocity ratio set, one can extract the time scale, in this case, given by the frequency ratio, λ ω .Taking into account the reduced frequency, κ = bω V , we have: From the mass ratio, µ = ρSb m , we obtain: and knowing that the baseline chord of both full-scale and scaled model is related by c m = h m b w c w and that the wing area is S = bc, one can relate the mass ratio, λ m , and air density ratio, λ ρ , as: For the present work, a λ ρ = 1 was considered; therefore, λ m = (125) −1 = 8 × 10 −3 .As noted in [35], if cruise conditions for the full-scale and scaled models were used for the calculation of the air density ratio, the resulting mass ratio would be unrealistically low and impractical for structural optimization.

Objectives and Constraints
This study aims to optimize both the static aeroelastic and modal similarity of the scaled model compared to the reference.The reference wing is the already optimized wing presented in Section 3.1, and it defines the scaled model's target responses modified by the parameters defined in Section 3.2.1.The framework utilizes a direct optimization of both the modal and static aeroelastic responses subjected to stiffness and mass distribution variables thus, constituting it as a multi-objective optimization.In the next section, a brief review of the objectives and constraints formulation and methodology is presented.

Static Aeroelastic Similarity
The problem of static aeroelastic disparity between the two wings arises from the inability to generate the same flow conditions and the difference in stiffness distribution.In order to eliminate any difference between the in-flight shape, we have to define an objective function that describes it.In that case, we will need the in-flight deflections of the reference aircraft (for a particular AoA), and to evaluate the fitness of a particular design we need to perform a complete static aeroelastic analysis for the given design variables and the flow parameters (i.e., the unmatched Mach number and the given AoA of the reference model for the critical loading condition) and then compare the in-flight shape to the scaled version of the reference one, as well as evaluating the acceleration and stress constraints.With this optimization problem, our goal is to find the wing design parameters that minimize the error in the in-flight shape with respect to the scaled in-flight shape of the reference aircraft.
As described in [35], the scaling factor and the air density ratio λ ρ = ρ m /ρ r determine the total mass of the scaled model: which will dictate the amount of lift that needs to be produced, and thus, the desired AoA.In order to maintain the +2.5g acceleration the AoA was adjusted and set at the value of 0.8 0 .Also, and since we want to ensure the structural feasibility of the most critical in-flight situation, we require the Von Mises stress in any point of the wingbox structure to be lower than the allowable stress of the structural material for the +2.5g acceleration condition.Thus, the proposed formulation of the optimization problem is the following: where {X c a } are the shape coordinates of the scaled model wing, X re f a are their counterpart of the reference aircraft, and λ l is the overall geometrical scale factor (i.e., wing span ratio).

Modal Similarity
We will now present and discuss the choices adopted for the modal matching criteria (i.c., what quantities we set as the objective function and constraints).Concerning the formulation of the optimization problem, we choose to use the MAC metric to represent the correlation between two vibration modes.As described in [36], if {ϕ ri } is one of the reference eigenvectors, and ϕ mj is one of the eigenvectors from the current model, we define a matrix whose elements are: which is the normalized dot product between two vectors {ϕ ri } and ϕ mj representing two modes.Indeed, a MAC value of 1 indicates that the two vectors represent the same mode, whereas a MAC value of 0 indicates that the two modes are orthogonal.Given a set of N reference modes and N modes that we want to evaluate, we define this N × N matrix containing the MAC value for each possible pair of modes between the reference ones and the ones being analyzed.Given the matrix [MAC] we define the objective function to be minimized as: where tr() indicates the trace of a matrix, [Φ r ] is the matrix containing the reference modes and [Φ m ] contains the modes being analyzed.Concerning the scaled frequency matching, we set an equality constraint for each scaled frequency as: where λ ω is the frequency ratio, as demonstrated in [13] Also, we set an equality constraint for the scaled mass where λ m is the mass ratio.By doing so, we aim to avoid having a multiplicity of solutions that closely match the reference modal shapes and frequencies.By using the mass constraint we choose the solution that matches the scaled mass, in addition to mode shapes and frequencies.
A problem encountered during the programming of the objective functions was the change in the coordinates of the mesh nodes compared to the reference wing.In fact, because of the re-meshing at each design point, the grid is different at each evaluation, resulting in the impossibility of directly calculating the error at all nodes.This problem was addressed by computing the displacements at nodes that are closer to the reference nodes defined by the full-scale model and then used directly either to compute the in-flight Root Mean Square (RMS) error or the [MAC] matrix.Overall, the scaling of the uCRM wing will be achieved by a single-step optimization algorithm that utilizes two objective functions, one for the in-flight shape difference and one for the modal similarity.Frequency, mass and stress constraints are also taken into account in order to achieve feasible solutions.The objectives and constraints of the problem are summarized in Table 10 and Table 11, respectively.The target frequency values are also presented in Table 12.

Variables Definition
As far as the design variables are concerned, they are stated at Section 2.5.Because of manufacturing constraints and for the avoidance of oversized design evaluations, the range of thickness values is constrained with upper and lower bounds.Beyond the geometrical and thickness variables, the addition of 10 lumped masses in order to control the dynamic response was deemed essential.Those 10 masses are equally distributed and adjusted according to their upper and lower bounds (Figure 20).The DV and their bounds are summarized in Table 13.The deformation of the scale wing was relatively low compared to the reference wing, and 532 the thickness variables converged towards their lower limit.So a change was made in the 533 material, and magnesium was chosen, whose properties are shown in the Table 14.

Quantity Value Units
Density 1800 kg/m 3  At first, a single-step optimization run was performed with the modal objective accounting for the 10 first modes, and an accuracy level of 0.5 was set for the constraints handling.The modal objective could not provide an acceptable level of matching, but the constraint violation converged to an allowable level.In order to have acceptable matching, the modal objective should be lower than 0.2-0.3.Also, the feasible points found did not exhibit an in-flight shape that matched the target shape with acceptable accuracy.The leading cause of the inability to converge to an acceptable value was the high elastic modulus of aluminum, in combination with the reduced aerodynamic loads.The deformation of the scale wing was relatively low compared to the reference wing, and the thickness variables converged toward their lower limit.A change was made in the material, and magnesium was chosen; the properties are shown in the Table 14.In order to reduce the computing time and converge to feasible solutions with acceptable accuracy, the optimization routine was modified to a serial optimization where the flight shape difference is first optimized.Then, the modal response of the wing is also optimized while the optimized stiffness parameters are kept constant.A downside of this choice is that the resulting optimized parameters are potentially sub-optimal to the scaling problem, as stated in [23].

Static Aeroelastic Response
Static aeroelastic similarity optimization was set and run twice with different FOCUS parameter values (Table 15).The tolerance for the violation of the Von Mises constraint was set to 2. In less than 500 iterations, the RMS error of the in-flight shape difference converged to an optimal value of 0.0587 m.After the first run, the focus parameter was set to 10, and the routine was recompiled and converged to an optimal value of 0.0549 m or RelativeError ∼ 10%.The optimal thickness DVs are depicted in Figure 21 and the geometric DVs are presented in Table 16.In Figure 22, we see the in-flight shape of the best-found point (with the specified tolerance) compared to the in-flight target shape scaled from the reference aircraft.The in-flight shape fits well overall along the wing, except for the nodes near the wing tip.In fact, as Figure 23 represents, an inaccuracy is observed at the AoA of the sections near the wing-tip.A solution to this problem would be the consideration of wing tip torsion equality between the scaled and the reference wing as a constraint.
Regarding the stress constraint, it is not active for the best-found point.As we can see in Figure 24, the maximum von Mises stress on the structure for the limit load condition is 149 MPa, compared to the defined allowable of 150 MPa.

Modal Response
The previously described mode similarity optimization problem was performed for the first N = 5 reference modes using the MIDACO optimizer.The frequency constraints converge to the required value of 0, with a maximum relative error of 0.63 % for the scaled natural frequencies.The scaled total mass constraint converges with an error of 0.19 %.The minimum of the objective function characterizing the mode shape similarity that satisfies the constraints within the mentioned errors is f min = 0.066, which implies that the average MAC value of the N = 5 modes considered for the comparison is 1 − f min = 0.934.Therefore, the results indicate a good matching for the first five modes according to [36]'s guidelines.However, mode shapes 6-10 are dissimilar, and the improvement of the framework is necessary to fully match the modal response.The mass DVs are also presented in Figure 25 and their corresponding values are those presented in Table 17.Finally, The eigenfrequency of each mode is presented in Table 18.The average error of the eigenfrequencies is 14.44%.An attempt for optimization of N=10 modes was also conducted.Nevertheless, the 5 optimizer did not converge to a satisfying objective value, rather focused on the demanding 5 constraint fulfilment.A multi-start optimization technique would give us a more correct 5 view of the convergence ability of the problem.As a result, it would be better to focus on 5  Figure 26 shows the comparison of the reference modes and the ones of the best design found by MIDACO after 170 iterations.
An attempt for the optimization of N = 10 modes was also conducted.Nevertheless, the optimizer did not converge to a satisfying objective value but rather focused on the demanding constraint fulfillment.A multi-start optimization technique would give us a more correct view of the convergence ability of the problem.As a result, it would be better to focus on the exact matching of a more limited number of modes, instead of trying to match relatively large numbers of modes (e.g., N = 10) and obtain a lower average result.

Flutter Speed
The results are depicted in Figure 27.Since flutter is defined as the point where the oscillations (naturally damped before the flutter speed) start to diverge, the flutter point is identified when the damping plotted against the airspeed turns from negative to positive values.In Figure 27 we find that the first mode turns from negative to positive at around 315 m/s.

Effectiveness of Geometrical Parameterization
The main novelty of this study is the introduction of geometric parameters in the optimization framework.The number of ribs is a parameter that has a significant role in stiffness and mass distribution.It affects the spanwise torsional stiffness of the wing, and thus, is a parameter that controls the wing's torsion, which has a significant role in static and dynamic aeroelastic phenomena.Moreover, the spars' position is an essential tool for the control of the bending stiffness of the wing.Adjusting the spars' position affects the moment of inertia in the x and z direction (because of the airfoil's thickness distribution), and consequently, the displacements in the z and x directions, thus, achieving a better matching between the two wings.Consequently, by adjusting those three variables, we can effectively affect the torsional stiffness, bending stiffness and mass distribution.
In Figures 28 and 29, the convergence of those variables is demonstrated.It is obvious that the algorithm tried to adjust the parameters appropriately so that there is a substantial reduction in torsional and bending stiffness in order to achieve greater similarity.The aerodynamic loads of the sub-scale model were noticeably lower, so reducing the number of ribs and placing the spars towards the trailing and leading edges reduced stiffness considerably.These parameters played a decisive role while the thickness parameters adjusted the spanwise stiffness distribution to achieve a similar deformation pattern.To fully comprehend the efficiency of the proposed aeroelastic scaling framework, the results of the optimized aeroelastic scaled model have been compared with results considering fixed geometrical internal parameters in terms of the MAC value.The results are summarized in the following Table 19: Overall, the above results indicate the superiority of the proposed framework with respect to a classic aeroelastic scaling approach considering fixed geometrical parameters.

Conclusions
This research study consists of the development of an aeroelastic optimization framework, which with appropriate configuration is useful for either mass or similarity optimization between a sub-scale model and a reference wing.The core of the framework consists of the aeroelastic analysis (static and dynamic) combined with the appropriate functions for the objectives and constraints, and an ACO algorithm (MIDACO).The parameterization of the wing includes geometric parameters and thicknesses for the stiffness distribution and 10 lumped masses for the mass distribution (for the scaling problem).
A mass optimization problem was initially set up for the uCRM wing.In this case, there were only geometric and thickness parameters (e.g., spar positions, number of ribs, number of stringers and thicknesses at specified control points).Stiffness, stress, and flutter constraints were included with target values being set based on the literature and the Civil Aviation Regulations FAR-25.The optimizer converged at a structural mass of 6902 kg while all the constraints were satisfied.Of interest are the optimized variables' values.The thickness parameters decrease gradually along the wing, except for the Yehudi break where maximum stress usually arises.The spars' position, as expected, converged around the position of the airfoil's maximum thickness in order to optimize the stiffness in the x-direction while the underestimated drag (because of the potential flow implemented by DLM) does not require so much stiffness in the z-direction.Finally, the ribs' number converged to 51, which is close to the optimal solution presented in [30].
As far as the scaling problem is concerned, the objectives and constraints were modified in order to meet similarity requirements.In fact, the problem's objectives were the RMS error of in-flight shape difference and the average MAC of the N modes that were selected to optimize.The constraints included were the equality of each eigenfrequency, the equality of the model's mass to a target value defined by the scaling parameters, and the avoidance of maximum stresses greater than the maximum allowable.Also, a constraint to avoid flutter at a certain speed for the optimal solution was considered.
The pre-described optimization problem was first run with a single-step approach, for N = 10 modes to optimize and the material was aluminum.The results were deemed unacceptable.The RMS error of the in-flight shape difference converged to the value of 0.45, indicating bad matching of the static aeroelastic response.The constraints of the modal response were satisfied but the modal similarity was quite low indicating no similarity at all.These problems were addressed by changing the material to magnesium, changing the modes to match N = 5 and the one-step approach to a two-step (serial) approach in order to reduce the computing time and obtain acceptable results quickly.After those alterations, the optimizer converged to acceptable values indicating a good matching of both the static aeroelastic and modal response.In fact, the RMS error of the in-flight shape difference of the optimal solution is 0.0549 m and the mean MAC value of the first five modes is 1 − f min = 0.934.The constraints of the problem were also satisfied but there is room for improvement mainly in terms of frequency error.As far as the flutter is concerned, the flutter speed of the optimal solution is beyond the diving speed, indicating no flutter.
Of interest is the location of the spars, which are quite close to the leading and trailing edges, in contrast to the reference wing.This is likely due to the increased stiffness presented by the sub-scale structure given the reduced aerodynamic loads.Placing the spars towards the edges reduces the stiffness around the x-axis, and thus, the structure becomes more flexible, avoiding the very low values of the thicknesses and the construction difficulties.The number of ribs was also reduced to 15 in order to reduce the torsional stiffness of the wing indicating the significance of the parameterization of ribs' number.Future implementations in the framework are directed towards the following:

•
The investigation of other types of materials (such as wood or composites) and the combination of them would greatly increase the design space and ensure the presence (or absence) of the effects of the full model (such as skin buckling).

•
The calculation of derivatives in order to include a gradient-based optimizer in the framework.It will give the opportunity for a more efficient optimization since we take advantage of a gradient-free optimizer.

•
The introduction of a high-fidelity aerodynamic solver such as Euler or RANS.As already mentioned, the potential theory upon which the DLM is based underestimates the drag distribution while overestimating the lift.A high-fidelity solver will avoid such problems and provide more accurate loading and flutter prediction.

•
The parameterization of the number of spars will give the opportunity to design wings with more than two spars, and it will greatly increase the design space.Also, it will allow the design of other types of wings, such as military aircraft wings, that incorporate multi-spar configurations.

•
The verification and validation of the efficiency and accuracy of the proposed aeroelastic scaling framework against already available frameworks.

Figure 2 .
Figure 2. Example of FEM mesh with shell and beam elements.

Figure 6 .
Figure 6.Example of coupling nodes in a wing.

Figure 13 .
Figure 13.3D view of the optimized geometry.

Figure 14 .
Figure 14.Optimized thickness distribution of the internal geometry (top), upper (middle) and lower (bottom) skins along with the corresponding control point values.

Figure 16 .
Figure 16.Wing static aeroelastic displacements in the x (top left), z (top right) direction and torsion angle (bottom).

Figure 17 .
Figure 17.Von Mises stress distribution of the uCRM wing 411

Figure 21 .Figure 22 .
Figure 21.Optimized thickness distribution of each component and its corresponding control point values for the scaling problem.

Figure 22 .
Figure 22.In-flight aerodynamic surface of the optimized design compared to the target shape.

Figure 23 .
Figure 23.Comparison of the scaled and reference wing in terms of deflection (left) and torsion angle (right).

Figure 24 .
Figure 24.Von Mises stress distribution of the scaled model.

Figure 25 .
Figure 25.Optimized lumped masses of the scaled model Figure 26 shows the comparison of the reference modes and the ones of the best design 5 found by MIDACO after 170 iterations.

Figure 25 .
Figure 25.Optimized lumped masses of the scaled model.

Figure 26 .
Figure 26.First eight eigenmodes of the scaled model compared to the modes of the reference wing (red indicates undeformed shape and blue indicates the deformed shape).

Figure 27 .
Figure 27.V−g (left) and V−f (right) plots of the optimized scaled wing.

Figure 28 .
Figure 28.Spars' position during the 1st run of static aeroelastic response optimization.

Figure 29 .
Figure 29.Number of ribs during the 1st run of static aeroelastic response optimization.

Table 1 .
Summary of the features of the publications referenced herein (✓ − Present, ✗ − Absent).

Table 4 .
CRM wing critical loading conditions.

Table 5 .
Objectives and constraints of the mass minimization problem.

Table 6 .
Design variables of the optimization problem.

Table 10 .
Objectives and constraints of the scaling problem.

Table 11 .
Objectives and constraints of the mass minimization problem.

Table 12 .
Target values of frequency for the optimization problem.

Table 13 .
Design variables of the optimization problem high elastic modulus of aluminium in combination with the reduced aerodynamic loads.531

Table 13 .
Design variables of the optimization problem.

Table 15 .
MIDACO parameters of the static aeroelastic response similarity optimization.

Table 17 .
Optimized mass values for the modal response similarity optimization problem.

Table 18 .
Target values of frequency for the optimization problem.

Table 19 .
Effectiveness of the proposed aeroelastic scaling framework.