Flutter and post-flutter constraints in aircraft design optimization

Abstract

Flutter is a dynamic aeroelastic instability driven by the interaction of inertial, elastic, and aerodynamic forces. It is an undesirable phenomenon in aircraft because it causes divergent oscillations that may lead to structural damage or failure, performance and ride comfort degradation, or loss of control. If flutter is discovered at the aircraft certification stage, costly redesign is required. Performing flutter analysis early in the design process can mitigate this problem. Furthermore, including flutter analysis as a constraint in multidisciplinary design optimization reduces the risk of costly modifications late in the design cycle. We review the methods for flutter analysis in the context of aircraft design optimization. We also include methods for predicting post-flutter limit cycle oscillations due to the increasing impact of nonlinear effects on future aircraft. While there has been extensive work in flutter and post-flutter analyses, developing design optimization constraints associated with these analyses has additional requirements, such as acceptable computational cost, function smoothness, robustness, and derivative computation. We discuss these requirements and review efforts in the development, implementation, and application of flutter and post-flutter constraints in aircraft design optimization. We conclude the paper by summarizing the current state of this field and the main open problems. Flutter constraints have been included in structural optimizations, but optimizing both the structural sizing and the aerodynamic shape remains a challenge due to the need to recompute the aerodynamic properties at each design iteration. Additional difficulties arise in the presence of large structural deflections and transonic flow conditions due to the dependency of the flutter point on the equilibrium state and the high cost of nonlinear computations. Post-flutter constraints have rarely been included into design optimization, but they are crucial in the prevention of undesirable limit cycle oscillations. Implementing such constraints requires the development of more efficient and robust prediction methods that can handle realistic configurations. While this paper focuses on flutter and post-flutter constraints for aircraft design optimization applications, the considerations and challenges are broadly applicable to the optimization of engineering systems including stability and post-critical dynamic constraints.

Introduction

Flutter is a dynamic aeroelastic instability that causes divergent oscillatory vibrations [1]. It is an undesirable phenomenon in aircraft because it can cause structural damage or failure, performance and ride comfort degradation, or loss of control. Flutter computations are typically performed only after an initial detailed design of the aircraft is completed, because they require the vehicle stiffness, mass, and aerodynamic models to be available [2]. If the design does not satisfy the flutter requirements at this stage, a redesign is necessary, which adds costs and delays to the aircraft development cycle. Thus, it is desirable to consider flutter concurrently with the aircraft design and the wing design in particular. Such a process would not only shorten the design cycle, but also allow for advantageous design trade-offs between the flutter requirements, the other constraints, and the aircraft performance.

Performing multidisciplinary design optimization (MDO) that considers both aerodynamic shape and structural sizing simultaneously while enforcing flutter constraints is a way to address this issue [3,4]. Structural optimization alone, even if including aerostructural analyses for enforcing flutter constraints, yields design solutions with suboptimal performance compared to the optimal designs resulting from MDO, where structural and aerodynamic sizing variables are optimized simultaneously [4,5]. MDO can minimize structural weight, fuel consumption, or a combination of these two objectives with respect to wing shape, internal structure arrangements, and sizing, while accounting for the interactions between aerodynamics, structures, and other disciplines, and satisfying various constraints. MDO with flutter constraints results in designs with optimal aeroelastic tailoring. Omitting flutter constraints in the MDO process when minimizing fuel consumption tends to yield light-weight, high-aspect-ratio wing (HARW) designs that despite being highly efficient may not be feasible [6,7].

After the aircraft has been designed and a prototype has been built, certification requires flight tests to demonstrate that the aircraft be free from flutter in the flight envelope with a 15% safety margin beyond the dive speed. If flutter is discovered at the flight test certification stage, it requires redesign to address it, incurring additional costs. The redesign effort typically increases the structural weight, reducing the performance originally anticipated for the aircraft.

The trend towards HARW aircraft is driven by better fuel efficiency, but their increased flexibility makes it all the more important to consider flutter accurately and early in the design process [8]. Another recent trend is the increasing use of control surfaces to suppress flutter. Active flutter suppression systems can be incorporated late in the design process when aeroelastic instabilities are encountered and a passive solution such as redesign is impractical and expensive [9,10]. Alternatively, MDO provides a way to obtain the best possible configuration by co-designing the wing shape and internal structure, which contribute to passive flutter suppression that can then be augmented with an active flutter control system.

While there has been extensive work in methods for flutter analysis, integrating flutter constraints into design optimization requires additional considerations. Models used for flutter prediction should capture the relevant physics with adequate accuracy to correctly drive the optimizer and inevitably there is a compromise between model fidelity and computational cost. To include flutter analysis in a numerical optimization cycle, speed of execution is particularly important to make sure that the overall optimization process is tractable.

Another important characteristic for integrating flutter analysis into the optimization process is the robustness of the flutter prediction method. Since the optimization process automatically samples the design space, it is likely to request for the analysis of designs that would normally not be chosen by a human designer. Thus, it is important that the flutter analysis converges for such designs so that the overall optimization process is not interrupted.

Gradient-based optimization algorithms are needed to optimize practical aircraft configurations parameterized with a large number of design variables [11]. When using gradient-based algorithms, it is important to consider the smoothness of the objective and constraint functions, as well as the accuracy and efficiency of the derivative computations.

Optimizing HARW configurations subject to flutter constraints is even more challenging because it requires capturing couplings between aeroelasticity and flight dynamics along with geometric nonlinearities that arise in the presence of low natural vibration frequencies and significant structural flexibility [12,13].

Nonlinearities in the structure (large deflections, free-play of control surfaces, follower loading) or the aerodynamics (shock waves and flow separation) can cause self-sustained oscillations of limited amplitude that remain constant in time, known as limit cycle oscillations (LCOs). For certain types of nonlinearities, LCOs may exist at flight conditions below the flutter point [14]. When nonlinear effects become important, post-flutter analysis should be integrated into the design process in the form of constraints to make sure that the optimal design be feasible. In the context of aircraft design practice, LCO typically refers to oscillations due to control-surface free-play nonlinearities. From a certification perspective, these are typically addressed by computing free-play tolerances. In this paper, we use the term LCO to refer more generally to a post-flutter response due to structural or aerodynamic nonlinearities. While free-play LCOs are a particular case of post-flutter response, the paper focuses on LCOs involving vehicle components or the entire airframe. Taking into account these global LCOs is challenging in design optimization due to the high cost of nonlinear post-flutter computations and the difficulties in formulating post-flutter constraints.

There have been several review papers and textbooks on flutter and post-flutter analysis. Livne [15,16] reviewed the state-of-the-art and future challenges in aeroelasticity of conventional and unconventional vehicles. A recent review by the same author focused on active flutter suppression control systems [10]. Friedmann [17] reviewed the general challenges in nonlinear aeroelasticity, where the applications focused on rotary wings. Dowell et al. [14] classified nonlinear aeroelastic behaviors and discussed theoretical, computational, and experimental analysis efforts. de C. Henshaw et al. [18] discussed traditional industrial linear flutter prediction and recent efforts for including nonlinear effects, particularly due to transonic flows. More recently, Afonso et al. [8] reviewed nonlinear aeroelasticity of HARWs. Dimitriadis [19] discussed nonlinear post-flutter behaviors in aeroelastic systems and the related analysis methods. Beran et al. [20] reviewed methods for uncertainty quantification in aircraft aeroelasticity and their application to formulate nondeterministic optimization problems. However, the field lacks a review on the integration of flutter and post-flutter analysis as constraints in aircraft design optimization.

In this paper, we address this shortcoming by reviewing methods for flutter and post-flutter prediction, and we discuss their advantages and disadvantages in the context of aircraft design optimization. First, we provide a brief background on multidisciplinary design optimization (Section 2) and on flutter and post-flutter modeling (Section 3). These sections emphasize the key aspects relevant to aircraft design and are beneficial for readers not familiar with either of these topics. Next, we present numerical methods for flutter and post-flutter analysis, which are summarized in Fig. 1. These methods and examples of their application in aircraft design optimization problems are discussed in Section 4 and 5. In reviewing previous work we focus on deterministic optimizations, but we also present few examples where constraints are formulated using failure probabilities (reliability-based design optimization) [20]. The paper concludes with remarks on the state of this field and the open challenges to be addressed for integrating flutter and post-flutter considerations into the optimal design of aircraft configurations. While the paper focuses on aircraft design optimization, the methods and challenges reviewed are more broadly applicable to the design optimization of engineering systems subject to stability and post-critical dynamic constraints.

Multidisciplinary design optimization couples the relevant disciplines of an engineering system and performing numerical optimization to aid the design of that system [21]. MDO considers several disciplines simultaneously such that their interactions can be leveraged, resulting in a better optimum than if each discipline were optimized sequentially [22]. Thus, considering MDO early in the design process allows engineers not only to improve the design, but also to minimize development time and cost of the overall design.

Performing MDO of aircraft configurations by describing its outer mold line (OML) and structural sizing with high fidelity requires a large number of design variables. Detailed aerodynamic optimization of wings requires hundreds of shape variables [23] and structural sizing based on a detailed finite-element wingbox model that is best utilized with an equally large number of sizing variables [24]. Gradient-based optimization methods are the feasible way to solve for high-dimensional problems within a reasonable computational time, especially when using high-fidelity analyses [11,25]. Gradient-based methods require derivatives of the objective and constraint functions with respect to the design variables to help the optimization algorithm find the most promising search directions and establish rigorous optimality conditions.

While gradient-free algorithms are typically more robust and some of them explore the design space more widely, their cost is prohibitive when the number of design variables is large. Although gradient-based methods only guarantee convergence to a local optimum, this can be mitigated by using a multi-start technique [26]. Furthermore, recent studies failed to find multiple local minima (multimodality) in airfoil and wing shape design optimization [11,23], except the case of planform optimizations where expected local minima were found related to choices such as upwards or downwards winglets [27].

The efficacy of gradient-based algorithms relies on accurate and efficient gradient computations. The gradient accuracy directly affects the ability to converge to the optimum with a specified tolerance and the order of convergence of the optimization. In the best case, inaccurate gradients increase the number of iterations required for convergence and in the worst case cause early stopping due to numerical issues. Efficiency gradient computation is also important because this computation is sometimes the bottleneck in the optimization cycle.

