Multidisciplinary Design Optimization of the NASA Metallic and Composite Common Research Model Wingbox: Addressing Static Strength, Stiffness, Aeroelastic, and Manufacturing Constraints

Abstract

This study explores the multidisciplinary design optimization (MDO) of the NASA Common Research Model (CRM) wingbox, utilizing both metallic and composite materials while addressing various constraints, including static strength, stiffness, aeroelasticity, and manufacturing considerations. The primary load-bearing wing structure is designed with high structural fidelity, resulting in a higher number of structural elements representing the wingbox model. This increased complexity expands the design space due to a greater number of design variables, thereby enhancing the potential for identifying optimal design alternatives and improving mass estimation accuracy. Finite element analysis (FEA) combined with gradient-based design optimization techniques was employed to assess the mass of the metallic and composite wingbox configurations. The results demonstrate that the incorporation of composite materials into the CRM wingbox design achieves a structural mass reduction of approximately 17.4% compared to the metallic wingbox when flutter constraints are considered and a 23.4% reduction when flutter constraints are excluded. When considering flutter constraints, the composite wingbox exhibits a 5.6% reduction in structural mass and a 5.3% decrease in critical flutter speed. Despite the reduction in flutter speed, the design remains free from flutter instabilities within the operational flight envelope. Flutter analysis, conducted using the p-k method, confirmed that both the optimized metallic and composite wingboxes are free from flutter instabilities, with flutter speeds exceeding the critical threshold of 256 m/s. Additionally, free vibration and aeroelastic stability analyses reveal that the composite wingbox demonstrates higher natural frequencies compared to the metallic version, indicating that composite materials enhance dynamic response and reduce susceptibility to aeroelastic phenomena. Fuel mass was also found to significantly influence both natural frequencies and flutter characteristics, with the presence of fuel leading to a reduction in structural frequencies associated with wing bending.

1. Introduction

Multidisciplinary design optimization (MDO) has emerged as a prominent methodology in aircraft design and mass estimation, enabling more efficient and integrated approaches to wing structure design. This approach allows for the concurrent application of diverse numerical optimization techniques across all relevant design disciplines. A comprehensive review of MDO methodologies and their applications is presented by Martins et al. [1], highlighting the advancements made in this area. Bindolina et al. [2] proposed a multilevel structural and multidisciplinary optimization procedure for the preliminary estimation of aircraft wingbox mass, where conventional empirical formulas and statistical methods may not provide sufficient reliability. Similarly, Hurlimann et al. [3] introduced a CAD/CAE-based multidisciplinary process for the mass estimation of transport aircraft wingbox structures, demonstrating the success of their approach in yielding satisfactory results. More recently, Gray et al. [4] developed OpenMDAO, an open-access platform for multidisciplinary design, analysis, and optimization. This platform incorporates advanced algorithms capable of addressing optimization problems with significant complexity and scale, thereby facilitating the resolution of intricate multidisciplinary engineering challenges. For large-scale, high-fidelity structural optimization problems involving a substantial number of design variables, gradient-based optimization algorithms prove to be particularly beneficial, as they reduce the number of function evaluations needed to reach a local optimum [5,6,7,8].

Weight optimization plays a crucial role in enhancing both the structural integrity and aeroelastic performance of an aircraft. By reducing wing mass, the stress on the overall structure is minimized, leading to improved fuel efficiency and overall performance. Achieving an optimal integration of structural, aeroelastic, and material considerations is vital, especially in the design of next-generation aircraft wings. The NASA Common Research Model (CRM), developed as a generic transport aircraft model, serves as a notable example [9,10]. The CRM was specifically designed to represent a long-range, twin-aisle transport aircraft configuration. Its geometry includes a fuselage, wing, horizontal tail, engine nacelle, and pylon, and it was originally intended to facilitate the verification and validation of computational fluid dynamics methodologies [11,12,13,14,15]. Since its introduction, the CRM has proven to be a valuable aerodynamic benchmark for predicting drag and optimizing aerodynamic design based on computational fluid dynamics. In recent years, there has been growing interest in expanding the CRM to include other types of studies. Various efforts have been made to develop a CRM wing structural model that accurately represents the structure, intended for use in aerostructural and aeroelastic design, analysis, and optimization studies [16,17]. Jutte et al. [10] developed an all-metallic, single-material baseline wingbox structure for the CRM consisting of wing skins, an orthogonal grid of primary ribs and spars, as well as secondary stringers and rib-stiffened structures. Their study investigates the use of tow-steered composite laminates, functionally graded metals (FGM), thickness distributions, and curvilinear rib/spar/stringer topologies for aeroelastic tailoring. One key recommendation from their work is the application of formal optimization, with gradient-based optimization being identified as the most cost-effective method, especially for systems with numerous design variables. This approach should be employed to gain a deeper understanding of the design space while ensuring careful design to prevent structural failure. In addition to static aeroelastic behavior, it is crucial to incorporate other load cases and design constraints into the optimization process. These additional load cases, which account for factors like engine mass, landing gear, leading and trailing edge masses, and fuel tanks, will impact the optimized weight, ultimately leading to an increase in weight as the structures must be capable of supporting all load cases.

A structural optimization of an aircraft wing was performed by [18] during the early design phase, aiming to reduce wing bending moments, minimize structural mass, and enhance fatigue life, subject to constraints on material strength, buckling stability, and static aeroelasticity. The implementation of load alleviation techniques led to a reduction in wingbox mass of 2.8% for the backward-swept wing configuration and 6.1% for the forward-swept configuration. With advancements in computational techniques, optimization methods have become extensively applied in aeroelastic optimization. Previous studies in aeroelastic optimization have primarily concentrated on the design of wing structural stiffness, as well as the development of multiscale aeroelastic optimization methodologies that simultaneously address structural and material design under representative flight loading conditions [19,20,21]. In their work, Kafkas et al. [22] proposed a multi-fidelity optimization framework for state-of-the-art composite aircraft wings, employing the Mixed Integer Distributed Ant Colony Optimization (MIDACO) algorithm. The primary case study focused on a modified CRM wing with planar geometry. Structural modeling was performed using either the Equivalent Plate Method (EPM) or a high-fidelity three-dimensional shell and rod representation, while aerodynamic loads were evaluated using a Vortex Lattice Method (VLM). Recommendations for future work included the incorporation of buckling constraints within the high-fidelity structural model and the integration of aeroelastic tailoring capabilities to enhance the fidelity and applicability of the framework. In their study, Wang et al. [23] employ a global–local optimization framework incorporating static failure, aeroelastic, buckling, and manufacturing constraints to determine optimal structural configurations for straight-fiber and variable-angle-tow (VAT) composite wingbox architectures. Localized buckling analyses are carried out on discrete wing subsections using high-fidelity finite element models, with sectional running loads extracted from a comprehensive aeroelastic simulation of the full wing structure. The optimization results demonstrate mass reductions of 12.5% and 13.2% for constant-thickness VAT and variable-thickness continuous tow shearing (CTS) designs, respectively, relative to a reference quasi-isotropic straight-fiber configuration. Dillinger et al. [24] investigated the effects of stiffness optimization, also known as aeroelastic tailoring, on the mass of forward-swept wings. The reader is encouraged to consult the reference list in [24] for additional studies on similar investigations. In their study, Handojo et al. [25] conducted a comprehensive structural optimization of the composite lifting surfaces using defined load cases while maintaining constant secondary structural masses and the existing wing planform geometry. The implementation of a suite of advanced load alleviation methods and prospective technologies resulted in a 26.5% reduction in wingbox mass, equating to a 4.4% decrease in the aircraft’s operating empty mass. In their work, the authors suggested that to achieve further mass reduction in future aircraft designs, additional design parameters should be reconsidered. Specifically, the imposed constraints on minimum structural thicknesses could be re-evaluated to assess the potential for further optimization.

Brooks et al. [26] defined an undeflected CRM to serve as a benchmark for aeroelastic analysis and design, addressing the challenges posed by 1g built-in deflections, which complicate aerostructural analyses of the CRM. The undeflected model is generated using an inverse design procedure, which integrates the outer mold line geometry of the wing with its corresponding internal wingbox structure. In their investigation, the wingbox components are optimized numerically, with all components modeled using aluminum alloy. Notably, the stiffeners are not explicitly included in the FEM model of the wingbox; instead, their stiffness contributions are incorporated by homogenizing them into the stiffness of the wingbox panels. The authors suggest that future research should focus on the development of variable fidelity models and the refinement of models to more accurately capture dynamic aeroelastic phenomena. Chauhan and Martins [27] employed a low-fidelity aerostructural design optimization approach, utilizing a simplified wingbox model, to optimize a wing based on the CRM wing. They compared their results with those of Brooks et al. [28], who applied high-fidelity models. The optimized wing mass reported by Chauhan and Martins was approximately 10% lower than the value obtained by Brooks et al. [28]. As anticipated by the authors in [27], differences in model fidelity led to variations in structural thickness distributions, optimized twist, and thickness-to-chord ratios, resulting in significant discrepancies in the optimized CRM wing masses.

Liu et al. [29] applied the integrated global–local optimization approach to the subsonic NASA CRM wing, demonstrating the scalability of the methodology in conjunction with medium-fidelity finite element method analysis. The optimization process aims to minimize wing weight while satisfying stress and buckling constraints. Specifically, the number of ribs and the configuration of stiffened panels are optimized to enhance structural efficiency. Compared to the unstiffened CRM wing, the integration of stiffened panels within the global–local optimization framework yields a substantial weight reduction of approximately 40%. However, the complexity of the design space associated with shape variables results in significant computational costs. To address this challenge, a parallel computing framework has been developed, substantially reducing computational time. Kennedy and Martins [30] conducted an aerostructural design optimization study on both metallic and composite aircraft wings, examining the impact of composite materials on the trade-off between structural mass and aerodynamic drag. Their structural model incorporates strength and buckling constraints, along with detailed laminate parametrization to determine the optimal stacking sequence while enforcing manufacturing constraints, such as matrix cracking and minimum ply-content requirements. Their findings indicate that composite wing designs achieve a weight reduction of 34% to 40% compared to their metallic counterparts. Kenway et al. [31] performed a high-fidelity aerostructural redesign of a subsonic aircraft wing for transonic flight conditions, employing a CAD-free geometry parametrization method. Their study captured the multidisciplinary trade-off between the wing’s structural mass and the sweep angle. The requisite level of finite element detail necessary to accurately characterize the structural geometry of aircraft wings for predicting the mass of the primary load-bearing wing structure, with an acceptable degree of accuracy using structural optimization techniques, has been investigated by the authors in [32]. Their findings demonstrate that although low-fidelity models are advantageous in terms of reduced computational time and minimal input requirements, they generally result in diminished accuracy in the estimation of wing mass, particularly when applied to composite wing structures. Thus, the application of MDO addresses the physics often omitted in simpler models. Weight optimization must be integrated with comprehensive structural analysis to ensure sufficient strength while minimizing weight, thereby requiring meticulous design to avoid structural failure.

This paper explores the multidisciplinary design optimization of the NASA Common Research Model wingbox, integrating a comprehensive set of realistic, engineering-relevant constraints. These include static stress, stiffness, aeroelastic flutter speed, and manufacturability considerations for both metallic and composite wing designs. Such constraints are essential for translating the optimization results into practical, real-world applications. The optimization process is applied to a high-fidelity finite element model of the CRM wingbox, with an emphasis on achieving structurally feasible and physically realizable solutions. The approach is driven by practical design considerations, rather than purely academic performance metrics, and it utilizes an efficient gradient-based optimization strategy. The study specifically evaluates the impact of integrating composite materials into the CRM wingbox design, focusing on the optimized reduction of structural mass relative to the metallic wingbox, as well as the effects on flutter speed while considering the effects of fuel mass under active flutter constraints. The paper is structured as follows. Section 2 provides an overview of the structural design and finite element modeling of the CRM wingbox. Section 3 outlines the aeroelastic modeling and flutter analysis methodology employed in this study. Section 4 defines the CRM wingbox optimization problem and its formulation. Section 5 presents the results and discussion of the CRM wingbox optimization case studies, emphasizing comparative reliability. Finally, Section 6 concludes the paper and offers recommendations for future work.

2. Structural Design and Modeling of the CRM Wingbox

2.1. Technical Data and Specification of CRM Wing

The publicly available NASA CRM aircraft is utilized in this study. The CRM represents a modern, single-aisle, transport-class aircraft configuration, developed as an open geometry for collaborative research within the aerodynamics community. The aircraft’s primary technical specifications typically encompass operating weights, performance characteristics, flight altitude, and powerplant details that are provided in Table 1. The characterization of the aircraft’s wing is primarily based on its cross-sectional shape (airfoil) and planform geometry. Key parameters that define the planform include the wing chord lengths at various positions along the semi-span, wing semi-span, leading-edge sweep angle, twist angle, dihedral angle, taper ratio, and aspect ratio. The original planform of the CRM wing corresponds to the 1g twisted flying shape under specific flight conditions. However, when this model is utilized for aerostructural analyses, the inherent 1g ‘built-in’ deflections present significant challenges. Consequently, it is essential to adopt the jig (undeflected) shape. The wing twist distribution from the side-of-body to the wingtip [16], along with the associated wing planform design parameters, are presented in Figure 1 and detailed in Table 1, respectively [33].

2.2. Structural and Finite Element Modeling of the CRM Wingbox

The external geometry of the wing is defined by CRM.65-BTE airfoil sections at three spanwise locations: the wing root, kink, and tip. The primary wing structure of the CRM is modeled to meet the minimum design criteria specified by the Federal Aviation Regulations (FAR) Part 25 [34]. A traditional two-spar wingbox architecture serves as the baseline design, with the wingbox derived from the wing surface model by positioning the front and rear spars at 12% and 71% of the local airfoil chord, respectively. The internal layout is defined by the stiffener pitch, rib pitch, and orientation, with values typical for large transport aircraft wings. Figure 2 illustrates the CRM wing surface model and the corresponding wingbox.

In the current study, the primary load-carrying wing structure is designed with high structural fidelity, as depicted in Figure 3. High-strength aluminum 7050-T7451 alloy [35] is used for the design of the upper skins, upper stiffeners, and spar caps of the wingbox, and 2024-T351 alloy [36] is used for the design of the lower skins, lower stiffeners, and ribs, as it is better suited for structures stressed by cyclic tension loads and therefore prone to fatigue damage. Composite materials made up of T300/N5208 carbon fibers and epoxy resin [37], which are widely used in the aircraft industry, are used as a second material choice for the wingbox structure’s design. For modeling the wingbox using a composite material, a symmetric and balanced laminate with ply orientation angles of [45/0/−45/90]s was created in order to obtain an orthotropic material. The aim of this design procedure was to avoid shear extension and membrane bending coupled behaviors.

The thin-walled structures of the wingbox configurations (skins, spar webs, and ribs) were modeled using two-dimensional quadrilateral and triangular shell elements (CQUAD4 and CTRIA3) with in-plane membrane and bending stiffness. On the other hand, stiffeners and spar caps were modeled using one-dimensional rod elements (CROD) with axial stiffness, as shown in Figure 4. The component masses of the leading and trailing edge devices are estimated based on the corresponding surface area using the semi-empirical and analytical equations of Torenbeek [38]. The inertial load impacts of these masses are modeled as lumped masses using point elements at the center of area of the leading and trailing edge devices, and they are attached to the front and rear spars of the wingbox along the span via multipoint constraint (MPC) non-stiffening rigid body elements (RBE3). The main assembly components of the landing gear in CRM transport aircraft are represented as concentrated lumped masses using point elements. These are integrated into the wingbox structure’s rear spar at the center of gravity using RBE3 elements. The engine mass is modeled as a concentrated lumped mass using point elements at the centre of gravity of the engine. For the engine pylon, a simple beam structure was created to realize a distributed engine pylon-to-wingbox connection. Rod elements (CROD), mounted on front and aft fittings, modeled using a combination of non-stiffening RBE3 and stiffening RBE2 rigid body-load-carrying elements, are used to model the engine and pylon structural components. The total fuel mass is modeled as concentrated lumped masses using point elements distributed along the wing fuel tanks at the center of gravity of all of the rib-bay volumes, depending on the user-defined filling levels and spanwise partitioning. The fuel lumped masses are connected to the wingbox’s lower skin using RBE3 elements. The breakdown of the masses per section between two rib locations is proportional to the volume of each section. The inertial load impact of the wingbox’s self-structural mass is derived by adding a downward gravitational acceleration (g = 9.81 m/s²) to the finite element model of the wingbox.

The aerodynamic loads were calculated for the symmetric pull-up maneuver at the limit load factor ($n$ = 2.5) at maximum take-off mass (260,000 kg) and design dive speed ($V_D$ = 221.7 m/s EAS, $M_D$ = 0.65) at sea-level conditions using ESDUpac A9510 utilizes steady lifting–surface theory based on the Multhopp–Richardson solution to calculate the spanwise loading of wings with camber and twist in subsonic attached flow [39]. This scenario was identified as the critical case for the design, analysis, and sizing optimization of the CRM wingbox, as well as for mass estimation. Consequently, the aerodynamic loads were distributed along the wing by computing the equivalent lift force and pitching moment components at the rib boundary locations, specifically at 25% of the local chord length. These loads were incorporated into the wingbox finite element model using RBE3 at the perimeter nodes of the ribs. Spring elements (CELAS1) combined with RBE2 elements were used to create realistic boundary conditions at the wingbox root at the aircraft centerline. The spring elements were attached to a fixed ground point. The translational and rotational stiffness properties were selected to result in end boundary conditions sufficiently close to the clamped case due to the lack of available data on wingbox root stiffness values for real aircraft structures in the open literature.

3. Aeroelastic Modeling and Flutter Analysis

The importance of aeroelasticity has been widely recognized by the aerospace industry. In recent years, the development of new types of aircraft has been explored, mainly by Boeing and Airbus. The new proposed aircraft are larger in size, quieter, greener, and more flexible than earlier generations of aircraft. By using lightweight composite materials, larger and more flexible wings have been created to generate lift with a minimum structural weight penalty. Aeroelastic effects due to wing flexibility may significantly alter the performance and safety of the new aircraft. Thus, considering aeroelastic effects at the early stages of the design is very important for producing a high-performance, competitive, and safe aircraft. The flutter phenomenon is generally accepted as a problem of primary concern for the design of current aircraft structures. Flutter is defined as a sustained oscillation of lifting surfaces, typically aircraft wings, and vertical and horizontal stabilizers, as a result of the high-speed passage of air along the lifting panels. It is related to the self-excited vibration present at certain forward flow speeds. For aircraft, flutter should only occur at speeds that are much higher than the operating speeds of the aircraft. EASA CS 25.629 and FAR Part 25.629 specify that the aircraft must be free from flutter with an appropriate margin of damping at all speeds up to 1.15 $V_D$, where $V_D$ is the design dive speed [34,40]. The total aerodynamic plus structural damping coefficient shall not be less than 3% ($g$ = 0.03) for any critical flutter mode at all altitudes and flight speeds from the minimum cruising speed up to the speed limit as defined by the flight envelope. Minimum damping and airspeed requirements are shown in Figure 5. For a given flight speed, damping is the modal parameter responsible for increasing the vibration amplitude, and it has been used as the index of flutter stability margins. If a mode exhibits damping characteristics similar to curve (1) in Figure 5, where the curve crosses the $g$ = 0 line above 1.15$V_D$, the curve is considered critical for flutter. The modes represented by curves (2) and (3) are classified as noncritical with respect to flutter. In the present study, the primary focus is on the critical mode shapes that are directly associated with the system’s flutter speed.