When it comes to methods for computing gradients, the finite-difference method is a popular choice because it is easy to implement and can always be used, even with black-box codes. The major drawbacks of the finite-difference method is that it is inaccurate and its computational cost scales poorly with the number of design variables. Unlike the finite-difference method, the complex-step method [28] is accurate, but its cost still scales unfavorably with the number of design variables, making it prohibitive for wing design applications. Furthermore, the complex-step method cannot be directly applied to programs that use complex numbers. Such programs need to be made complex-step compatible by modifying them so that the real and imaginary part of the complex number are represented as two real numbers. Another option for computing gradients is automatic differentiation (AD), which uses a software tool to parse the code of an analysis to produce a new code that computes derivatives of that analysis [29,30]. Although AD can scale well with the number of variables, it does not handle iterative simulations efficiently in general. Finally, analytic methods are the most desirable because they are both accurate and efficient, especially for iterative simulations [31]. However, they require significant implementation effort. There are two main approaches within the analytic methods: the direct approach and the adjoint approach. The adjoint approach is especially attractive because the computational cost is dependent of the number of outputs of interest (objectives and constraints) but independent of the number of design variables [31,32]. A coupled-adjoint approach can solve static aeroelastic problems [5,6] and can be generalized to any multidisciplinary problem [33,34].

In the context of aircraft design optimization subject to flutter or post-flutter constraints, most of the early efforts used gradient-based optimization with gradients computed with either finite differences or the direct analytic method. However, recent efforts implemented the more efficient adjoint approach, and some also used AD techniques. These applications are further discussed in Section 4.2 for flutter constraints and Section 5.3 for post-flutter constraints.

For aircraft designs to be useful and practical, the underlying models used in the flutter and post-flutter analysis need to capture the correct physics involved. However, a simplification of the phenomena is often necessary to make problems tractable to solve. Therefore, the choice of model should balance the fidelity needed to obtain accurate predictions and the mathematical or computational tractability for design applications.

This section highlights the modeling aspects to be considered in flutter and post-flutter analysis, which are discussed in more detail in Section 4 and 5 for obtaining meaningful results in a design optimization. By flutter, we mean the onset of divergent oscillations as the flight conditions of aircraft cross the critical stability boundary (flutter boundary).

Mathematically, flutter occurs at a Hopf bifurcation point [35] beyond which the system is in the post-flutter regime. Several post-flutter behaviors are possible, as discussed in detail by Dimitriadis [19] for a two-dimensional aeroelastic system with stiffness and damping nonlinearities. Among these behaviors, we are particularly interested in self-sustained oscillations with limited amplitude that remains constant in time, known as LCOs.

LCOs typically develop beyond the flutter boundary; however, for certain types of nonlinearities, they can also occur before reaching the flutter boundary [12]. Integrating post-flutter analyses into the design process can prevent this undesirable situation. This places additional challenges due to the difficulties in modeling nonlinear aeroelastic systems accurately and the inherent high cost of nonlinear aeroelastic computations.

Flutter is defined as a self-exciting dynamic instability that is associated with the interaction of inertial, elastic, and aerodynamic forces [1]. At the onset of flutter, this aeroelastic instability can be physically described as an oscillation with a small amplitude that is constant in time triggered by a small-amplitude disturbance, as shown in Fig. 2. The flight condition in which the system damping vanishes, resulting in this self-sustained oscillation, represents the flutter point (or flutter boundary). For linear systems, the flutter point is defined as the minimum dynamic pressure at which at least one of the modes becomes unstable [1]. The dynamic pressure can be replaced by equivalent airspeed, and is a function of altitude and Mach number.

Past the flutter point, in the absence of restraining nonlinearities from the aerodynamics, the structure, or both, the amplitude of the oscillations grow exponentially. Fluid-structure interactions may also result in a static instability called divergence [1], which is not associated with oscillations. As for flutter, the structural response grows unbounded past the onset point, eventually reaching a limited-amplitude oscillation if restraining nonlinearities are present.

In the following discussion, we focus mainly on flutter phenomena, because for many practical configurations flutter occurs before divergence. However, accounting for divergence and the associated post-critical response in the design process shares many of the modeling and analysis aspects associated with flutter. Furthermore, some of the analysis methods and constraints discussed in Section 4 are applicable to divergence as well as flutter. Moreover, in this paper we focus on global wing or component flutter rather than localized effects such as panel flutter that typically occurs at supersonic flow conditions.

The possible flutter characteristics are illustrated in Fig. 3, which shows the variation of the modal damping with flight speed at a fixed altitude. This is known as $V-g$ diagram, which is a classical tool used in linear flutter analyses for determining the flutter point and interpreting the flutter characteristics. A similar representation can be obtained by varying dynamic pressure at a fixed Mach number.

Damping changes with flight speed in different ways among different designs, leading to different flutter behaviors. Soft flutter occurs when the damping decreases gradually with increasing flight speed, while hard flutter occurs when this decrease is abrupt. Another possibility is that there is a gradual decrease in damping with increasing flight speed, all the way to cross the zero value and beyond, followed by a damping increase. This phenomenon is known as a hump mode. These concepts are important when considering how to formulate a smooth and continuous flutter constraint for gradient-based optimization and are discussed in Section 4.

In flutter analysis, the physics described above is often represented by less expensive linear models due to the number of conditions to be considered for certification. However, nonlinear structural and aerodynamic effects or the interaction between elastic and rigid-body degrees of freedom (DOF), which become important in the presence of low structural vibration frequencies, may significantly impact the flutter point. Therefore, the flutter prediction accuracy depends on the appropriate modeling of nonlinear effects and boundary conditions.

Furthermore, nonlinear effects impact not only the models used in flutter analysis, but also the analysis methods themselves. For linear systems, flutter characteristics do not depend on the deformation state. Therefore, flutter is typically analyzed by considering the unloaded and undeformed structure. For nonlinear systems, stability characteristics vary with the deformation configuration. Therefore, flutter analysis must be performed by computing the eigenvalues of the aeroelastic system linearized around equilibrium points for each flight condition to identify at what point the damping vanishes [12] as shown in Fig. 4. The eigenvalues can be computed by considering both the elastic and rigid-body DOFs (flutter in free flight) or by retaining only the elastic DOFs (traditional flutter) or the rigid-body DOFs (flight dynamic stability).

Many possible sources of nonlinearities can be present simultaneously in aeroelastic systems [8,14,19]. In this paper, we focus on aerodynamic nonlinearities due to transonic flow regimes and geometric structural nonlinearities due to large deflections, both of which are critical in the design of next-generation transport aircraft.

Aerodynamic nonlinearities due to shock waves and flow separation significantly impact the flutter speed. This decreases dramatically in the transonic regime, a phenomenon known as the transonic dip [[36], [37], [38], [39], [40]] illustrated in Fig. 5. Low-order, linear unsteady aerodynamic models commonly used in flutter analysis are in general accurate enough for subsonic and supersonic flows, but they severely overestimate the flutter speed for transonic conditions [[41], [42], [43]].

As shown in Fig. 5 for a hypothetical wing, linear theory is non-conservative when compared to nonlinear viscous models. Nonlinear inviscid models based on Euler equations can capture shock waves but they still fail to accurately predict the flutter boundary [44]. In many cases, the nonlinear inviscid theory predicts a highly conservative flutter speed at the dip, even though it is generally closer to viscous theory predictions. Depending on the severity of shocks, models based on Navier–Stokes equations (which include viscous and turbulence effects like boundary layer thickening, flow separation, and interactions between shocks and regions of separated flow) are necessary to obtain accurate flutter points [42]. Studies on various geometries demonstrated that taking into account viscous phenomena in the transonic flow regime improves the numerical prediction of transonic dip [[45], [46], [47], [48]].

A common approach to improve the accuracy of transonic flutter computations while minimizing the increase in computational cost is to use numerical or experimental corrections applied to potential-flow linear models [49]. However, the correction data may not be available for optimization, either because it requires high-fidelity computations that are too expensive or because it is obtained from wind-tunnel measurements. This problem can be addressed by analyzing flutter using time-accurate dynamic simulations and higher-fidelity aerodynamic models. On the other hand, flutter prediction based on time-accurate computational fluid dynamics (CFD) is a challenge even for just analyzing the final configuration and is currently prohibitive for design space exploration.

Methods exist that try to preserve the computational efficiency of lower-fidelity methods while retaining the nonlinear physics modeled by the higher-fidelity CFD methods. One possibility is to use time-linearized transonic small disturbance (TSD) equations. Linear small-disturbance theory is inadequate for capturing strong transonic shocks, but small-disturbance solutions about the steady nonlinear background flow computed using high-fidelity CFD can provide acceptable performance and accuracy [50]. Another possibility is to use the transonic equivalent strip (TES) theory [51] and a provided pressure distribution from either experimental data or a high-fidelity CFD code to compute the small-disturbance transonic aerodynamic loads for flutter analysis [52]. Furthermore, several efforts have applied the time-linearization directly to CFD solvers. This approach is appealing since they retain nonlinear effects, provide accurate results in the transonic regime, and intrinsically account for geometric properties of the body such as thickness and camber.

Kreiselmaier and Laschka [53] developed an unsteady method based on the small-disturbance Euler equations, which was later extended to small disturbance Navier–Stokes equations to include viscous effects [54]. The proposed method produced good results in the transonic flow regime [55,56]. Thormann and Winghalm [57] developed a linear frequency domain (LFD) solver taking advantage of preconditioned Krylov GMRES [58] solution method. Later, Widhalm and Thormann [59] improved the algorithm and provided the analytic derivatives needed in the solution, improving the solver efficiency. The method was shown to be accurate when compared to full unsteady time-marching solutions.

These approaches consider unsteady aerodynamic models linearized about nonlinear equilibrium states and thus can capture the impact of static nonlinear effects on flutter. Moreover, they demonstrate computational savings well beyond an order of magnitude compared to fully unsteady time-marching solutions [53,57]. However, computing derivatives of such methods for optimization is challenging due to the need for second-order derivative information.

Motivated by the interest in capturing key transonic flow physics with low computational cost, recent efforts also developed low-order unsteady transonic aerodynamic models suitable for integrating transonic flutter analyses into aircraft design.

Skujins and Cesnik [60] proposed a reduced-order unsteady aerodynamic model for multiple Mach regimes based on linear convolution with a nonlinear static correction. The methodology included error estimation capabilities based on the newly developed method of segments, which represents a flexible wing as a collection of rigid spanwise segments subject to local angle of attack and Mach number conditions. The method of segments was also applied to transonic flutter analysis of a transport vehicle by Kitson and Cesnik [61].