The equation of motion of a multi degree of freedom, discrete, and damped aeroelastic system can be derived based on the dynamic equilibrium of forces. The corresponding time-domain formulation, expressed in matrix form, is provided here for reference and further detail. The equation of motion of a multi degree of freedom, discrete, and damped aeroelastic system can be derived based on the dynamic equilibrium of forces. The time-domain equation of motion in matrix form is given as

(1)$\left[\mathit{M}\right]\ddot{\mathit{x}}\left(t\right)+\left[\mathit{C}\right]\dot{\mathit{x}}\left(t\right)+\left[\mathit{K}\right]\mathit{x}\left(t\right)=\mathit{F}(\mathit{x},t),$

where $\left[\mathit{M}\right]$, $\left[\mathit{C}\right]$, and $\left[\mathit{K}\right]$ represent the mass, damping, and stiffness matrices of the system, respectively, and $\mathit{x}\left(t\right)$ is the structural deformation vector. In general, the applied aerodynamic load vector $\mathit{F}\left(\mathit{x},t\right)$ is a time function of the structural deformation and the free stream Mach number $M_{\mathrm{\infty }}$ which is defined as

(2)$M_{\mathrm{\infty }}={\displaystyle \frac{V_{\mathrm{\infty }}}{a_s}},$

where $V_{\mathrm{\infty }}$ is the free stream velocity and $a_s$ is the speed of sound, which is a function of the flow temperature and density. The applied aerodynamic load vector $\mathit{F}(\mathit{x},t)$ consists mainly of two parts, as given by

(3)$\mathit{F}\left(\mathit{x},t\right)={\mathit{F}}_e\left(t\right)+{\mathit{F}}_a\left(\mathit{x}\left(t\right)\right),$

where ${\mathit{F}}_e\left(t\right)$ represents the externally applied non-aeroelastic forces, such as gust and control surface loads, and ${\mathit{F}}_a\left(\mathit{x}\left(t\right)\right)$ represents the aeroelastic forces, which are the induced aerodynamic forces due to the deformation of the structure. The aeroelastic forces are functions of flight speed and altitude, and the calculation of them relies on theoretical prediction methods that require unsteady aerodynamic computations. Because the aeroelastic forces are also functions of the structural deformation, Equation (1) can be written as

(4)$\left[\mathit{M}\right]\ddot{\mathit{x}}\left(t\right)+\left[\mathit{C}\right]\dot{\mathit{x}}\left(t\right)+\left[\mathit{K}\right]\mathit{x}\left(t\right)-{\mathit{F}}_a\left(\mathit{x}\left(t\right)\right)={\mathit{F}}_e\left(t\right).$

For aeroelastic stability analysis, the non-aeroelastic forces are ignored, resulting in the following equation:

(5)$\left[\mathit{M}\right]\ddot{\mathit{x}}\left(t\right)+\left[\mathit{C}\right]\dot{\mathit{x}}\left(t\right)+\left[\mathit{K}\right]\mathit{x}\left(t\right)-{\mathit{F}}_a\left(\mathit{x}\left(t\right)\right)=0.$

Equation (5) is a nonlinear time-domain equation that defines an aeroelastic structural system that can be self-excited in nature and gives rise to aeroelastic stability problems such as flutter and divergence.

The frequency domain p-k method, which is widely used in the aircraft industry, is utilized to perform flutter analysis using the MSC Flight Loads and Dynamics module using MSC Nastran flutter solution (Sol 145). One of the advantages of using the p-k method over other methods is that it produces results directly for the given velocity values, and the damping found from the p-k method equation represents a more realistic estimate of physical damping when compared to other methods. It also allows for the use of flutter results as the design responses for aeroelastic optimization. The fundamental matrix equation of the modal flutter solution is formulated as

(6)$\left[{\mathit{M}}_{hh}{p}^2+\left({\mathit{B}}_{hh}-{\displaystyle \frac{1}{4}}{\displaystyle \frac{\rho c{V}_{\infty }{\mathit{Q}}_{hh}^{Im}}{k}}\right)p+\left({\mathit{K}}_{hh}-{\displaystyle \frac{1}{2}}\rho V_{\infty }{\mathit{Q}}_{hh}^{Re}\right)\right]{\mathit{u}}_h=0,$

where ${\mathit{M}}_{hh}$, ${\mathit{B}}_{hh}$, and ${\mathit{K}}_{hh}$ are the modal mass, damping, and stiffness matrices, respectively, and ${\mathit{u}}_h$ is the modal amplitude vector. The unsteady aerodynamic loads are induced into the damping and stiffness matrices. The aerodynamic matrices are dependent on the reduced frequency $k$ but at a slow rate only. All of the matrices in Equation (6) are real; ${\mathit{Q}}_{hh}^{Re}$ and ${\mathit{Q}}_{hh}^{Im}$ are the real and imaginary parts of the complex aerodynamic matrix ${\mathit{Q}}_{hh}$, respectively. ${\mathit{Q}}_{hh}$ is dependent on the Mach number and the reduced frequency ${\mathit{Q}}_{hh}(M_{\infty }, k)$. In Equation (6), $k$ is the reduced frequency parameter and defined in terms of the system’s angular frequency, $\omega$, the reference chord length, $c$, and the selected free stream velocity $V_{\infty }$ as

(7)$k={\displaystyle \frac{\omega c}{2V_{\infty }}}.$

The eigenvalue $p$ parameter is defined in terms of the angular frequency $\omega$ and the coefficient of the transient decay rate $\gamma$ as

(8)$p=\omega \left(\gamma \mp i\right), \mathrm{where} i=\sqrt{-1 }.$

The transient decay rate is related to the structural damping coefficient $g$ by the following relation:

(9)$\gamma ={\displaystyle \frac{g}{2}}.$

Equation (6) is solved at a set of user-specified free stream velocities $V_{\infty }$ and air densities $\rho$, for the complex roots of the eigenvalue $p$ parameter with the modes of interest.

For the solution of the p-k method, the baseline equation, Equation (6), of the system can be represented in the state-space form as

(10)$\left(\left[\mathit{A}\right]-\left[\mathit{I}\right]p\right) {\mathit{u}}_h=0,$

where $\left[\mathit{A}\right]$ is the real matrix, $\left[\mathit{I}\right]$ is the identity matrix, and the vector ${\mathit{u}}_h$ includes both the modal amplitudes and the velocities:

(11)$\left[\mathit{A}\right]=\left[\begin{array}{cc}0& \mathit{I}\\ -{\mathit{M}}_{hh}^{-1}\left({\mathit{K}}_{hh}-{\displaystyle \frac{1}{2}}\rho {V_{\infty }}^2{\mathit{Q}}_{hh}^{Re}\right)& -{\mathit{M}}_{hh}^{-1}\left({\mathit{B}}_{hh}-{\displaystyle \frac{1}{4}}{\displaystyle \frac{\rho c{V}_{\infty }{\mathit{Q}}_{hh}^{Im}}{k}}\right)\end{array}\right].$

Equation (6) is solved at several given values of velocities $V_{\infty }$ and air densities $\rho$ for the complex roots of the eigenvalue parameter $p$ associated with the modes of interest. This is achieved through an iterative solution that matches the reduced frequency $k$ to the imaginary part of the eigenvalue parameter $p$ for every structural mode. Plots of $V_{\infty }$ versus $g$ can then be used to determine the flutter speed. Flutter occurs for values of $M_{\infty }$, $k$, and $\rho$ for which $g$ = 0, where $g$ goes from negative to positive values, indicating instability.

The natural frequency values of interest are the input values for the calculation of reduced frequencies k for the flutter analysis. MSC Nastran solution Sol 103 is used to perform natural frequency and normal mode analysis. In MSC Nastran, the Lanczos method is suggested to perform real eigenvalue extraction. It is the preferred method for most medium-to-large-sized finite element models. The next step is to calculate the aerodynamic matrices ${\mathit{Q}}_{hh}(M_{\infty }, k)$. This requires the calculation of the unsteady aerodynamic forces. These are calculated in the frequency domain for a discrete set of reduced frequencies based on the assumption of the undamped harmonic motion using the Doublet Lattice Method (DLM). The DLM is based on the linearized aerodynamic potential theory of subsonic flow and offers low-order models for unsteady aerodynamics in the subsonic regime, which lead to aeroelastic models suitable for flutter analysis. The DLM models the lifting surface, such as the wing of an aircraft, as an aerodynamic flat panel parallel to the flow. It does not account for the airfoil thickness, camber or pre-twist along the span of the wing. The aerodynamic panel is divided into small boxes arranged in strips parallel to the free stream velocity $V_{\mathrm{\infty }}$ with surface edges on the box boundaries. The aerodynamic mesh must comply with the specific criteria and guidelines as defined in [41,42] in order to achieve a sufficient level of accuracy, including requirements for the number of spanwise strips, the magnitude of the box aspect ratios and the numbers of chordwise boxes determined by the reduced frequency. The DLM mesh is composed of 15 x 80 boxes in chordwise and spanwise directions, as shown in Figure 6. The structural finite element model of the CRM wing is connected to an aerodynamic model by means of splines, as described in [41,42]. The Finite Plate Spline (FPS) fixes points on the structural model to points on the aerodynamic model so that the loads and deformations can be transferred at those points. Figure 7 and Figure 8 show the selected spline points of the structural model and the coupling between the aerodynamic and structural grids achieved using surface splines, respectively. The final step is the flutter analysis, carried out using a non-matched flutter analysis approach. This type of non-matched analysis is usually applied for the certification of subsonic aircraft. In this approach, the aerodynamic matrices ${\mathit{Q}}_{hh}\left(M_{\infty }, k\right)$ are calculated for one reference value of the Mach number $M_{Ref}$ and for a set of reduced frequencies at fixed altitude. The reference Mach number is set to $M_D$, which is the maximum Mach number that occurs in the aircraft flight envelope. Therefore, the analysis velocities do not match the reference Mach number and the velocity results over $V_D$ are artificial. These results represent the flutter stability rate of reserve with respect to the flutter speed requirements.

4. Definition of the CRM Wingbox Optimization Problem

The structural optimization of the CRM wingbox presented in this study is primarily focused on property optimization. In this context, the geometry of the wing, including the rib count and the positions of the spars and stiffeners, is held fixed throughout the sizing optimization process, with no shape optimization being performed. The primary objective of the study is to minimize the masses of the metallic and composite configurations of the CRM wingbox while adhering to a set of design constraints, which are detailed in the following subsections The optimization process utilizes the MSC Nastran gradient-based Sol 200 optimizer, which employs gradient-based algorithms to solve the problem. A key advantage of using gradient-based methods is their efficiency in handling large design spaces, where the number of design variables significantly exceeds the number of objectives and constraints. Additionally, these methods are computationally efficient, offering rapid convergence with clear convergence criteria. However, a notable limitation of gradient-based approaches is the potential for multiple local optima, making it difficult to guarantee global optimality. To address this, global optimality is pursued by initiating random searches from various starting points within the design space. Nevertheless, for large nonlinear optimization problems that combine continuous and discrete variables, this approach can become slow and computationally inefficient, particularly when a single optimizer run fails to converge to a feasible solution. To address this challenge, the practical optimization procedure developed by Dababneh et al. [8] in an industrial context was integrated into the sizing optimization process. This procedure incorporates enhanced features for solving large-scale nonlinear structural optimization problems, with a key innovation being an efficient method for generating initial design variable values. For additional details, the reader is encouraged to refer to [8].

4.1. Objective Function

The objective function is the structural mass of the CRM wingbox excluding any non-structural masses, like fuel mass, landing gear mass and engine mass. The objective function can be represented by

(12)$\mathrm{minimize} M\left(\mathit{x}\right)=\sum _{l=1}^{N element}{\rho }_l{V}_l\left(\mathit{x}\right).$

Here, the objective function $M\left(\mathit{x}\right)$ represents the wingbox structural mass, while $N element$ is the number of elements in the finite element model, $V_l$ is the volume of the $l\mathrm{t}\mathrm{h}$ element, ${\rho }_l$ is the corresponding material density, and $\mathit{x}$ is the design variable vector.

4.2. Design Variables

For the optimization problem, considering the wingbox construction material to be a metallic material, one design variable per design field is defined, as described in Figure 3, Section 2.2, where the model is discretized into components which act as design optimization zones along the span, totaling 1870 design zones. The chordwise design zones are prescribed by the stringer pitch, while in the spanwise direction the design zones are limited by the rib spacing. In the finite element model, each design field consists of a number of finite elements that all comprise the same thicknesses/cross-sectional areas and stiffness properties. The design variables include the thicknesses of the wingbox skins, spar webs and ribs, as well as the cross-sectional areas of the wingbox spar caps and stiffeners. A minimum gauge thickness of 2 mm and a cross-sectional area of 144 mm² are specified for the design variables. The limits on the design variables are defined as follows:

(13)$2.0 \le t_{metallic}\le 13.5, \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e} t_{metallic} \mathrm{i}\mathrm{s} \mathrm{t}\mathrm{h}\mathrm{e} \mathrm{t}\mathrm{h}\mathrm{i}\mathrm{c}\mathrm{k}\mathrm{n}\mathrm{e}\mathrm{s}\mathrm{s} \mathrm{i}\mathrm{n} \mathrm{m}\mathrm{m},$

(14)$144.0\le a_{metallic}\le 972.0, \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e} a_{metallic}\mathrm{i}\mathrm{s} \mathrm{t}\mathrm{h}\mathrm{e} \mathrm{a}\mathrm{r}\mathrm{e}\mathrm{a} \mathrm{i}\mathrm{n} {\mathrm{m}\mathrm{m}}^2.$

On the other hand, considering the wingbox construction material to be a composite material, the corresponding design variables for the wingbox skins, spar webs and ribs are the thicknesses of each ply or lamina in the composite laminate associated with each design field. The cross-sectional areas of the composite spar caps and stiffeners are also treated as individual design variables for each design zone. For modeling the wingbox using a composite material, a symmetric and balanced laminate with ply orientation angles of [45/0/−45/90]s was created to get an orthotropic material and Laminate design rules were applied were angle differences between plies should be minimized to minimize interlaminar shear stresses caused by shear coupling. The number of plies in each ply orientation angle is determined based on the optimization constraints. These ply thicknesses are varied discretely during the optimization process. The schematic of the composite laminate is shown in Figure 9 for completeness.

The minimum ply thickness is taken to be 0.127 mm, while a 3 mm minimum gauge laminate thickness is recommended to maintain an adequate level of laminate damage tolerance. The laminate ply thicknesses are treated as individual design variables and a count is made of the required number of plies in each ply orientation angle. The limits on the number of plies in each ply orientation angle are given as

(15)$3 \le n_{ply}\le 60, \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e} n_{ply} \mathrm{i}\mathrm{s} \mathrm{t}\mathrm{h}\mathrm{e} \mathrm{n}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r} \mathrm{o}\mathrm{f} \mathrm{p}\mathrm{l}\mathrm{i}\mathrm{e}\mathrm{s}.$

Minimum cross-sectional areas of 216 mm² for the composite spar caps and stiffeners are specified and the limits on the design variables are defined as follows:

(16)$216.0\le a_{composite}\le 972.0, \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e} a_{composite} \mathrm{i}\mathrm{s} \mathrm{t}\mathrm{h}\mathrm{e} \mathrm{a}\mathrm{r}\mathrm{e}\mathrm{a} \mathrm{i}\mathrm{n} {\mathrm{m}\mathrm{m}}^2.$

4.3. Static Strength Constraints

For metallic skin panels, spar webs and ribs, the von Mises stress is checked against the material allowable stress as defined in the following equation:

(17)${\sigma }_{von Mises}\le {\sigma }_{allowable}$

For composite skin panels, spar webs, and ribs, the Tsai–Wu criterion [43,44,45] is used to predict the strength of the composite laminate in terms of the failure index ($FI$). For orthotropic plate analysis, under the plane stress state, the Tsai–Wu strength theory predicts that a lamina will undergo failure when the following inequality in satisfied:

(18)$FI=F_1{\sigma }_1+F_2{\sigma }_2+F_{11}{\sigma }_1^2+2F_{12}{\sigma }_1{\sigma }_2+F_{22}{\sigma }_2^2+F_{66}{\sigma }_6^2\ge 1.$

The coefficients $F_1$–$F_{66}$, with the exception of $F_{12}$, are described in terms of strengths in the principal material directions. $F_{12}$ accounts for the interaction between normal stresses, ${\sigma }_1$, and ${\sigma }_2$.

The principal strains in each ply are also checked against the material allowable strain to ensure the integrity of the plies and failure-free laminates. The allowable strain value of 3500$\mu \epsilon$ includes the margins due to fatigue and damage tolerance, assuming that the allowable strains are identical in terms of tension and compression. Thus, the following constraint is placed on the strain value used for sizing the structure:

(19)${\epsilon }_{principal}\le {\epsilon }_{allowable}$

The spar caps and the longitudinal stiffeners are designed to carry axial stress only. Therefore, they are designed according to their stress state against the allowable stress of the material as defined in the following equation:

(20)${\sigma }_{axial}\le {\sigma }_{allowable}$

4.4. Static Stiffness Constraints

The flexural stiffness of the wingbox is controlled by limiting the vertical displacement of the wingtip leading edge [46,47]. The wingtip deflection ${\delta }_{tip\left(Z\right)}$ for the CRM wing at a 2.5 g pull-up maneuver is assumed to be 17% of the wing semi-span b.

(21)${\delta }_{tip\left(Z\right)}\le 17\%·b, \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e} b \mathrm{i}\mathrm{s} \mathrm{t}\mathrm{h}\mathrm{e} \mathrm{w}\mathrm{i}\mathrm{n}\mathrm{g} \mathrm{s}\mathrm{e}\mathrm{m}\mathrm{i}-\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}.$

The torsional stiffness, which is necessary to counteract the twisting of the wing under aerodynamic loads and thus prevents flutter, is controlled by constraining the twist angle at the tip chord of the wing. The angular deformation at the wingtip chord is constrained by limiting it to a value of 6° to ensure sufficient torsional stiffness and thus an adequate aeroelastic response [48]. The twist angle constraint is defined using the vertical displacements at the wingtip chord ends. Equation (17) shows that the twist angle at the wingtip should not exceed 6°. ${\left(\delta \right)}_{max}^{+}$ and $({\delta )}_{max}^{-}$ are the maximum vertical displacements in positive and negative directions of the z-coordinate, respectively. Here, $C$ is the wing chord length at the required location:

(22)${\Theta }_{tip}\le 6°, \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e} \Theta =\arctan \left({\displaystyle \frac{{\left(\delta \right)}_{max}^{+}-({\delta )}_{max}^{-}}{C}}\right).$

4.5. Manufacturing Constraints

The design and optimization process of the CRM wingbox incorporates practical design rules and manufacturing constraints. To facilitate manufacturing, the thicknesses of the wingbox’s thin panels are selected from a discrete set of values ranging from the lower to upper bounds, with increments of 0.1 mm. Specifically, the metallic panel thickness is chosen from the set as:

(23)$t_{metallic}\in \left\{2.0,2.1,\dots ,20\right\}, \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e} t_{metallic} \mathrm{i}\mathrm{s} \mathrm{t}\mathrm{h}\mathrm{e} \mathrm{t}\mathrm{h}\mathrm{i}\mathrm{c}\mathrm{k}\mathrm{n}\mathrm{e}\mathrm{s}\mathrm{s} \mathrm{i}\mathrm{n} \mathrm{m}\mathrm{m}.$