Mallik et al. [62] developed a reduced-order model for HARW configurations by combining time-linearized Leishman indicial functions [63] with lift-curve and moment-curve slopes obtained by solving the RANS equations around airfoils at various Mach number and angle of attack conditions and for various thickness ratios. They obtained a state-space formulation for the airfoil unsteady aerodynamics to be used for eigenvalue analysis, which was extended to three-dimensional HARW discretized in spanwise strips by accounting for sweep correction. Flutter results were compared with wind-tunnel experiments for a truss-braced wing (TBW) configuration. The low-order model captured the transonic dip that was not predicted by potential-flow theory and presented good agreement with experiments at significantly lower computational cost compared to unsteady RANS simulations. These results showed the method suitability for conceptual HARW aircraft design including transonic flutter constraints.

Opgenoord et al. [64] developed a physics-based two-dimensional low-order model for transonic airfoils using the perturbations of the lowest-order volume-source and vorticity moments with respect to a known nonlinear background flow solution as the states. Evolution equations for these perturbations were derived and calibrated using data from high-fidelity Euler CFD simulations. A state-space unsteady aerodynamic model was obtained for airfoil flutter analysis which was later extended to three-dimensional HARW configurations [65] using strip theory and sweep correction, as done by Mallik et al. [62]. The method was applied in conceptual design and optimization problems including transonic flutter considerations [65,66].

In addition to capturing transonic effects, a more recent flutter modeling need is to take into account geometric structural nonlinearities. These are particularly important in the analysis of HARW configurations, which achieve higher energy efficiency at the cost of increased structural flexibility and thus experience large deflections under normal operating load conditions. The changes in geometric shape and stiffness properties due to these deflections significantly affect the flutter boundary [13]. When deformations are large, traditional flutter analysis based on the vehicle undeformed shape does not capture the actual behavior of the aircraft during flight.

Studies on isolated HARWs [67], high-altitude long-endurance (HALE) flying-wing configurations [12,68,69], and commercial transport vehicles [61] pointed out the need to analyze very flexible aircraft in statically deformed configuration at trim, which varies with the flight condition. Including structural nonlinearities in flutter prediction is challenging for both analysis and design due to the high computational cost of nonlinear aeroelastic simulations and the flutter boundary dependency on the deformation state, which is not considered in linear approaches.

Finally, classical wing flutter analyses typically assume the vehicle to be clamped at the wing root. While this may be an acceptable simplification for some vehicles, it does not reflect the vehicle behavior in free flight [15,16]. For some configurations, simply including rigid-body DOFs influences the flutter solution substantially [12]. This occurs due to the coupling between rigid-body motion and structural dynamics that arise in the presence of low natural vibration frequencies. These interactions usually result in lower flutter points than the cantilevered configurations or different flutter mechanisms like body-freedom flutter (BFF) [70]. Therefore, it is imperative to understand the effect of boundary conditions and state variables on flutter prediction.

Mazidi et al. [71] investigated the effect of engine placement and roll maneuver on flutter results. They observed that the roll maneuver has a destabilizing effect on the flutter boundary dependent on the bank angle. Additionally, the location of the engine or external store greatly affects the flutter boundary and the roll-induced effects. Nearly all vehicles perform roll maneuvers during turns, making the inclusion of these conditions relevant to the aircraft design process.

There has been further work on the effect of rigid-body DOFs on the flutter problem. Niblett [70] investigated the causes of BFF for a conventional wing-fuselage configuration using a linear analytical flutter method. Su and Cesnik [12] investigated the flutter behavior of a blended wing-body (BWB) aircraft both for cantilevered and free-flying conditions using the University of Michigan's Nonlinear Aeroelastic Simulation Toolbox (UM/NAST). They observed a reduction in the flutter speed when the rigid body DOFs were included compared to the cantilever case. Moreover, the flutter mode changed to include pitch and plunge motions, resulting in BFF. Similarly, Jones and Cesnik [72] investigated the BFF characteristics of the X-56A experimental aircraft, describing the entire modeling process used for the flutter prediction. Cesnik and Su [73] analyzed the University of Michigan's X-HALE very flexible aeroelastic testbed [74] and observed that significant wing deformations can drive lateral BFF due to the coupling of the Dutch roll and asymmetric wing bending modes. Su and Cesnik [69] investigated the stability and dynamic response of a highly flexible flying wing for different payload configurations and gust disturbances. They found that wing deflections can lead to an unstable phugoid mode and an aperiodic short-period mode. Similar behaviors were observed by Patil and Hodges [68] and Patil and Taylor [75]. Richards et al. [76] also analyzed the coupled flight dynamics and aeroelasticity of flying wings. They noted that BFF occurred due to a coupling of the short period pitching mode and the first elastic bending mode. They compared different inertial configurations of the aircraft and noted that BFF depends largely on the inertia about the pitch axis. They found a boundary value for pitch inertia that uncoupled the pitch and bending modes, thus replacing BFF with a more conventional flutter. After parameter studies, they concluded that BFF was caused by low fuselage inertia, which could be mitigated by redistributing the fuselage mass.

BFF is not exclusive to flying-wing configurations; conventional tube-and-wing aircraft can also encounter this type of instability. Cavallaro et al. [77] investigated the flutter behavior of a Prandtl Plane (boxed wing aircraft) including rigid-body DOFs using MSC Nastran. Similarly to Richards et al. [76], they observed a dependence on the fuselage mass in causing BFF: while the baseline configuration encountered flutter without rigid-body contributions, increasing the fuselage mass resulted in BFF. Therefore, capturing these phenomena in flutter analyses is advisable even for conventional configurations and imperative when investigating non-conventional ones [15,16].

When aircraft flutter, disturbances perturbing an equilibrium state cause an oscillatory response with divergent amplitude. This may ultimately result in structural failure or, if the amplitude growth is gradual enough, the aircraft can be brought back at a stable operating point by pilot or control action. Dynamic nonlinearities in the aerodynamics or the structure, however, usually play a stabilizing role and restrain divergent motions, leading to a self-sustained oscillatory response of limited amplitude that remains constant in time, known as LCO. Unfortunately, nonlinear effects tend to be destabilizing and cause LCOs even before reaching the flutter point [14]. Furthermore, LCOs may or may not develop for a given flight condition depending on the amplitude of disturbances (initial conditions) and previous motion history (hysteresis) [78]. Therefore, when nonlinear effects become significant, post-flutter phenomena need to be considered in addition to the flutter boundary in the design process.

As mentioned in Section 1, in aircraft aeroelastic qualification, the term LCO is typically associated with control-surface free-play. In this paper, LCO has the more general meaning of post-flutter response leading to the development of self-sustained oscillations in vehicle components or in the entire structure. When nonlinear behavior is anticipated in aircraft configurations, optimization should include constraints on the post-flutter response, in addition to flutter constraints, to avoid the occurrence of undesirable LCOs at the operating conditions.

Post-flutter behaviors are typically classified using bifurcation diagrams. These diagrams represent the LCO amplitude variation with a chosen control parameter, called the bifurcation parameter. For aeroelastic systems, this is usually the flight speed or the dynamic pressure [79]. Two qualitative post-flutter scenarios are possible, as illustrated in Fig. 6.

Fig. 6 (a) shows the bifurcation diagram for a system with supercritical post-flutter behavior. When the system is perturbed before the flutter point, the response always converges towards the initial equilibrium state regardless of the magnitude of the disturbance. When the perturbation is applied beyond flutter, the response grows unbounded until dynamic stabilizing nonlinearities lead to an LCO with limited amplitude that keeps constant in time. Each point on the bifurcation diagram gives the LCO amplitude as function of the flight parameter. In the case of supercritical behavior, the LCO amplitude increases smoothly with the flight parameter and with a decreasing slope as stabilizing nonlinearities become stronger [14].

The bifurcation diagram of a system with subcritical post-flutter behavior is illustrated in Fig. 6(b). In this case, an unstable bifurcation branch (dashed line) exists below the flutter point. As a result, when disturbances are sufficiently large, the initial equilibrium configuration is lost. If higher-order stabilizing nonlinearities are present, both an unstable and a stable (solid line) solution branches exist for a range of flight conditions. In this case, the response grows until restraining nonlinearities lead to an LCO with an amplitude given by the upper stable bifurcation branch. Supercritical bifurcation is characterized by discontinuous behavior with respect to the amplitude of initial disturbance due to destabilizing nonlinear effects.

Supercritical LCOs are undesirable in aircraft, since they cause structural fatigue or failure, performance and ride comfort degradation, or loss of control. However, if the flutter boundary is satisfied with the required safety margins, LCOs never develop since the aircraft always operates below the flutter speed by certification. Additionally, enforcing flutter constraints during design can guarantee that flutter does not occur and therefore, supercritical LCOs are not encountered.

However, subcritical post-flutter behavior is a more serious problem, because instabilities may occur at flight conditions that are expected to be stable according to a conventional flutter analysis. Depending on the type of nonlinearities, these flight conditions could lie within the flight envelope. This considerably reduces the effectiveness of flutter constraints enforced during design.

When higher-order stabilizing nonlinearities are present, as shown in Fig. 6(b), the system suddenly responds with a large-amplitude LCO when perturbations before the flutter point pass the amplitude limit given by the unstable solution branch. Similarly, the system suddenly responds with a large amplitude LCO when it is perturbed beyond the flutter point. Therefore, subcritical bifurcations are characterized by discontinuous behavior with respect to both the initial conditions and the flight parameters.

In contrast, LCO amplitude grows gradually and with much smaller amplitudes past the flutter point in supercritical systems (see Fig. 6(a)). Furthermore, subcritical post-flutter response is characterized by hysteresis. Fig. 6(b) shows that once a stable LCO develops, the system does not recover the initial equilibrium state when the flight parameter is reduced below the critical (flutter) value.

To eliminate the LCO, it is necessary to reduce the flight parameter below the value at the folding point connecting the unstable and stable solution branches. Depending on the type of nonlinearities present in the system, this may be significantly lower than the flutter value. The above considerations make clear that when nonlinear effects become significant, design optimization should include post-flutter constraints for obtaining feasible optimal designs that do not experience subcritical LCOs.

The post-flutter behavior of wings with structural and aerodynamic dynamic nonlinearities was analyzed by various researchers [12,78,[80], [81], [82]]. This topic was also included in review papers [8,14,18]. However, only a few efforts integrated post-flutter analyses in design optimization as a constraint, as discussed in Section 5.4.

Because flutter is a safety-critical phenomenon, analyses and experimental investigations are required for vehicle certification. Analyzing a configuration for flutter late in the development cycle is likely to result in an inefficient design solution or pose challenges mitigating unexpected instabilities, resulting in performance decrease, financial losses, or both. For this reason, flutter should be integrated into the design process early in the form of a constraint.