In this study, rod elements were used as a simplified means to represent axial load-carrying members primarily for the purpose of evaluating relative sizing effects and distribution within the optimization framework. The composite material properties assigned to these elements represent equivalent axial stiffness values based on the dominant fiber direction, but do not capture the full complexity of the laminate behavior or local effects, such as buckling or delamination. Regarding the stiffener cross-section, the axial stiffness of the rod elements was updated during the optimization process, reflecting changes in the cross-sectional area. However, the detailed geometry was not explicitly modeled, as the focus of this work was on overall structural trends and load path optimization rather than detailed stringer design. For the metallic and composite spar caps and stiffeners, rod elements are used in the modeling process. To simplify manufacturing, the rods are sized from a discrete set of cross-sectional areas, allowing for the production of various flange shapes, such as L, T, and Z. The cross-sectional areas for the metallic and composite rod elements are defined as:

(24)$a_{metallic}\in \left\{144,158,\dots ,972\right\}, \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e} a_{metallic} \mathrm{i}\mathrm{s} \mathrm{t}\mathrm{h}\mathrm{e} \mathrm{a}\mathrm{r}\mathrm{e}\mathrm{a} \mathrm{i}\mathrm{n} {\mathrm{m}\mathrm{m}}^2,$

(25)$a_{composite}\in \left\{216,223,\dots ,972\right\}, \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e} a_{composite} \mathrm{i}\mathrm{s} \mathrm{t}\mathrm{h}\mathrm{e} \mathrm{a}\mathrm{r}\mathrm{e}\mathrm{a} \mathrm{i}\mathrm{n} {\mathrm{m}\mathrm{m}}^2.$

In addition, the number of plies in each ply orientation is constrained to discrete integer values, ranging from a lower bound to an upper bound and incremented by 1 as

(26)$n_{ply}\in \left\{3,4,\dots ,60\right\}, \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e} n_{ply} \mathrm{i}\mathrm{s} \mathrm{t}\mathrm{h}\mathrm{e} \mathrm{n}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r} \mathrm{o}\mathrm{f} \mathrm{p}\mathrm{l}\mathrm{i}\mathrm{e}\mathrm{s}.$

Furthermore, practical design criteria are applied to the optimization of the composite laminate wingbox structures. The ply orientation percentages within a laminate are constrained between 10% and 60% to avoid matrix-dominated behaviors. An additional optimization constraint ensures that the thicknesses of the +45° and −45° layers are identical, promoting a balanced laminate structure and minimizing the introduction of manufacturing stresses, such as torsion. A maximum property drop-off rate criterion is applied. It aims on the one hand at avoiding delamination and, on the other hand, at obtaining ply layouts that can actually be manufactured. The property drop-off rate between neighboring elements/panels is evaluated according to the following equation [49]:

(27)$Property drop-off rate={\displaystyle \frac{{prop}_1-{prop}_2}{distance\left(d\right)}}\le 20\%,$

where ${prop}_{\mathrm{i}}$ is the element/panel property value of the parent (1) or adjacent (2) element/panel and the distance $d$ is computed along the element/panel surfaces between adjacent centroids, as shown in Figure 10 [49].

4.6. Aeroelastic Stability Constraints

The dynamic aeroelastic stability constraint (flutter) is imposed by constraining the damping rather than the flutter speed. Defining the constraint in such a way eliminates the need for computation of the flutter speed. Exact computation of this speed is a computationally expensive task. Treating the aeroelastic constraints in this manner was first proposed by Hajela [50] and has become a standard process in MSC Nastran [51]. Thus, the flutter constraint is defined as follows:

(28)$g={\displaystyle \frac{{\gamma }_{jl}-{\gamma }_{jREQ}}{GFACT}}\le 0,j=\mathrm{1,2},\dots ,nv,l=\mathrm{1,2},\dots ,nroot$

where${\gamma }_{jl}$ is the damping for the $l\mathrm{t}\mathrm{h}$ root calculated at the $j\mathrm{t}\mathrm{h}$ velocity and ${\gamma }_{jREQ}$ is the user-defined required damping value at the $j\mathrm{t}\mathrm{h}$ velocity, typically 0.03. $GFACT$ is a scale factor that converts the damping numbers into a range of suggested values, typically in the range of 0.1–0.5, which is consistent with other constraints in the design task. A typical $GFACT$ value of 0.1 is used in this study [51,52]. The flutter constraint is to be satisfied at a series of velocities up to, and may be above, the required flutter speed.

The p-k method of flutter analysis produces solutions only at the velocities of interest. The evaluation of the flutter constraint is performed at a number of velocities to handle the development of ‘hump’ modes, as at velocities lower than the required flutter speed these modes could become critical. In the context of flutter analysis, hump modes are identified by local maxima in flutter speed versus frequency curves, as documented in [28,29]. These modes are of particular interest as they may indicate conditions under which the structure exhibits heightened sensitivity to aeroelastic instabilities at specific frequencies. Analytical identification of hump modes typically involves evaluating the system’s frequency response and locating points where such local maxima occur, which may signal potential non-physical numerical artifacts or critical dynamic behavior depending on the modeling approach. Hence, in the current study, if the flutter analysis is performed at flight speeds that are 1.1, 1.0, 0.9, 0.75, and 0.5 times the required speed, the results should be sufficient to prevent this unwanted behavior.

5. Results and Discussion of CRM Wingbox Optimization Case Studies

To assess the impact of integrating the dynamic aeroelastic stability constraint (flutter) during the preliminary design stages, the CRM wing is initially optimized with static strength, stiffness, and manufacturing constraints. It is then re-optimized to meet dynamic aeroelastic stability in addition to the static and manufacturing constraints. A comparative analysis of the optimization results, including the optimized masses and flutter speeds, is performed to evaluate the impact of this integration.

5.1. Optimization Results of the CRM Wingbox Subject to Static Strength, Stiffness and Manufacturing Constraints

The optimized mass values of the CRM wingbox for both metallic and composite configurations, as well as for continuous and discrete optimization approaches, are presented in Table 2 and Table 3. The entries in bold within these tables indicate the local minimum feasible optimized solutions, which were derived using discrete values for the design variables. However, it is important to note that only one feasible discrete design solution was identified for the composite CRM wing, with the remaining discrete design solutions being infeasible. The primary cause of the infeasibility in the discrete designs was traced to the response output associated with the maximum constraint violation. This was determined through an analysis of the design response results, which revealed that the minimum ply orientation percentage for several design variables fell below the minimum allowable threshold of 10%. Consequently, these designs were deemed infeasible. This outcome underscores the critical importance of employing diverse initial guesses for the design variables to attain feasible local optimum solutions in gradient-based optimization techniques.

Table 4 presents a summary of the maximum tip displacement values for the optimized CRM wingbox models, comparing both metallic and composite configurations, along with the change in the wingtip twist angle relative to the baseline wing model. As indicated in Table 4, the use of metallic materials in the design of the CRM wing results in a significantly higher wingtip deflection compared to the composite material configuration. This observation suggests that the flexural stiffness of the metallic wing is lower than that of the composite wing, leading to greater deflection in the metallic wing. This higher deflection raises concerns regarding the performance of the wing’s ailerons, the dynamic aeroelastic effects such as flutter, and potential negative impacts on aerodynamic performance, including lift reduction. For both wingbox configurations, the change in wingtip twist from the baseline model to the optimized model indicates an upward twist of the wingtip. Notably, the composite wingbox exhibits a greater twist than the metallic wingbox, suggesting that the torsional stiffness of the composite wingbox is lower than that of the metallic wingbox. This issue could be mitigated by incorporating a higher proportion of 45° fibers in the laminate, and by imposing a lower bound for the minimum 45° ply orientation percentage in the laminate, which could exceed the 10% value used in the optimization study.

5.2. Flutter Analysis the Optimized CRM Wingbox Subject to Static Strength, Stiffness and Manufacturing Constraints

A free vibration analysis is conducted to determine the natural frequencies and mode shapes of the optimized CRM wingbox using the normal mode analysis module in MSC Nastran with the Lanczos method [53]. The resulting global mode shapes and natural frequencies for the optimized metallic and composite CRM wingboxes are presented in Table 5 and Table 6, respectively. The analysis is performed under both full fuel and zero fuel loading conditions. The comparison of global mode shapes shows that the natural frequencies of the optimized composite wingbox are higher than those of the metallic wingbox. Moreover, the presence of fuel masses in the inboard and outboard regions reduces the structural frequency associated with wing bending, compared to the empty fuel tanks condition. Overall, fuel quantity significantly influences both flutter modes and speed, highlighting the need to assess the most critical conditions.

An aeroelastic stability analysis is performed for the optimized CRM wingbox using the p-k method under sea-level conditions. Flutter instability in the subsonic regime is identified through non-matched flutter analysis. Aerodynamic forces are computed for a reference Mach number $M_{Ref}$ and a set of reduced frequencies at a fixed altitude. The reference values for the flutter critical flight altitude and Mach number are typically selected at sea level and at the maximum design dive speed ($V_D$ = 221.7 m/s EAS), respectively. The Doublet Lattice Method integrated in MSC Nastran is used to calculate the unsteady aerodynamic forces. For the optimized metallic CRM wingbox, the flutter speed with full fuel tanks is determined to be 370 m/s, with a corresponding flutter frequency of 4.42 Hz, while with zero fuel, the flutter speed is 306 m/s and the flutter frequency is 4.18 Hz. The damping and frequency curves for the metallic CRM wingbox are shown in Figure 11 and Figure 12, respectively. In contrast, for the optimized composite CRM wingbox, the flutter speed with full fuel tanks is 428 m/s, with a corresponding flutter frequency of 6.42 Hz, and with zero fuel, the flutter speed is 394 m/s, with a corresponding flutter frequency of 5.80 Hz. The damping and frequency curves for the composite CRM wingbox are presented in Figure 13 and Figure 14, respectively. The results demonstrate that both the optimized metallic and composite CRM wings are free from flutter instabilities, with flutter speeds exceeding ($V_{Flutter}$ > 256 m/s). Additionally, the flutter speed of the optimized composite wing is found to be higher than that of the metallic wing. These flutter speed results, which are influenced by static displacement, align with the displacement data presented in Table 4. The higher displacement values provide insight into the flexural stiffness of the wingbox. Moreover, the flutter speed is shown to be highly sensitive to the fuel mass in the wing fuel tanks. The lowest flutter speed, which is critical for design, occurs under empty fuel tank conditions. The flutter analysis confirms that the stiffness constraints are appropriately applied, ensuring the wing remains free from flutter. The wingbox stiffness was controlled by limiting both the displacement and twist of the wingtip, thereby ensuring an adequate aeroelastic response.

5.3. Optimization Results of the CRM Wingbox Subject to Static Strength, Stiffness, Aeroelastic and Manufacturing Constraints

To enable a direct comparison with the results presented in the previous section, the optimization was conducted using the same initial conditions that produced a feasible optimal discrete solution. The CRM wing was optimized to simultaneously satisfy the static strength and stiffness constraints, the dynamic aeroelastic constraint, and the manufacturing constraints. Table 7 gives the optimized masses of the CRM wingbox for both the metallic and composite configurations and for both the continuous and discrete optimization solutions. During the optimization process, hard convergence is obtained along with a hard feasible discrete design solution.

Table 8 summarizes the maximum tip displacement value and the change in wingtip twist angle from the baseline wing model in the optimized CRM wingbox model for both the metallic and composite configurations.

The results from Table 7 show that the mass of the optimized metallic CRM wingbox has increased as compared to the results obtained by only considering static strength, stiffness and manufacturing constraints only. This is not the case for the optimized composite CRM wingbox. The deformation results presented in Table 8 show that the wingtip deflection as well as the twist angle of the optimized metallic CRM wingbox has decreased and as a result the wingbox structural mass has increased. On the other hand, while the wingtip deflection of the optimized composite has decreased, the wingtip twist angle has increased. The increase in the wingtip twist indicates that the torsional stiffness of the optimized composite CRM wingbox has decreased and as a result reduced the structural mass of the composite wingbox.

Figure 15 and Figure 16 present fringe plots of the von Mises stress and axial stress distributions, respectively, for the optimized metallic CRM wingbox. The stress contours in both figures indicate that all stress values lie within the allowable limits defined within the design optimization domain, thereby satisfying the stress constraints imposed during the structural optimization process. Figure 15 reveals localized peak stress concentrations in the vicinity of the rear spar, specifically near the attachment points of concentrated lumped masses. As established in finite element analysis literature [54], the accuracy of stress approximations is generally higher at the element centroids compared to the nodal regions. Consequently, these localized peak stresses are attributed to numerical artifacts arising from mesh discretization and are not indicative of physically meaningful structural behaviour; thus, they are excluded from the assessment of global structural performance. Figure 16 further identifies the locations of maximum axial stresses, both tensile and compressive. Elevated compressive stresses are observed in the upper flanges, while moderately lower tensile stresses appear in the lower flanges. This distribution is consistent with the expected structural response of the wingbox under the applied loading conditions, reflecting proper load transfer and structural integrity.

The principal strain and axial stress distributions were analyzed for each optimized composite laminate layer to evaluate structural performance and ensure compliance with design constraints. Figure 17 and Figure 18 illustrate the fringe plots of the maximum and minimum principal strain distributions, respectively, for the optimized composite CRM wingbox model. As observed, the principal strain values in both plots remain below the prescribed upper limit of 3500 με, confirming that the strain constraint is satisfied throughout the structure.

The resulting global mode shapes and the corresponding natural frequencies of the optimized metallic and composite CRM wingboxes are presented in Table 9 and Table 10 and Figure 19 and Figure 20, respectively. The analysis was performed on the CRM wingbox models with full fuel and zero fuel loading.

An aeroelastic stability analysis is conducted for the optimized CRM wingbox models to ensure that the design is free from flutter. For the metallic CRM wingbox, the flutter speed with full fuel tanks is 399 m/s, corresponding to a flutter frequency of 6.87 Hz. With zero fuel, the flutter speed is 351 m/s, with a flutter frequency of 5.00 Hz. This results in a 14.7% increase in the critical flutter speed at zero fuel, achieved at the cost of a 1.8% increase in the total wingbox structural mass. The damping and frequency curves for the optimized metallic CRM wingbox at the lowest flutter speed are presented in Figure 21.

For the optimized composite CRM wingbox, the flutter speed with full fuel tanks is 432 m/s, with a flutter frequency of 8.10 Hz. With zero fuel, the flutter speed is 373 m/s, and the flutter frequency is 6.04 Hz. This represents a 5.6% reduction in the total wingbox structural mass, accompanied by a 5.3% decrease in the critical flutter speed at zero fuel. Despite the reduction in flutter speed, the wing remains free from flutter. The composite material’s added versatility has allowed the optimization process to tailor the stiffness distribution, preventing flutter and achieving a minimum mass design by controlling the thickness of each lamina in the wingbox structural panels. The damping and frequency curves for the optimized composite CRM wingbox at the lowest flutter speed are shown in Figure 22.

6. Concluding Remarks

Based on the results presented in this study, the following points could be concluded:

The optimization of the CRM wingbox, which first incorporated static strength, stiffness, and manufacturing constraints and subsequently re-optimized to include dynamic aeroelastic stability, provides key insights into the significance of aeroelastic constraints in the preliminary design phase. The study demonstrates that integrating dynamic stability constraints, such as flutter, profoundly impacts the optimization process, particularly in terms of wing mass and flutter speeds.

The feasibility of incorporating composite materials into the CRM wingbox design demonstrates that the composite wing results in a structural mass reduction of approximately 17.4% relative to the metallic wingbox when flutter constraints are considered, compared to a 23.4% mass reduction when flutter constraints are not imposed.

When accounting for flutter constraints, the findings indicate a 5.6% reduction in the structural mass of the composite CRM wingbox, accompanied by a 5.3% decrease in the critical flutter speed. Despite the reduction in flutter speed, the design remains free from flutter instabilities within the operational flight envelope. Furthermore, the free vibration and aeroelastic stability analyses indicate that the composite wingbox exhibits higher natural frequencies than its metallic counterpart, suggesting that composite materials enhance the dynamic response and reduce the susceptibility to aeroelastic phenomena.

The flutter analysis, performed using the p-k method under sea-level conditions, confirmed that both the metallic and composite wingboxes are free from flutter instabilities, with flutter speeds exceeding the critical threshold of 256 m/s. The composite wingbox demonstrated a higher flutter speed compared to the metallic wingbox, consistent with its increased stiffness, as indicated by the displacement data. Additionally, the analysis identified that the lowest flutter speed occurs when the fuel tanks are empty, representing the most critical condition for flutter design.

Fuel mass was found to have a significant impact on both the natural frequencies and flutter characteristics, with the presence of fuel resulting in a reduction of structural frequencies associated with wing bending. This highlights the importance of incorporating varying fuel loads into the design process, as fuel mass can influence both aeroelastic performance and flutter speed.

7. Future Work Recommendations

Building on the findings, several recommendations are made to enhance the accuracy of wing mass estimation through the integration of multidisciplinary analysis and design optimization techniques. Specifically, the use of advanced optimization algorithms is advised to overcome the issue of local optima, facilitating a broader search for solutions and minimizing the risk of converging on suboptimal results. Future work will focus on deriving spanwise bending and torsional rigidity properties, specifically the EI and GJ distributions, to accurately capture stiffness variations along the wingspan. Plotting these variations will provide valuable insights into the structural behavior for aeroelastic modeling. Additionally, the development of a reduced-order stick-beam model will streamline the structure while preserving essential flexibility characteristics, enhancing the accuracy of aeroelastic simulations.

The optimization model could be extended to include additional constraints, such as fundamental frequency, local buckling, divergence speed, and thermal stresses. Fatigue failure can occur at stress levels significantly below the material’s yield or ultimate strength. Therefore, the application of a knock-down factor is appropriate to account for fatigue effects under different flight conditions, ensuring a more conservative and fatigue-tolerant structural design. Future investigations could explore the application of MDO to a hybrid CRM wingbox structure that combines both composite and metallic materials in the primary wing components. Additionally, the optimization model could be enhanced to incorporate multi-objective aerodynamic and structural optimization, including constraints on the wing’s leading edge sweep angle, taper ratio, and aspect ratio. This extension would enable the identification of multiple Pareto-optimal solutions, thereby providing valuable design alternatives for the preliminary design phase of the CRM wing.

Author Contributions

O.D. conducted all simulations and drafted the manuscript, while T.K. and J.F.W. provided guidance and feedback on the presented results. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 Comparison of the twist distribution between the CRM and the uCRM (undeflected shape) [16].

Figure 2 Surface and wingbox model of the CRM wing.

Figure 3 High-fidelity wingbox structural model of the CRM wing.

Figure 4 Finite element model of the CRM wingbox structure.

Figure 5 Minimum damping and airspeed requirements.

Figure 6 CRM wing DLM aerodynamic mesh.

Figure 7 CRM wing structural model selected spline points.

Figure 8 Structural and aerodynamic meshes coupling. Structural (green) and aerodynamic mesh (blue).

Figure 9 Schematic of the composite laminate.

Figure 10 Control of property taper rate.

Figure 11 Velocity versus damping and velocity versus frequency plots of the optimized metallic CRM wingbox—full fuel.

Figure 12 Velocity versus damping and velocity versus frequency plots of the optimized metallic CRM wingbox—zero fuel.