In this section, we review flutter prediction methods and previous research that addressed optimization subject to flutter constraints. Past work has focused primarily on flutter analyses using linear structures and linear aerodynamic models. Further work considered transport aircraft operating in the transonic regime and thus used linear structures and nonlinear aerodynamics. The recent trend towards more flexible aircraft has led to research on constraining flutter for geometrically nonlinear structures.

Due to the sheer number of design variables typically used in practical aircraft optimization, most of the previous work reviewed in this paper used gradient-based methods. However, studies on geometrically nonlinear configurations were limited to simple structures parametrized by few design variables. For this reason, efforts that considered gradient-free optimizations are also included.

Flutter computations are typically performed in the frequency domain by solving an eigenvalue problem. Well-established eigenvalue-analysis methods exist for flutter analysis of linear aeroelastic systems (e.g., k-, p-, $pk$- and g-method) [[83], [84], [85], [86], [87]]. These methods are also applicable to nonlinear systems by linearizing the equations of motion about the nonlinear equilibrium configuration at each flight condition for capturing static nonlinearities due to the structure (large deflections) or the aerodynamics (background transonic flow). Direct methods, based on the Hopf-bifurcation theory, can also be used to predict the flutter point of nonlinear aeroelastic systems directly in the frequency domain [18].

It is also possible to predict flutter in the time domain, but this incurs a much higher computational cost. When analyzing the stability of a system in the time domain, the flutter point is typically evaluated by perturbing equilibrium configurations at different flight conditions and by time-marching the equations of motion in order to verify the decay or growing of the response, or to extract damping values in a post-processing stage [[88], [89], [90], [91], [92], [93], [94], [95]].

The computational cost of time-domain flutter analysis based on transient simulations is currently prohibitive for optimization, particularly in the presence of aerodynamic or structural nonlinearities. The high computational cost is due to the large number of computations required to evaluate the flutter speed by means of flight speed (or dynamic pressure) sweep or bisection. Additionally, ascertaining the stability close to the flutter point requires long integration times due to small damping values.

For gradient-based optimization, another challenge is the efficient computation of derivatives for time-marched systems. Adjoint methods are advantageous for optimizations with many design variables, but they are computationally expensive when applied to time marching solvers [96]. This is because the adjoint solution requires a reverse time integration following the forward integration for computing the system response [97], which results in large computational times and storage memory requirements [98]. Due to the above limitations, we focus on frequency-domain flutter prediction methods, but some fully coupled unsteady aeroelastic adjoint implementations do exist in literature [[99], [100], [101], [102], [103], [104]]. Eigenvalue-based methods are discussed in Section 4.1.1 and direct methods in Section 4.1.2.

The discrete equations of motion for a generic linear aeroelastic system can be written as

$$\mathbf{M}\ddot{\mathbf{u}}\left(t\right)+\mathbf{K}\mathbf{u}\left(t\right)-{\mathbf{F}}_{\text{aero}}\left(t\right)=0\text{,}$$

where $N_s$ is the number of structural DOFs (nodal displacements and rotations), $\mathbf{M},\mathbf{K}\in {\mathrm{ℝ}}^{N_s\times N_s}$ are the mass and stiffness matrices resulting from the finite-element discretization, $\mathbf{u}\in {\mathrm{ℝ}}^{N_s}$ is the nodal displacement vector, and ${\mathbf{F}}_{\text{aero}}\in {\mathrm{ℝ}}^{N_s}$ is the nodal aerodynamic load vector. With no loss in generality, we omit the structural damping in Eq. (1) along with other external loads acting on the structure. For the sake of conciseness, Eq. (1) omits the dependency of the structural matrices and aerodynamic loads on the design variables $\mathbf{d}\in {\mathrm{ℝ}}^{N_d}$. Therefore, flutter and post-flutter analysis methods are presented for a fixed design and all the structural and aerodynamic quantities are updated at each optimization step based on the current values of the design variables.

Due to the sheer number of structural DOFs, it is common to reduce the computational effort by rewriting Eq. (1) in terms of a reduced set of $N_r$ generalized coordinates (modal amplitudes) where $N_r\ll N_s$. The displacement field can be then approximated by

$$\mathbf{u}(x,y,z;t)\approx {\mathbf{\Phi }}_r(x,y,z)\mathbf{q}\left(t\right)\text{,}$$

where $\mathbf{q}\in {\mathrm{ℝ}}^{N_r}$ is the vector of retained generalized coordinates and ${\mathbf{\Phi }}_r\in {\mathrm{ℝ}}^{N_s\times N_r}$ is a matrix whose columns contain the corresponding eigenvectors (mode shapes). Substituting Eq. (2) into Eq. (1) and pre-multiplying by ${\mathbf{\Phi }}_r^T$ yields

$${\mathbf{M}}_r\ddot{\mathbf{q}}\left(t\right)+{\mathbf{K}}_r\mathbf{q}\left(t\right)-{\mathbf{F}}_{\text{aero},r}\left(t\right)=0\text{,}$$

where ${\mathbf{F}}_{\text{aero},r}={\mathbf{\Phi }}_r^T{\mathbf{F}}_{\text{aero}}\in {\mathrm{ℝ}}^{N_r}$ is the vector of generalized aerodynamic forces and, assuming that the eigenvectors are $\mathbf{M}$-orthonormal, the generalized mass and stiffness matrices take the form

$${\mathbf{M}}_r={\mathbf{\Phi }}_r^T\mathbf{M}{\mathbf{\Phi }}_r={\mathbf{I}}_r\in {\mathrm{ℝ}}^{N_r\times N_r},$$$${\mathbf{K}}_r={\mathbf{\Phi }}_r^T\mathbf{K}{\mathbf{\Phi }}_r=\text{diag}\left\{{\omega }_i^2\right\}\in {\mathrm{ℝ}}^{N_r\times N_r},$$

where ${\omega }_i$ is the i-th natural angular frequency ($i=1,\dots ,N_r$).

Assuming linear aerodynamics, the generalized aerodynamic forces are written in the Laplace domain as

$${\hat{\mathbf{F}}}_{\text{aero},r}\left(s\right)=q_{\infty }{\mathbf{Q}}_r\left(s\right)\hat{\mathbf{q}}\left(s\right)$$

where s is the Laplace variable, $q_{\infty }=\rho U^2/2$ is the flight dynamic pressure corresponding to the density ρ and flight speed U, ${\mathbf{Q}}_r\left(s\right)$ is the generalized aerodynamic force (GAF) matrix, and $\hat{\mathbf{q}}\left(s\right)$ is the Laplace transform of $\mathbf{q}\left(t\right)$. Taking the Laplace transform of Eq. (3) yields

$$\left[s^2{\mathbf{M}}_r+{\mathbf{K}}_r-q_{\infty }{\mathbf{Q}}_r\left(s\right)\right]\hat{\mathbf{q}}\left(s\right)=0\text{,}$$

or

$$\left[{\left(\frac{U}{b}\right)}^2{p}^2{\mathbf{M}}_r+{\mathbf{K}}_r-q_{\infty }{\mathbf{Q}}_r\left(p\right)\right]\hat{\mathbf{q}}\left(p\right)=0\text{,}$$

where $p=sb/U=g+ik$ is the nondimensional Laplace variable, b is the reference half chord, g is the nondimensional damping, and k is the reduced frequency [83].

Eq. (7) requires computing the GAF in the Laplace domain. However, the GAF is typically given as a transcendental function, ${\mathbf{Q}}_r\left(ik\right)$, of the reduced frequency. Two approaches are used for overcoming this problem [10]. One approach is to approximate ${\mathbf{Q}}_r\left(p\right)\approx {\mathbf{Q}}_r\left(ik\right)$ and solve the flutter equation, Eq. (6), or an equivalent form by computing the GAF matrix in the reduced frequency domain while enforcing the condition $\mathrm{\Im }\left(p\right)=k$. An alternate approach is to obtain a rational function approximation (RFA) of ${\mathbf{Q}}_r\left(ik\right)$ and use analytic continuation [105] to extend its domain from the imaginary axis (reduced frequency) to the entire complex plane (nondimensional Laplace variable). The aeroelastic system can be recast in state-space form by introducing additional aerodynamic states, such that flutter can be analyzed using a standard eigenvalue analysis.

When the GAF is represented as a transcendental function of k, flutter analysis is performed using iterative or non-iterative methods that either compute the true damping only at the flutter point (k-method) or at all flight speeds or dynamic pressure values (root locus, p-, $pk$-, and g-methods) [[83], [84], [85],87,106]. As an example, using the $pk$-method, assuming ${\mathbf{Q}}_r\left(p\right)\approx {\mathbf{Q}}_r\left(ik\right)$, Eq. (7) is rewritten as [83],

$$\left[{\left(\frac{U}{b}\right)}^2{p}^2{\mathbf{M}}_r+{\mathbf{K}}_r-q_{\infty }{\mathbf{Q}}_r\left(ik\right)\right]\hat{\mathbf{q}}\left(p\right)=0\text{.}$$

This is a second-order nonlinear eigenvalue problem, where the nonlinearity stems from the dependency of ${\mathbf{Q}}_r$ on the imaginary part of p.

Several decompositions of the GAF ${\mathbf{Q}}_r$ in Eq. (8) exist in the literature. Stanford [107] summarized these decompositions and showed that they predict the same flutter speed, but that the mode migration and characteristics may be very different. This directly impacts the optimization, resulting in different optimal designs.

As an example, a possible approach is to split the GAF into its real and imaginary parts as ${\mathbf{Q}}_r={\mathbf{Q}}_r^R+i{\mathbf{Q}}_r^I$. Assuming small damping, which gives $p/\left(ik\right)\approx 1$ [108], Eq. (8) can be rewritten in first-order form as the generalized nonlinear eigenvalue problem [84,109,110],

$$p\left[\begin{array}{cc}{\mathbf{I}}_r& 0\\ 0& {\left(\frac{U}{b}\right)}^2{\mathbf{M}}_r\end{array}\right]\left\{\begin{array}{c}\hat{\mathbf{q}}\\ p\hat{\mathbf{q}}\end{array}\right\}-\left[\begin{array}{cc}0& {\mathbf{I}}_r\\ -\left({\mathbf{K}}_r-q_{\infty }{\mathbf{Q}}_r^R\right)& \frac{q_{\infty }}{k}{\mathbf{Q}}_r^I\end{array}\right]\left\{\begin{array}{c}\hat{\mathbf{q}}\\ p\hat{\mathbf{q}}\end{array}\right\}=0\text{,}$$