Figure 13 Velocity versus damping and velocity versus frequency plots of the optimized composite CRM wingbox—full fuel.

Figure 14 Velocity versus damping and velocity versus frequency plots of the optimized composite CRM wingbox—zero fuel.

Figure 15 von Mises (MPa) distributions on the thin panels of the optimized metallic CRM wingbox mode.

Figure 16 Axial stress (MPa) distributions on flanges of the optimized metallic CRM wingbox model.

Figure 17 Maximum principal strain distribution (με)—layer 1.

Figure 18 Minimum principal strain distribution (με)—layer 1.

Figure 19 Global mode shapes of the optimized metallic CRM wingbox (flutter constraint included)—full fuel.

Figure 20 Global mode shapes of the optimized composite CRM wingbox (flutter constraint included)—Full Fuel.

Figure 21 Velocity versus damping and velocity versus frequency plots of the optimized metallic CRM wingbox (flutter constraint included).

Figure 22 Velocity versus damping and velocity versus frequency plots of the optimized composite CRM wingbox (flutter constraint included).

Table 1

Relevant data with respect to NASA CRM aircraft [33].

Description	Value
Max. take-off mass	260,000 kg
Max. zero fuel mass	19,500 kg
Main landing gear mass	9620 kg
Engine mass (2×)	15,312 kg
Max. fuel mass	131,456 kg
Wing gross area	383.7 m2
Wingspan	58.76 m
Aspect ratio	9.0
Root chord	13.56 m
Tip chord	2.73 m
Taper ratio	0.275
Leading edge sweep	35.0°
Cruise speed	193.0 m/s EAS
Dive speed	221.7 m/s EAS
Cruise altitude	10,668 m

Table 2

Optimized masses (kg) of the metallic CRM wingbox.

Design Variable Initial Values
Min	25% Max	50% Max	75% Max	Max
Continuous Solution
13,168	13,422	13,303	13,388	13,342
Discrete Solution
13,265	13,523	13,401	13,486	13,441

Table 3

Optimized masses (kg) of the composite CRM wingbox.

Design Variable Initial Values
Min	25% Max	50% Max	75% Max	Max
Continuous Solution
10,420	10,528	10,869	11,031	11,256
Discrete Solution
10,837	10,954	11,289	11,465	11,665

Table 4

Deformation values of the optimized CRM wingbox models.

Deformation	Metallic	Composite
Displacement ${\mathit{\delta }}_{\mathit{t}\mathit{i}\mathit{p}\left(\mathit{Z}\right)}$	4.99 m	3.48 m
Twist ${\mathit{\Theta }}_{\mathit{t}\mathit{i}\mathit{p}}$	0.42°	3.77°

Table 5

Global mode shapes and the associated frequencies of the optimized metallic CRM wingbox.

Mode Number	Mode Shape Description	Natural Frequency [Hz](Full Fuel)	Natural Frequency [Hz](Zero Fuel)
1	1st Out-of-Plane Bending	1.07	1.65
2	2nd Out-of-Plane Bending	3.12	4.28
3	1st Lateral Bending	3.39	5.51
4	3rd Out-of-Plane Bending	4.57	8.79
5	1st Torsion	19.02	19.19

Table 6

Global mode shapes and the associated frequencies of the optimized composite CRM wingbox.

Mode Number	Mode Shape Description	Natural Frequency [Hz](Full Fuel)	Natural Frequency [Hz](Zero Fuel)
1	1st Out-of-Plane Bending	1.39	2.28
2	2nd Out-of-Plane Bending	3.72	5.63
3	1st Lateral Bending	4.44	7.52
4	3rd Out-of-Plane Bending	6.21	11.01
5	1st Torsion	23.40	24.79

Table 7

Optimized masses (kg) of the CRM wingbox (flutter constraint included).

Design Variable Initial Values
Metallic (Min)	Composite (25% Max)
Continuous Solution
13,393	9931
Discrete Solution
13,498	10,340

Table 8

Deformation values of the optimized CRM wingbox models (flutter constraint included).

Deformation	Metallic	Composite
Displacement ${\mathit{\delta }}_{\mathit{t}\mathit{i}\mathit{p}\left(\mathit{Z}\right)}$	4.97 m	3.45 m
Twist ${\mathit{\Theta }}_{\mathit{t}\mathit{i}\mathit{p}}$	0.0°	4.19°

Table 9

Global mode shapes and the associated frequencies of the optimized metallic CRM wingbox (flutter constraint included).

Mode Number	Mode Shape Description	Natural Frequency [Hz](Full Fuel)	Natural Frequency [Hz](Zero Fuel)
1	1st Out-of-Plane Bending	1.07	1.63
2	2nd Out-of-Plane Bending	3.04	4.55
3	1st Lateral Bending	3.47	5.49
4	3rd Out-of-Plane Bending	6.83	8.83
5	1st Torsion	18.68	19.30

Table 10

Global mode shapes and the associated frequencies of the optimized composite CRM wingbox (flutter constraint included).

Mode Number	Mode Shape Description	Natural Frequency [Hz](Full Fuel)	Natural Frequency [Hz](Zero Fuel)
1	1st Out-of-Plane Bending	1.40	2.31
2	2nd Out-of-Plane Bending	3.75	4.98
3	1st Lateral Bending	3.99	7.40
4	3rd Out-of-Plane Bending	6.25	10.96
5	1st Torsion	23.09	24.32

References