where ${\mathbf{I}}_r\in {\mathrm{ℝ}}^{N_r\times N_r}$ is an identity matrix. For any given value of k, Eq. (8) (or Eq. (9) with the above specific form of ${\mathbf{Q}}_r$) can be solved using standard eigenvalue routines such as LAPACK [111] to compute the eigenvalues and eigenvectors. However, valid roots need to satisfy the equivalence $\mathrm{\Im }\left(p\right)=k$. One approach to this problem is to use an iterative solution algorithm as suggested by Hassig [83], similar to what is implemented in MSC Nastran [108]. As a disadvantage, iterative methods may converge to incorrect roots or not converge at all [112,113].

In gradient-based optimization, robustness of analysis methods, continuity, and smoothness are of utmost importance. For this reason, noniterative methods [87,113] that mitigate the aforementioned issues were more recently proposed and applied in optimization settings [109,110,114].

On the other hand, the GAF can be approximated by numerically fitting data in the reduced frequency domain using a RFA. This includes polynomial terms up to the quadratic term for representing aerodynamic stiffness, damping, and mass effects along with a rational term representing unsteady aerodynamic lag effects. Once a RFA is computed, its validity can be extended from the reduced frequency domain (imaginary axis of the complex plane) to the nondimensional Laplace domain (entire complex plane) using analytic continuation [105].

Several RFA approaches exist in the literature. Roger [115] proposed the following widely used form (given in nondimensional Laplace domain below),

$$\mathbf{Q}\left(p\right)={\mathbf{A}}_0+p{\mathbf{A}}_1+p^2{\mathbf{A}}_2+\sum _{l=3}^{N_l+3}\frac{p{\mathbf{A}}_l}{p+{\beta }_{l-2}}\text{,}$$

where $N_l$ is the number of lag terms (aerodynamic roots) ${\beta }_{l-2}$, which are specified by the user based on the reduced frequency range of interest. A common choice is to assume $N_l=3$ [116]. The coefficient matrices ${\mathbf{A}}_0,\dots ,{\mathbf{A}}_n\in {\mathrm{ℝ}}^{N_r\times N_r}$ are determined using a least-square fitting of the GAF as a function of k [115], whose accuracy depends on the number of aerodynamic roots, as well as their values.

Karpel [116] reviewed early RFA methods and proposed a new technique based on the minimum number of aerodynamic states ($N_m$), corresponding to the approximation

$${\mathbf{Q}\left(p\right)={\mathbf{A}}_0+p{\mathbf{A}}_1+p^2{\mathbf{A}}_2+\mathbf{D}[p\mathbf{I}-\mathbf{R}]}^{-1}\mathbf{E}p\text{,}$$

which includes Eq. (10) as a special case [117]. The interpolating matrices ${\mathbf{A}}_0,{\mathbf{A}}_1,{\mathbf{A}}_2\in {\mathrm{ℝ}}^{N_r\times N_r}$, $\mathbf{D}\in {\mathrm{ℝ}}^{N_r\times N_m}$, $\mathbf{R}\in {\mathrm{ℝ}}^{N_m\times N_m}$, and $\mathbf{E}\in {\mathrm{ℝ}}^{N_m\times N_r}$ are again computed by fitting the GAF data at different reduced frequencies. However, the least-square fitting requires to solve a system of nonlinear algebraic equations due to unknown coefficients in the denominator, which results in an iterative process. Morino et al. [118] later proposed a simpler and easier approach to RFA based on an interpolative structure similar to the one proposed by Karpel [116] but requiring the solution of a linear system. Furthermore, they also proposed higher-order RFA models. More recently, Ripepi and Mantegazza proposed an improved RFA interpolation structure for state-space aeroelastic stability and gust response analysis based on three nonlinear least-squares identification techniques combined with a model-order reduction [119].

Eversman and Tewari [120] discussed methods to optimize the aerodynamic lag using gradient-free or gradient-based approaches for improving the fitting accuracy. They also argued that the p in the numerator of the rational term in Eq. (10) can be omitted, resulting in the following alternate form

$$\mathbf{Q}\left(p\right)={\mathbf{A}}_0+p{\mathbf{A}}_1+p^2{\mathbf{A}}_2+\sum _{l=3}^{N_l+3}\frac{{\mathbf{A}}_l}{p+{\beta }_{l-2}}.$$

They optimized the aerodynamic lag parameters and minimized the number of these parameters for a given level of accuracy using a gradient-free optimization method. Pasinetti and Mantegazza [121] also used optimization to obtain a state-space aeroelastic model reduced to a minimum number of states. More recently, Nissim [122,123] formulated a minimum-state approximation as a nonlinear optimization problem considering the aerodynamic lag terms and the two matrices that operate on them as design variables.

Using a RFA of the GAF and analytic continuation, the aeroelastic system Eq. (7) can be rewritten in the state-space form $\dot{\mathbf{x}}=\mathbf{A}\mathbf{x}$. The state vector ${\mathbf{x}=[{\mathbf{q}}_r^T,{\dot{\mathbf{q}}}_r^T,{\mathbf{a}}^T]}^T\in {\mathrm{ℝ}}^N$ includes the generalized coordinates ${\mathbf{q}}_r$, their derivatives ${\dot{\mathbf{q}}}_r$, and $N_a$ additional aerodynamic states $\mathbf{a}$ (whose number varies depending on the type of RFA used), for a total number of $N=2N_r+N_a$ states. Representing the aeroelastic system in terms of a linear time-invariant model with RFAs allows the aeroelastic eigenvalues to be computed using standard linear eigenvalue analysis techniques (root locus) and it facilitates the application of control theory for active flutter suppression [10]. Moreover, it enables state-space time-domain aeroservoelastic response analysis and facilitates synthesizing load alleviation control laws. RFAs have also been used in aeroservoelastic stability and response sensitivity computation [[124], [125], [126], [127], [128]].

The above treatments are typically applied to linear systems. For these systems, stability is analyzed about the undeformed configuration by describing the structure using a detailed linear finite element model (FEM) and by computing the GAF using linear potential-flow lifting-surface methods like the doublet-lattice method (DLM) [129] or panel methods [8,18]. In the presence of nonlinear aerodynamic effects, GAF matrices are computed by considering small perturbations with respect to the nonlinear background flow solution for each flight condition.

One common approach to this problem is to use DLM [129] augmented with higher-fidelity numerical data or experimental measurements [49]. An alternate, more accurate approach is to use time-linearized CFD [57] to identify the GAF matrix while capturing static nonlinear aerodynamic effects directly. If structural nonlinearities are moderate, as in the case of modern transonic transport vehicles, this statically nonlinear aerodynamic description can be still combined with linear FEM structural models. However, when nonlinear effects in the structure become significant, the equations of motion need to be linearized about the aeroelastic equilibrium configuration at each flight condition to retain the effects of large static deflections or background nonlinear steady flow. The eigenvalues of the linearized aeroelastic system are then computed for each flight condition to determine the flutter boundary, that is, the flight conditions where damping vanishes [12].

Due to the high computational cost of nonlinear aerostructural solutions and the number of flight conditions to be analyzed in the flutter search, structural nonlinearities are frequently captured with beam models in nonlinear flutter analyses [7,12,61,68,69,130]. These models incur low computational cost, while still being able to capture the behavior of very flexible HARW vehicles. Low-fidelity aerodynamic formulations typically used in geometrically nonlinear flutter analyses are potential-flow theories with tip and compressibility corrections [7,12,68,69,130], vortex-lattice methods [131], or surrogate models [61].

The previous section was a brief introduction to eigenvalue-analysis methods. We now discuss the critical aspects to be considered when using these methods for constraining flutter in design optimization. As previously mentioned, we mainly refer to gradient-based optimization methods, which require the evaluation of the constraint values and their derivatives with respect to the design variables.

The process of solving for the eigenvalues and computing derivatives by differentiating the eigenvalue problem analytically is widely used [124,[132], [133], [134]] and it is included in several commercial analysis and design software packages [[135], [136], [137]]. Depending on the form of the flutter equation, different levels of accuracy are considered when taking the derivative of the eigenvalue problem. If mode shapes are used to reduce the model order, the fixed-mode approximation [134] is often used, which neglects the derivatives of the mode shapes with respect to the design variables. Such practice is commonly applied for the sake of efficiency as these derivatives, depending on the number of structural DOFs and mode shapes, can be very expensive to compute. Examples of this approximation are found in recent literature for structural optimization (with no planform changes) [109,138,139] and topology optimization [114,140]. However, ignoring the derivatives of the mode shapes can result in inaccurate gradients [134].

Another aspect that is critical to the efficacy of gradient-based optimization is the smoothness of the constraint function, which should be at least $C^1$ continuous. Several considerations are necessary to formulate a smooth flutter constraint function. A naive approach would be to specify the flutter point directly. However, this may introduce discontinuities in the constraint value between two consecutive design iterations, ${\mathbf{d}}^i$ and ${\mathbf{d}}^{i+1}$, due to mode switching or to a hump mode becoming active at a significantly lower speed [3]. Therefore, additional steps are required to ensure a continuous constraint.

Mode switching occurs when the mode that first becomes unstable (that is, at the lower speed) changes between two consecutive design iterations ${\mathbf{d}}^i$ and ${\mathbf{d}}^{i+1}$. An example of mode switching is shown in Fig. 7(a), where the hypothetical damping of a system with two modes is plotted with respect to speed. Mode switching causes a $C^1$ discontinuity of the flutter point, which poses a challenge to gradient-based optimizers [141]. The imaginary part of the eigenvalue (frequency) also switches, and in many cases the frequencies coalesce, causing a mode to become unstable. A more serious problem is when a hump mode is present in design ${\mathbf{d}}^i$ and it becomes the critical mode in the new design, ${\mathbf{d}}^{i+1}$, as shown in Fig. 7(b). The constraint value experiences a $C^0$ discontinuity, which is even more challenging for gradient-based optimizers.

Techniques exist to mitigate these problems and they are summarized by Stanford et al. [138]. Frequency-separation constraints proposed by Langthjem and Sugiyama [142] and also by Odaka and Furuya [143] can prevent mode switching by enforcing the critical mode to remain the same. This approach is illustrated in Fig. 8 (left) for a hypothetical case. The disadvantage of this method is that $N_m-2$ constraints are needed for a case with $N_m$ modes (with the expectation that two modes coalesce, and hence no constraint is needed for two modes). Furthermore, specifying the unstable mode could over-constrain the optimization process. A more serious flaw with this approach is that a hump mode is still possible.

Other techniques exist to handle both mode switching and hump modes. One approach is to enforce the real part of each eigenvalue to remain below a preset bounding curve. Such a bounding curve is proposed by Stanford et al. [144], which is a modified version of the one proposed by Ringertz [141]. The boundary ($G\left(U\right)$ in Fig. 8) spans the operating conditions of interest from wind-off to some maximum speed. This approach mitigates the issues mentioned above and has further benefits: 1) No constraint is placed on the flutter point itself, so there is no need to compute it explicitly; 2) Modes that become unstable abruptly (hard flutter) and have steeper slope than the boundary are handled as well.

Although the above approach mitigates the discontinuity issues, it requires one to apply a constraint for each speed increment and mode. Thus, $N_U$ speed increments and $N_r$ modes results in $N_r{N}_U$ constraints, dramatically increasing the total number of constraints. To mitigate the high cost of computing derivatives for so many constraints, the active set method can be used. This considers the full set of points, but reduces them to a smaller set before evaluating derivatives based on the constraint value. This approach was applied by Ringertz [141] and was also applied in the time domain by Kang et al. [145]. An alternate approach to reduce the number of constraints was proposed by Haftka [146], who suggested replacing parametric constraints by minimum-value constraints. This reduces the number of constraints to the total number of modes for the entire flight envelope. The approach can further benefit from the aggregation of the set of constraints for the individual modes into a single constraint using for example the Kreisselmeier–Steinhauser (KS) function [[147], [148], [149]].

Other classical constraint aggregation techniques can be used, such as the p-norm function or induced aggregates [150]. Kennedy and Hicken [150] compared the accuracy of classical aggregation functions with their proposed induced aggregates. They demonstrate that the induced aggregates provide better accuracy and properties compared to the classical aggregation functions (KS and p-norm). Lambe et al. [24] studied the accuracy and efficiency of the above aggregation methods for structural mass minimization of a wingbox subject to failure constraints. The induced exponential aggregation was shown to be more accurate compared to the classical constraint aggregation functions [24].

Eigenvalue-analysis methods were originally introduced for linear flutter analysis and later applied to systems with structural or aerodynamic nonlinearities. Due to the growing interest in analyzing nonlinear aeroelastic systems, methods based on Hopf bifurcation theory that compute the flutter point directly in the frequency domain while enforcing the nonlinear aeroelastic equilibrium equations were more recently used. Approaches were also developed to compute derivatives for direct methods to use them for the formulation of constraints in gradient-based optimization (see Section 4.2).

Here, we introduce the formulation of the Hopf-bifurcation approach and highlight recent work that contributed to its development. In particular, we focus on efforts that compute the Hopf point directly [151]. This is in contrast to indirect methods that monitor eigenvalues, look for sign changes, and then perform a local search (using a Newton–Raphson, secant, bisection, or interpolation method) [19], similar to what was discussed in Section 4.1.1.

Consider a generic nonlinear aeroelastic system with the dynamics governing equations written as

$$\dot{\mathbf{x}}=\mathbf{f}\left(\mathbf{x};U\right)\text{.}$$

The state vector ${\mathbf{x}=[{\mathbf{q}}^T,{\dot{\mathbf{q}}}^T]}^T\in {\mathrm{ℝ}}^N$ includes the generalized coordinates assumed to describe the system and their derivatives, for a total number of N states, $\mathbf{f}\in {\mathrm{ℝ}}^N$ is a nonlinear vector field, and U is the flight speed (the bifurcation parameter). The system Jacobian matrix is $\mathbf{A}:=\partial \mathbf{f}/\partial \mathbf{x}\in {\mathrm{ℝ}}^{N\times N}$. A Hopf bifurcation point (flutter point) occurs when the flight speed reaches the critical value $U_F$, which corresponds to an equilibrium state ${\mathbf{x}}_F$ such that $\mathbf{f}\left({\mathbf{x}}_F;U_F\right)=0$. The Jacobian matrix evaluated at the Hopf point $\mathbf{A}\left({\mathbf{x}}_F;U_F\right)$ has a pair of purely imaginary eigenvalues $\pm i{\omega }_F$, where ${\omega }_F$ is the flutter angular frequency. At this point, they cross the imaginary axis, causing the Hopf bifurcation. As explained in Section 3.2, this leads to an LCO.

Let the right complex eigenvectors associated with the bifurcating eigenvalues be $\mathbf{u}={\mathbf{u}}_r+i{\mathbf{u}}_i$, and $\overline{\mathbf{u}}={\mathbf{u}}_r-i{\mathbf{u}}_i$, where ${\mathbf{u}}_r,{\mathbf{u}}_i\in {\mathrm{ℝ}}^N$. The critical eigenvectors are computed by solving

$$\left(\mathbf{A}-i\mathbf{I}{\omega }_F\right)\mathbf{u}=0,$$

where $\mathbf{I}\in {\mathrm{ℝ}}^{N\times N}$ is an identity matrix, and using the normalization condition ${\mathbf{c}}^T\mathbf{u}=i$, where $\mathbf{c}\in {\mathrm{ℝ}}^N$ is arbitrary and constant. The flutter point is given by the solution of the nonlinear algebraic system,