1. Martins, J.R.R.A.; Lambe, A.B. Multidisciplinary Design Optimization: A Survey of Architectures. AIAA J.; 2013; 51, pp. 2049-2075. [DOI: https://dx.doi.org/10.2514/1.J051895]

2. Bindolino, G.; Ghiringhelli, G.; Ricci, S.; Terraneo, M. Multilevel Structural Optimization for Preliminary Wing-Box Weight. J. Aircr.; 2010; 47, pp. 475-489. [DOI: https://dx.doi.org/10.2514/1.41552]

3. Hurlimann, F.; Kelm, R.; Dugas, M.; Oltmann, K.; Kress, G. Mass estimation of transport aircraft wingbox structures with a CAD/CAE-based multidisciplinary process. Aerosp. Sci. Technol.; 2010; 15, pp. 323-333. [DOI: https://dx.doi.org/10.1016/j.ast.2010.08.005]

4. Gray, J.S.; Hwang, J.T.; Martins, J.R.R.A.; Moore, K.T.; Naylor, B.A. OpenMDAO: An open-source framework for multidisciplinary design, analysis, and optimization. Struct. Multidiscip. Optim.; 2019; 59, pp. 1075-1104. [DOI: https://dx.doi.org/10.1007/s00158-019-02211-z]

5. Hwang, J.T.; Ning, A. Large-scale multidisciplinary optimization of an electric aircraft for on-demand mobility. Proceedings of the 2018 AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, AIAA, AIAA; Kissimmee, FL, USA, 8–12 January 2018; [DOI: https://dx.doi.org/10.2514/6.2018-1384]

6. Hwang, J.T.; Lee, D.Y.; Cutler, J.W.; Martins, J.R.R.A. Large-scale multidisciplinary optimization of a small satellite’s design and operation. J. Spacecr. Rocket.; 2014; 51, pp. 1648-1663. [DOI: https://dx.doi.org/10.2514/1.A32751]

7. Gray, A.C.; Riso, C.; Jonsson, E.; Martins, J.R.R.A.; Cesnik, C.E.S. High-Fidelity Aerostructural Optimization with a Geometrically Nonlinear Flutter Constraint. AIAA J.; 2023; 61, pp. 2430-2443. [DOI: https://dx.doi.org/10.2514/1.J062127]

8. Dababneh, O.; Kipouros, T.; Whidborne, J. Application of an Efficient Gradient-Based Optimization Strategy for Aircraft Wing Structures. Aerospace; 2018; 5, 3. [DOI: https://dx.doi.org/10.3390/aerospace5010003]

9. Klimmek, T. Parametric Set-Up of a Structural Model for FERMAT Configuration for Aeroelastic and Loads Analysis. ASD J.; 2014; 3, pp. 31-49.

10. Jutte, C.V.; Stanford, B.K.; Wieseman, C.D. Internal Structural Design of the Common Research Model Wing Box for Aeroelastic Tailoring; NASA/TM–2015-218697 NASA: Hampton, VA, USA, 2015.

11. Vassberg, J.C.; DeHaan, M.A.; Rivers, S.M.; Wahls, R.A. Development of a common research model for applied CFD validation studies. Proceedings of the 26th AIAA Applied Aerodynamics Conference; Honolulu, HI, USA, 18–21 August 2008; AIAA 2008-6919.

12. Vassberg, J. A unified baseline grid about the common research model wing/body for the fifth AIAA CFD drag prediction workshop (invited). Proceedings of the 29th AIAA Applied Aerodynamics Conference; Honolulu, HI, USA, 27–30 June 2011.

13. Kenway, G.K.W.; Martins, J.R.R.A. AIAA ADODG Case 5: CRM Wingbody-Tail Optimization at Flight Reynolds Number; Technical Report AIAA: Reston, VI, USA, 2015.

14. Kenway, G.K.W.; Martins, J.R.R.A. Multipoint aerodynamic shape optimization investigations of the Common Research Model wing. Proceedings of the AIAA Science and Technology Forum and Exposition (SciTech); Kissimmee, FL, USA, 5–9 January 2015.

15. Lyu, Z.; Kenway, G.K.; Martins, J.R.R.A. Aerodynamic shape optimization investigations of the Common Research Model wing benchmark. AIAA J.; 2015; 53, pp. 968-985. [DOI: https://dx.doi.org/10.2514/1.J053318]

16. Kenway, G.; Kennedy, G.; Martins, J.R.R.A. Aerostructural optimization of the Common Research Model configuration. Proceedings of the 15th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference; Atlanta, GA, USA, 16–20 June 2014; [DOI: https://dx.doi.org/10.2514/6.2014-3274]

17. Kenway, G.K.; Martins, J.R. Multi-Point High-Fidelity Aerostructural Optimization of a Transport Aircraft Configuration. J. Aircr.; 2014; 51, pp. 144-160. [DOI: https://dx.doi.org/10.2514/1.C032150]

18. Handojo, V. Investigation of load alleviation in aircraft pre-design and its influence on structural mass and fatigue. Aerosp. Sci. Technol.; 2022; 122, 107405. [DOI: https://dx.doi.org/10.1016/j.ast.2022.107405]

19. Li, K.; Yang, C.; Wang, X.; Wan, Z.; Li, C. Multiscale Aeroelastic Optimization Method for Wing Structure and Material. Aerospace; 2023; 10, 866. [DOI: https://dx.doi.org/10.3390/aerospace10100866]

20. Baker, M.L.; Mendoza, R.; Hartwich, P.M. Transonic aeroelastic analysis of a high speed transport wind tunnel model. Proceedings of the AIAA Conference; St. Louis, MO, USA, 12–15 April 1999; AIAA-99-1217.

21. Riso, C.; Cesnik, C.E.S. Investigation of geometrically nonlinear effects in the aeroelastic behavior of a very flexible wing. Proceedings of the 2023 AIAA SciTech Forum, Online; National Harbor, MD, USA, 23–27 January 2023; [DOI: https://dx.doi.org/10.2514/6.2023-0759]

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

23. Wang, Z.; Wan, Z.; Groh, R.M.J.; Wang, X. Aeroelastic and local buckling optimisation of a variable-angle-tow composite wing-box structure. Compos. Struct.; 2021; 258, 113201. [DOI: https://dx.doi.org/10.1016/j.compstruct.2020.113201]

24. Dillinger, J.K.S.; Abdalla, M.M.; Klimmek, T.; Gürdal, Z. Static aeroelastic stiffness optimization and investigation of forward swept composite wings. Proceedings of the World Congress on Structural and Multidisciplinary Optimization; Orlando, FL, USA, 19–24 May 2013.

25. Handojo, V.; Himisch, J.; Bramsiepe, K.; Krüger, W.R.; Tichy, L. Potential Estimation of Load Alleviation and Future Technologies in Reducing Aircraft Structural Mass. Aerospace; 2022; 9, 412. [DOI: https://dx.doi.org/10.3390/aerospace9080412]

26. Brooks, T.R.; Kenway, G.K.W.; Martins, J.R.R.A. Benchmark Aerostructural Models for the Study of Transonic Aircraft Wings. AIAA J.; 2018; 56, pp. 2840-2855. [DOI: https://dx.doi.org/10.2514/1.J056603]

27. Chauhan, S.S.; Martins, J.R.R.A. Low-Fidelity Aerostructural Optimization of Aircraft Wings with a Simplified Wingbox Model Using OpenAeroStruct. Proceedings of the 6th International Conference on Engineering Optimization, EngOpt 2018; Lisbon, Portugal, 17–19 September 2018; Springer: Cham, Switzerland, 2018; pp. 418-431.

28. Brooks, T.R.; Kenway, G.K.W.; Martins, J.R.R.A. Undeflected common research model (uCRM): An aerostructural model for the study of high aspect ratio transport aircraft wings. Proceedings of the 18th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference; Denver, CO, USA, 5–9 June 2017; [DOI: https://dx.doi.org/10.2514/6.2017-4456]

29. Liu, Q.; Jrad, M.; Mulani, S.B.; Kapania, R.K. Global/Local Optimization of Aircraft Wing Using Parallel Processing. AIAA J.; 2016; 54, pp. 3338-3348. [DOI: https://dx.doi.org/10.2514/1.J054499]

30. Kennedy, G.J.; Martins, J.R.R.A. A Comparison of Metallic and Composite Aircraft Wings Using Aerostructural Design Optimization. Proceedings of the 12th AIAA Aviation Technology, Integration, and Operations (ATIO) Conference and 14th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference; Indianapolis, Indiana, 17–19 September 2012; [DOI: https://dx.doi.org/10.2514/6.2012-5475]

31. Kenway, G.K.; Kennedy, G.J.; Martins, J.R.R.A. A CAD-Free Approach to High-Fidelity Aerostructural Optimization. Proceedings of the 13th AIAA/ISSMO Multidisciplinary Analysis Optimization Conference; Fort Worth, TX, USA, 13–15 September 2010; [DOI: https://dx.doi.org/10.2514/6.2010-9231]

32. Dababneh, O.; Kipouros, T. Influence of high fidelity structural models on the predicted mass of aircraft wing using design optimization. Aerosp. Sci. Technol.; 2018; 79, pp. 164-173. [DOI: https://dx.doi.org/10.1016/j.ast.2018.05.043]

33. Dababneh, O. Multidisciplinary Design Optimisation for Aircraft Wing Mass Estimation. Ph.D. Thesis; Cranfield University: Bedford, UK, 2016.

34. FAA. FAR 25, Airworthiness Standards: Transport Category Airplanes (Title 14 CFR Part 25). Available online: http://flightsimaviation.com/data/FARS/part_25.html (accessed on 15 March 2025).

35. ASM. ASM Aerospace Specification Metals, Aluminum 2024-T3. 1978; Available online: http://asm.matweb.com/search/SpecificMaterial.asp?bassnum=%20MA2024T3 (accessed on 10 March 2025).

36. ASM. ASM Aerospace Specification Metals, Aluminum 7050-T7451. 1978; Available online: http://asm.matweb.com/search/SpecificMaterial.asp?bassnum=MA7050T745 (accessed on 12 March 2025).

37. Soni, S.R. Elastic Properties of T300/5208 Bidirectional Symmetric Laminates–Technical Report Afwal-Tr-80-4111; Materials Laboratory–Air Force Wright Aeronautical Laboratories–Air Force Systems Command: Dayton, OH, USA, 1980.

38. Torenbeek, E. Development and Application of a Comprehensive Design Sensitive Weight Prediction Method for Wing Structures of Transport Category Aircraft; Delft University of Technology: Delft, The Netherlands, 1992.

39. ESDU. Computer Program for Estimation of Spanwise Loading of Wings with Camber and Twist in Subsonic Attached Flow. Lifting-Surface Theory. 1999; Available online: https://www.esdu.com/cgi-bin/ps.pl?t=doc&p=esdu_95010c (accessed on 12 March 2014).

40. EASA. Certification Specifications and Acceptable Means of Compliance for Large Aeroplanes CS-25, Amendment 16. 2015; Available online: https://www.easa.europa.eu/en/document-library/certification-specifications/cs-25-amendment-16 (accessed on 21 March 2025).

41. Rodden, W.P.; Johnson, E.H. MSC Nastran Version 68 Aeroelastic Analysis User’s Guide; MSC Software Corporation: Newport Beach, CA, USA, 1994.

42. MSC. Flight Loads and Dynamics User’s Guide Version 2006; MSC Software Corporation: Los Angeles, CA, USA, 2006.

43. Jones, R.M. Mechanics of Composite Materials; 2nd ed. Taylor & Francis: Abingdon, UK, 1999.

44. Tsai, S.W.; Hahn, H.T. Introduction to Composite Materials; Technomic Publishing Co.: Chicago, IL, USA, 1980.

45. Kassapoglou, C. Review of Laminate Strength and Failure Criteria, in Design and Analysis of Composite Structures: With Applications to Aerospace Structures; John Wiley & Sons Ltd.: Oxford, UK, 2013.

46. Starnes, J.R., Jr.; Haftka, R.T. Preliminary Design of Composite Wings for Buckling, Stress and Displacement Constraints. J. Aircr.; 1979; 16, pp. 564-570. [DOI: https://dx.doi.org/10.2514/3.58565]

47. Oliver, M.; Climent, H.; Rosich, F. Non Linear Effects of Applied Loads and Large Deformations on Aircraft Normal Modes. Proceedings of the RTO AVT Specialists’ Meeting on Structural Aspects of Flexible Aircraft Control; Ottawa, ON, Canada, 18–20 October 1999.

48. Liu, Q.; Mulani, S.; Kapani, R.K. Global/Local Multidisciplinary Design Optimization of Subsonic Wing, AIAA 2014-0471. Proceedings of the 10th AIAA Multidisciplinary Design Optimization Conference–AIAA SciTech; National Harbor, MD, USA, 13–17 January 2014; AIAA, Inc.: Reston, VA, USA, 2014.

49. Barker, D.K.; Johnson, J.C.; Johnson, E.H.; Layfield, D.P. Integration of External Design Criteria with MSC Nastran Structural Analysis and Optimization. Proceedings of the Worldwide Aerospace Conference & Technology Showcase; Toulouse, France, 8–10 April 2002; MSC Software Corporation: Newport Beach, CA, USA, 2002.

50. Hajela, P. A Root Locus-Based Flutter Synthesis Procedure. J. Aircr.; 1983; 20, pp. 1021-1027. [DOI: https://dx.doi.org/10.2514/3.48206]

51. MSC Nastran 2012 Design Sensitivity and Optimization User’s Guide; MSC Software Corporation: Santa Ana, CA, USA, 2012.

52. Neill, D.J.; Johnson, E.H.; Canfield, R. ASTROS—A Multidisciplinary Automated Structural Design Tool. J. Aircr.; 1990; 27, pp. 1021-1027. [DOI: https://dx.doi.org/10.2514/3.45976]

53. McLean, M.D. MSC Nastran 2012 Linear Static Analysis User’s Guide; MSC Software Corporation: Santa Ana, CA, USA, 2012.

54. Herrmann, L.R. Improved Stress Calculations for Simple Quadrilateral Elements. Comput. Struct.; 1976; 6, pp. 141-148. [DOI: https://dx.doi.org/10.1016/0045-7949(76)90064-X]