$$\{\begin{array}{l}\mathbf{f}\left({\mathbf{x}}_F;U_F\right)=0\hfill \\ \mathbf{A}{\mathbf{u}}_r+{\omega }_F{\mathbf{u}}_i=0\hfill \\ \mathbf{A}{\mathbf{u}}_i-{\omega }_F{\mathbf{u}}_r=0\hfill \\ {\mathbf{c}}^T{\mathbf{u}}_r=0\hfill \end{array},$$

where the $3N+2$ unknowns are the components of ${\mathbf{x}}_F$, ${\mathbf{u}}_r$, and ${\mathbf{u}}_i$, the flutter speed $U_F$, and the angular frequency ${\omega }_F$. This system can be re-written in the compact form as

$$\mathbf{F}\left(\mathbf{X}\right)=0,\text{where}\mathbf{X}={\left[\begin{array}{ccccc}{\mathbf{x}}_F^T& {\mathbf{u}}_r^T& {\mathbf{u}}_i^T& U_F& {\omega }_F\end{array}\right]}^T,$$

and solved iteratively via the Newton–Raphson method. Starting from a first guess ${\mathbf{X}}^{\left(0\right)}$ for the solution vector, at the nth iteration this is updated by solving

$${\mathbf{J}}^{\left(n\right)}\text{Δ}{\mathbf{X}}^{\left(n\right)}=-{\mathbf{F}}^{\left(n\right)}$$

where $\mathbf{J}=\partial \mathbf{F}/\partial \mathbf{X}$ is the Jacobian of the extended system and $\text{Δ}{\mathbf{X}}^{\left(n\right)}$ is the Newton–Raphson update. Once the linear system, Section 17, is solved, the solution vector at the next iteration is updated as ${\mathbf{X}}^{\left(n+1\right)}={\mathbf{X}}^{\left(n\right)}+\text{Δ}{\mathbf{X}}^{\left(n\right)}$ and the process is repeated until convergence is reached.

While the Newton–Raphson method usually exhibits fast convergence, it is not without some disadvantages. In particular, this method requires a suitable initial guess ${\mathbf{X}}^{\left(0\right)}$ that is sufficiently close to the Hopf point in order to converge successfully [152]. Furthermore, this method may converge to other Hopf points, such as divergence [153], or to a zero velocity.

A common strategy to overcome these issues is to solve the eigenvalue problem using large $\text{Δ}U$ steps and track when one eigenvalue crosses the real axis. Once the crossing has been identified, the direct Newton method is applied to identify the Hopf point, which typically converges in few iterations [154,155].

Computing the Jacobian $\mathbf{J}$ accurately is also a challenge. Various techniques can be employed to approximate the Jacobian when solving Eq. (17): approximating it entirely using finite differences, approximating or omitting individual entries at the expense of convergence, or reducing the size of the Jacobian by taking advantage of symmetries in the problem [152,153,155,156]. For practical aircraft configurations modeled by CFD and FEM, the resulting Jacobian can be a large sparse system (which requires a preconditioning strategy) that is also ill-conditioned [153,157].

Hui and Tobak [158] were among the first authors to apply the Hopf-bifurcation method to investigate the flutter behavior of nonlinear aeroelastic systems. They analyzed the stability of pitch motion about a large mean angle of attack and showed that the aerodynamic damping vanishes when the mean angle of attack exceeds a critical value.

Morton and Beran [152] studied a two-degree-of-freedom airfoil undergoing plunge and pitch in inviscid transonic flow. The aerodynamics was modeled with an Euler CFD solver, while the structural model consisted of linear and torsional springs for the plunge and pitch DOFs, respectively. They computed the flutter point directly using the Hopf-bifurcation method and validated it with time-accurate analyses. The Jacobian was computed using second-order finite differences and the system, Eq. (16), was solved using a Newton–Raphson method. They showed that the approach is computationally efficient compared to full time-accurate simulations, and that it produces accurate results for inviscid transonic flow.

Badcock et al. [159] extended their previous work [153] and applied a modified Newton–Raphson method to solve the Hopf-bifurcation problem. The method was demonstrated on the AGARD 445.6 aeroelastic wing [160] by computing the full flutter boundary in an inviscid transonic flow. For this configuration, the direct approach was up to two orders of magnitude more efficient than time marching.

Woodgate and Badcock [161] extended this work and analyzed three geometries: the Goland wing [162], a transonic wing-body configuration [163], and the AGARD 445.6 wing. The direct method reduced the computational cost by one to three orders of magnitude compared to time-accurate simulations.

More recently, Badcock and Woodgate [157] further improved the approach for large-scale CFD-based systems coupled with FEM structures. They used a Schur-complement eigenvalue formulation to treat parts of the system Jacobian in a parallel computing environment. This procedure eliminated the process of solving large sparse systems, which are almost singular, further increasing the method efficiency. This method shows similarities to the eigenvalue methods applying mode tracking discussed in Section 4.1.1. The approach was applied to the Goland wing, a supersonic transport arrow rudder buzz in transonic flow [164], the wing-body configuration mentioned above [163], and a generic fighter configuration based on publicly available F-16 data [157,165]. The computational cost of the proposed method varied from 5 to 60 times that of a steady-state solution.

In the context of gradient-based optimization, formulating a successful flutter constraint using the direct method needs careful consideration. Possible discontinuities due to competing flutter mechanisms in the design process, such as mode swapping or hump modes, described in Section 4.1.2, may confuse the gradient-based optimizer. However, techniques addressing these concerns in the context of direct methods are deficient or nonexistent.

Despite early work by researchers such as Haftka [3,146] and Hajela [166] optimization subject to flutter constraints is still not a standard design practice. More recently, several authors have integrated flutter constraints into design and investigated their effect on the optimal solutions. These efforts are summarized in Table 1 and reviewed below. They differ in the use of eigenvalue analysis (EV), direct (Hopf), or time-domain (TD) prediction methods for flutter analysis, the fidelity of structural and aerodynamic models, and the optimization problem formulation in terms of objective, type and number of design variables, and use of gradient-based (GB) or gradient-free (GF) solution algorithms. Efforts that used gradient-based algorithms also differ in the methods used for computing derivatives with respect to design variables. Finally, previous work shown in Table 1 differs in the types and number of design variables that were considered in the optimization problem. Most of the efforts included only structural sizing variables without considering changes in the aerodynamic shape. Fewer efforts optimized the airfoil or planform shapes. Including planform shape variables is challenging because changes in the mode shapes and the corresponding natural frequencies need to be considered when computing derivatives, incurring additional computational cost.

In one of the early efforts to constrain flutter, Turner [167] formulated a mass minimization problem by considering distribution of material rather than the structure topology to meet a specified flutter speed. Later, Bhatia and Rudisill [168] developed a numerical procedure to minimize wing mass while satisfying a flutter constraint. They applied the procedure to a uniform cross-section box beam consisting of three bays in order to minimize the mass while maintaining the flutter speed. Rudisill and Bhatia [169] improved the rate of convergence of this procedure by computing second-order derivatives of the eigenvalues and of the flutter speed with respect to the design variables. Gwin and Taylor [170] developed the method of feasible directions for the mass minimization of a structure subject to a minimum flutter speed. They were able to handle up to 33 structural design variables in a supersonic problem.

Due to the limited computational capabilities of the time, these early efforts used simple aerodynamic and structural models and employed a similar strategy to enforce the flutter constraint. They all formulated the flutter problem in the frequency domain as an eigenvalue problem, which was then differentiated with respect to the design variables. The derivatives of the eigenvalues were obtained using the left and right eigenvectors and one of these efforts also considered second-order derivatives to better guide the optimization process and improve the rate of convergence [169].

Ringertz [141] applied the k-method to minimize the weight of a cantilevered wing in incompressible flow subject to flutter and divergence constraints. The structure was modeled as a composite FEM model, while the unsteady aerodynamic loads were computed using the DLM. The eigenvalue problem was analytically differentiated in the modal space. A continuous flutter constraint was formulated using a boundary to constrain the damping values, which resulted in a large number of constraints. The method was demonstrated on a rectangular wing and on a swept wing with taper, where the objective was to minimize weight with respect to element group thicknesses subject to flutter and divergence constraints. In both cases, a considerably lighter design was achieved.

Mallik et al. [171] investigated the impact of a previously developed flutter constraint [174] on the MDO of a TBW aircraft. The flutter speed was computed using the k-method applied to a linear pre-stressed structural model. The unsteady aerodynamic loads were computed using Theodorsen theory [175] with a Prandtl–Glauert compressibility correction. They implemented an iterative procedure to ensure consistency of the flutter speed, Mach number, and altitude, and optimized a representative TBW configuration for minimum takeoff weight and fuel burn using a genetic algorithm. The flutter constraint was formulated by enforcing a minimum flutter Mach number and was added to several other mission constraints. Comparing the optimization results with those obtained by removing the flutter constraint, they showed that this is necessary to obtain a flutter-free optimal solution.

Jonsson et al. [110] developed a flutter constraint suitable for high-fidelity gradient-based optimization including wing planform variables for the first time. They developed an efficient non-iterative root finding method to compute the flutter eigenvalues, based on the $pk$ flutter equation that includes a scheme to track the mode migration between subsequent reduced frequencies and designs.

The flutter-free flight envelope was implicitly defined using a constraint curve for the damping values to prevent discontinuities in the flutter constraints and to allow the computation of derivatives. The differences between damping values and the constraint curve were aggregated into a single value using the KS function [147]. Accurate and efficient derivatives of the aggregated flutter constraint value were computed with respect to both structural sizing and wing planform variables. The derivatives were computed using a combination of analytic and automatic differentiation methods in reverse (adjoint) mode and validated against the complex-step method.

The flutter constraint was used to maximizing the range of an idealized rectangular wing (flat plate) by varying structural thickness and planform design variables. The optimal design had an increased aspect ratio and an improved range while remaining flutter and divergence free. This approach provides a way to include airfoil shape variables in the future to perform full aerostructural design optimization [176] subject to flutter constraints.

Several authors have focused on including aerodynamic nonlinearities in flutter analysis to optimize transonic configurations. Stanford et al. [144] evaluated six different novel tailoring schemes employed in mass minimization optimization. They analyzed the flutter characteristics using the $pk$ method and the commercial software ZTRAN [50] to retain aerodynamic nonlinearities. The nonlinear higher-fidelity Euler code, ZEUS, was used to compute steady background flow at multiple transonic Mach numbers for a fixed cruise shape. Using these steady-state CFD solutions as an input, the linearized unsteady loads were computed for a range of reduced frequencies using time-linearized transonic small disturbance (TSD) analyses about the equilibrium solutions.

The system damping values were forced to be under a stability boundary, similarly to the approach by Ringertz [141]. The transonic aerodynamic loads were computed offline before the optimization and their derivatives were obtained by differentiating the eigenvalue problem [134]. Flutter, stress, and buckling constraints computed in this study were all aggregated using a KS function [147]. They considered the fixed-mode derivatives to improve computational efficiency, an approach that does not allow for shape changes.

They studied the undeflected Common Research Model (uCRM) wing [176] and obtained six different optimal wing structures corresponding to different tailoring schemes, all for the same operating condition and setup. The six tailoring schemes considered for the structural design were metallic thickness variations, functionally graded materials, balanced or unbalanced composite laminates, curvilinear tow steering, and distributed trailing edge control surfaces. While there was a structural wing mass reduction for every case, many of the lighter designs had an active flutter constraint, while the buckling constraint was active for the heavier cases.

Chen et al. [137] extended their previous work by [136] computing derivatives of flutter constraints with respect to shape variables using ZEUS coupled with a boundary layer code. The derivatives with respect to shape were computed using the complex-step approach [28,31], which is numerically exact. The flutter constraint was formulated using the g-flutter method [87], which was differentiated analytically with respect to the design variables. The approach was verified for a cantilevered planform, similar to the F-5 geometry [177,178], where a structure consisting of 10 spars, 10 ribs, and upper and lower skins was modeled using MSC Nastran. While no optimization results were presented, the flutter derivatives were verified against exact results.

Bartels and Stanford [139] proposed an approach to enforce a CFD-based flutter constraint for gradient-based structural optimization of transonic vehicles. The flutter analysis was performed as a standard eigenvalue analysis on the state-space representation of the aeroelastic system obtained via RFA [115]. The GAF matrix of the baseline structure was identified from time-linearized unsteady RANS simulations about nonlinear steady-state solutions [179]. The GAF matrix of the updated design was then computed by projecting the updated modes onto the baseline modes. The methodology was used to minimize the uCRM mass subject to structural and aeroelastic constraints. The flutter was constrained by requiring that the real part of the system eigenvalues be below a bounding curve [166]. The optimization assumed a fixed-mode approximation. The authors compared the optimal solutions obtained using unsteady RANS and DLM aerodynamics in the optimization loop. They showed that the DLM-based solution was not conservative and had a significantly different thickness distribution compared to the CFD-based optimal design.

Opgenoord et al. [65] extended a low-order two-dimensional transonic flutter prediction model [64] to wings and implemented the model into a conceptual aircraft design tool to investigate the impact of geometric parameters and Mach number on the flutter boundary. Furthermore, they optimized the D8.0 configuration by minimizing the maximum take-off weight and fuel burn with and without a transonic flutter constraint. Enforcing the flutter constraints resulted in lower optimal aspect ratio and a weight penalty or lower fuel burn reduction compared to the case when the flutter constraint was omitted. Opgenoord et al. [66] also used the developed flutter model to optimize the internal lattice structure of a wing by minimizing weight with and without enforcing a flutter constraint in addition to stress and buckling constraints. The optimal design including the flutter constraint showed only a slight mass increase thanks to the appropriate aeroelastic tailoring of the lattice structure.

While several authors have performed flutter-constrained optimizations using nonlinear aerodynamic models, examples of flutter constraints considering nonlinear structures are more rare due to the more recent interest in optimizing very flexible aircraft. Variyar et al. [7] developed a framework for MDO of geometrically nonlinear aircraft subject to flutter constraints by coupling the SUAVE design tool [180] with the ASWING nonlinear aeroelastic solver [181]. They developed an interface to convert the arbitrary aircraft designs output by SUAVE into equivalent nonlinear beam models for ASWING to drive structural and aeroelastic analyses and to post-process the results. The flutter speed was obtained iteratively by computing the eigenvalues of the statically deformed aircraft at different flight conditions.

This MDO framework was used to optimize the Sugar VOLT strut-braced aircraft [182] for minimum fuel burn subject to mission, maneuver, gust, and flutter constraints. The flutter constraint was implemented by imposing a minimum flutter speed, and derivatives with respect to the design variables were obtained via finite differences. The authors performed three MDO cycles by adding the maneuver, gust, and flutter constraints to the mission constraints. Despite having better performance, the optimal solution achieved with only maneuver and gust constraints experienced flutter within the flight envelope, highlighting the need for a flutter constraint in the design process.

Lupp and Cesnik [130] studied the effect of a flutter constraint including geometrical nonlinearities on the design of a BWB. They extended the University of Michigan's Nonlinear Aeroelastic Simulation Toolbox (UM/NAST) to utilize AD in reverse mode to determine coupled aeroelastic derivatives including geometrical nonlinearities. They proposed an algorithm to increase the computational efficiency of the gradient evaluation for a geometrically nonlinear aeroelastic analysis. The sample optimization formulation included a geometrically nonlinear beam-based flutter constraint using a KS aggregation [110] to obtain a scalar constraint for the entire flight envelope. The authors ran three fuel burn minimization cases: with a strength constraint, with a linear flutter constraint, and with a geometrically nonlinear flutter constraint with the wing chord distribution, wing box size, and wing box thickness as design variables. While the linear flutter constraint became active over the strength constraint, it was not conservative compared to the geometrically nonlinear constraint. They concluded that a geometrically nonlinear flutter constraint is needed for very flexible aircraft.

Xie et al. [131] used a flutter constraint to minimize the weight of a very flexible wind tunnel model. The purpose of this constraint was to ensure flutter within the wind tunnel speed range. The optimization problem coupled a geometrically nonlinear beam solver with a vortex lattice code for the static aeroelastic analysis, while a doublet lattice code was used to compute the unsteady loads. Since they use a gradient-free algorithm, no derivatives were computed. They compared optimization results based on linear and geometric nonlinear beam models subject to a flutter constraint. In addition to flutter, tip displacement and torsion angle constraints were also enforced. The linearly optimized configuration resulted in a wing lighter than the optimal solution obtained with the nonlinear flutter constraint. Furthermore, the flutter and displacements constraints were violated when the linear optimized configuration was analyzed considering geometric nonlinear effects. The difference in the results highlighted the need for a flutter constraint when optimizing very flexible vehicles and the importance of using nonlinear flutter prediction methods not only for analysis, but in design optimization as well.

Bhatia and Beran [172] developed a framework to optimize thermally stressed nonlinear structures subject to transonic flutter constraints. They showed that including aerothermoelastic static nonlinearities in the flutter analysis impacts the optimal design. These effects are particularly important for high-speed vehicles, which are subject to significant thermoelastic stresses when flying through the transonic regime during reentry [183]. The authors optimized a skin panel with respect to the thickness and density distributions to minimize the mass subject to a flutter constraint. The structure was modeled using a Timoshenko beam finite element with nonlinear von Kármán strain, while the transonic flow was solved via a finite-element discretization of the Euler equations. The structure was linearized around the static thermoelastic response and the aerodynamics was linearized around the background steady transonic flow past the baseline geometry. Flutter was analyzed in the frequency domain as an iterative V-g solution [1] considered the system linearized around the nonlinear equilibrium configuration for each operating condition. The flutter speed was constrained directly. The optimal solution obtained by linearizing around the thermally stressed equilibrium configuration showed a mass lower than the result using an unstressed analysis.

The Hopf-bifurcation method has also been applied in design optimization. Kennedy et al. [156] proposed a variant of previous bifurcation approaches [157,159] to optimize an aeroelastic system subject to a flutter constraint, which was formulated in terms of flight speed rather than damping. The proposed method had the benefit of not requiring a search of the flutter point (which may be located well outside of the flight envelope) at each design iteration. The KS function [147] was used to smooth the effect of mode switching in the constraint value by aggregating the less damped modes, which yielded smooth gradients. The authors presented preliminary results for a medium fidelity three-dimensional aerodynamic panel code coupled with the structural finite-element code TACS [184]. Although no detailed optimization was performed, they performed a preliminary study on the uCRM benchmark [176].

Beran et al. [173] developed a fast adjoint method to compute derivatives of flutter points computed via the Hopf-bifurcation method for gradient-based MDO. The approach was applied to the highly flexible cantilevered wing studied by Tang and Dowell [81,82]. Both geometric nonlinearities due to the structure and aerodynamic effects described using the ONERA stall model [82] were considered when computing the flutter point and its derivatives with respect to aerodynamic and structural design variables. While no optimization study was presented, the authors outlined the future work required to apply the methodology: validate the derivatives, assess the computational cost compared to alternate time- and frequency-domain flutter prediction methods, and develop the handling of mode switching to avoid discontinuous flutter points.

Some rare examples of optimizations with flutter constraints based on time-accurate analyses are also found in the literature, but they have been restricted to simple problems. In an early work, Holden [185] applied a colocation method to constrain the aeroelastic response envelope for a wing optimization problem. More recently, Mani and Mavriplis [99] were among the first to present a fully coupled unsteady adjoint for aeroelastic optimization. They demonstrated their method successfully in a shape optimization of two-dimensional airfoil section to suppress flutter, using a total of 32 design variables in the form of Hicks–Henne bump functions [186]. Later, Mishra et al. [100] extended previous work and presented the fully coupled unsteady adjoint for three-dimensional aeroelastic problems, which was demonstrated in a shape optimization of a flexible rotorcraft configuration.

Zhang et al. [101] developed a coupled adjoint method for unsteady aerostructural problems solved via time simulations. The method was applied to an airfoil shape optimization problem with the goal of suppressing flutter. The aerodynamics was computed with an Euler CFD code coupled with a boundary layer code to account for viscous effects. Both the continuous and discrete coupled adjoint were developed for steady-state analyses. The discrete approach proved more promising and only this version was developed for unsteady cases. A damping objective function was proposed that used a Hilbert transform [92] of the nonlinear unsteady time-history. To achieve the required flutter boundary, the damping objective function was minimized to obtain a neutral response, indicating the flutter point. The authors applied the methodology to the optimization of the two-dimensional (2D) Isogai airfoil [37,38] to suppress flutter. Only derivatives with respect to shape variables were computed and the airfoil shape was parametrized using 48 Hicks–Henne bump functions [186]. A neutrally stable (zero-damping) configuration was obtained for a given flight condition.

Zhang et al. [102] extended their previous work [101], where the coupled adjoint was developed for time-marching simulations. Two objective functions were used: one that maximized the flutter boundary and another that matched a given flutter boundary. To maximize the flutter boundary, they minimized the squared and time averaged history of the lift coefficient. To match a given flutter boundary, the damping value of a given time history was minimized to obtain a neutral response. B-spline curves were chosen as a parametrization method due the large number of Hicks–Henne functions that were previously needed. The authors presented steady-state optimization results for a 2D airfoil optimized to match a given pressure distribution and for the 3D ONERA M6 case [187], where the objective was a composite function of lift and drag. Configurations considered in the unsteady optimization consisted of the 2D Isogai airfoil and the Goland wing as modeled by Kurdi et al. [188]. For the unsteady 2D airfoil, two optimization cases were considered: flutter margin maximization and a flutter boundary matching (i.e., a neutral response for the given flight condition). The design variables were the plunge and pitch stiffness values. For the Goland wing, they optimized the aerodynamic shape to maximize the flutter speed with respect to 120 shape variables. A structural optimization of the Goland wing was also performed to maximize the flutter speed with respect to the skin thickness. No optimization was performed using simultaneously structural sizing and aerodynamic shape variables.

Kiviaho et al. [189] developed a flutter constraint using their previously developed unsteady aeroelastic framework with adjoint sensitives [103,104]. The flutter constraint is based on a matrix pencil method [90] applied to a time-history, which estimates the damping based on most critical aeroelastic modes. Two methods are proposed, a direct flutter point evaluation which finds the flutter point where the dynamic pressure lower bound is specified as the design flight condition, and a flutter margin or clearance approach where the dynamic pressure times some tolerance is fixed. The direct method was demonstrated in a single design variable optimization of an airfoil were the dynamic pressure is minimized subject to the flutter constraint in order to identify the flutter point.

Now that we reviewed methods for flutter prediction and their applications to optimization (see Table 1), we can summarize the open problems in goal of integrating flutter constraints into aircraft design optimization.

Gradient-based optimization is the preferred choice for optimizations with respect to large numbers of design variables. When enforcing a flutter constraint in a gradient-based optimization, a serious challenge is varying structural and aerodynamic design variables simultaneously to optimize the aircraft external shape, planform, and internal sizing. Most of previous work optimized only structural sizing, and only a few efforts included airfoil shape and wing planform design variables as well. Furthermore, simplifying assumptions like the fixed-mode approximation were frequently used when computing derivatives to limit computational cost [109,138,139]. These assumptions are adequate for a structural optimization, but they can lead to inaccurate results when varying aerodynamic properties because this can cause significant changes in the mode shapes at each design iteration. Therefore, gradient-based optimizations with respect to structural, planform, and shape variables need approaches that take into account the derivatives of the mode shapes when computing the derivatives of the flutter constraints. Few examples of these approaches applied to simplified configurations or using lower-fidelity models are available in the literature [7,110,130]. However, they still have to be demonstrated on practical configurations parametrized by large number of structural and aerodynamic design variables.

A second major challenge is developing efficient flutter analysis models and methods applicable in the presence of aerodynamic or structural nonlinearities. In previous aircraft optimizations including linear flutter constraints, the natural choice was to analyze flutter in the frequency-domain due to the availability of well-established and computationally efficient eigenvalue-analysis methods. The computational cost of these methods is a big challenge when including aerodynamic or structural nonlinearities because the flutter point depends on the equilibrium state. For each flight condition considered in the flutter search, the steady background flow solution (in the case of aerodynamic nonlinearities) or the coupled aerostructural equilibrium (in the case of both aerodynamic and structural nonlinearities) needs to be determined first. Next, the linearized model about each equilibrium point must be identified for computing the aeroelastic eigenvalues and determine at which point modal damping vanished. This process must be repeated for each flight condition, while linear methods analyze flutter by considering the undeformed configuration at zero angle of attack for each flutter search point. Moreover, in the presence of nonlinear effects multiple equilibrium points may also exist for each flight condition, which further increases complexity and computational cost. Stability must be analyzed about all equilibrium points, otherwise critical constraint values may be missed.

For large high-fidelity models with both structural and aerodynamic nonlinearities, computing the aerostructural equilibrium points and the corresponding linearized systems may be computationally prohibitive. For moderately flexible configurations, structural nonlinearities can be neglected, so eliminating the need to solve a nonlinear static aeroelastic problem at each flight condition. However, transonic aerodynamic nonlinearities still require identifying a linearized unsteady aerodynamic model for each steady background flow solution. Computing and retaining aerodynamic nonlinearities to accurately predict the flutter point in the transonic regime remains a big challenge due to the large computational cost associated with CFD. To address this problem, some efforts tried to preserve the computational efficiency of lower-fidelity methods while retaining the nonlinear physics modeled by the higher-fidelity CFD methods [50,52]. Other works proposed time-linearized CFD solvers [53,57] or reduced- or low-order models calibrated using CFD [62,64]. Despite these progresses, optimizing aircraft subject to transonic flutter constraints is still an open problem, particularly when using gradient-based methods and when seeking both structural sizing and aerodynamic shape changes.

Due to these challenges, previous work mainly optimized linear aircraft configurations or included only aerodynamic nonlinearities while limiting to structural sizing. Aerostructural optimizations considering both aerodynamic and geometric structural nonlinearities used low-order models or optimized simple configurations to limit the computational cost [7,130]. No previous work considered both structural and aerodynamic nonlinearities in a high-fidelity MDO setting. Additionally, rigid-body DOFs were never included in the flutter constraints, which may lead to unfeasible designs for configurations prone to coupled rigid-elastic instabilities.

Few studies used alternatives to frequency-domain eigenvalue analysis methods for nonlinear flutter prediction, like the Hopf-bifurcation method or time-domain simulations. Such studies are rare because the high computational cost of time-accurate analyses makes them prohibitive to optimize complex configurations. For this reason, these applications were limited to simple configurations and frequently two-dimensional problems.