Multiﬁdelity Sensitivity Study of Subsonic Wing Flutter for Hybrid Approaches in Aircraft Multidisciplinary Design and Optimisation

: A comparative sensitivity study for the ﬂutter instability of aircraft wings in subsonic ﬂow is presented, using analytical models and numerical tools with different multidisciplinary approaches. The analyses build on previous elegant works and encompass parametric variations of aero-structural properties, quantifying their effect on the aeroelastic stability boundary. Differences in the multiﬁdelity results are critically assessed from both theoretical and computational perspectives, in view of possible practical applications within airplane preliminary design and optimisation. A robust hybrid strategy is then recommended, wherein the ﬂutter boundary is obtained using a higher-ﬁdelity approach while the ﬂutter sensitivity is computed adopting a lower-ﬁdelity approach.


Introduction
Within aircraft multidisciplinary design and optimisation (MDO) [1,2], efficient methods and robust tools are highly sought after as effective reduced-order models (ROMs) [3][4][5] for fast parametric aeroelastic analyses [6][7][8], especially for sensitivity and uncertainty evaluation purposes at the preliminary design stage [9][10][11]. Simplified semi-analytical formulations for the bending-torsion instabilities of flexible wings have been available for a long time but have inherent limitations [12][13][14], whereas complex fluid-structure interaction (FSI) [15] procedures based on the finite element method (FEM) [16] and computational fluid dynamics (CFD) [17] have recently been used to improve accuracy and generality, but remain computationally expensive and require special care [18,19] with significant effort to pre-and post-process the numerical simulations.
Continuing previous studies on multilevel techniques for practical aeronautical applications [20,21], this conceptual work investigates a robust hybrid strategy where the flutter boundary is accurately obtained using a higher-fidelity approach while the flutter sensitivity is efficiently computed by adopting a (tuned) lower-fidelity approach as an effective ROM. A comparative sensitivity study is hence presented for the aeroelastic stability boundary of a uniform cantilever wing in subsonic flow as the standard benchmark [22]. The Typical Section idealisation [23,24] is employed as the analytical lower-fidelity model, whereas the numerical higher-fidelity model couples a beam-based FEM with panel-based CFD. Considering Loring's wing [25] with a NACA0002 airfoil [26], the flutter analyses build on earlier publications [27,28] and encompass parametric variations of wing properties, such as its structural inertia, stiffness and aspect ratio. The effects of the latter on the aeroelastic instability are quantified and differences between lower-fidelity and higher-fidelity results are critically assessed from both theoretical and computational perspectives, in order to study the necessary trade-off between complexity and costs of model fidelity for intensive industrial MDO activities [29].

Problem Formulation
Uniform cantilever wings have long been used as the standard benchmark for fundamental studies on aeroelastic analyses and parametric sensitivities for the bending-torsion instabilities of flexible wings, both numerically and experimentally [12,22]. They feature constant material properties (i.e., structural density ρ m , Young's elastic modulus E and shear modulus G), chord c, mass m and moment of inertia µ per unit length at the inertial axis x CG (i.e., the line of sectional centres of gravity); bending stiffness EI and torsion stiffness GJ at the elastic axis x EA ≡ 0 (i.e., the line of sectional shear centers) along the semi-span l; structural damping is typically small and conservatively neglected [30].
With θ(y, t) and h(y, t) the pitch and plunge motion of the flexural axis at the time t, respectively, the wing vertical displacement w(x, y, t) is given as w = h + xθ directly, x and y being the chordwise and spanwise directions. Assuming an Euler-Bernoulli beam model [31], the linear system of coupled partial differential equations (PDEs) for the dynamic aeroelastic response of wing bending and torsion undergoing small deformations reads: with ∆L(y, t) and ∆M(y, t) being the sectional lift and pitching moment, respectively; the governing equations are then completed by both geometrical and natural boundary conditions, namely: As for the case of Loring's wing, the problem formulation assumes an isotropic material without loss of generality; an anisotropic material (e.g., for the case of composite wings [32]) would feature elastic coupling between bending and torsion (due to different mechanical characteristics in respective directions [12]) but would not introduce conceptual changes to the overall methodology and main conclusion of this work.

Uncoupled Natural Vibration Modes
Assuming h ≈ φ h χ h and θ ≈ φ θ χ θ , where χ(ω, t) are the generalised coordinates of the uncoupled vibration mode shapes φ(y) with natural frequencies ω, the latter for the cantilever wing's bending and torsion dynamics can independently be obtained by solving the relative homogeneous PDEs via separation of variables with the respective boundary conditions, namely [12]: which form a complete set of modal bases for the generalised solution of the aeroelastic equations; yet, note that these uncoupled bending and torsion modes are inherently orthogonal within their own bases but not to one another.

Unsteady Aerodynamic Model
According to thin aerofoil theory for incompressible potential flow, the sectional unsteady aerodynamic loads due to the wing motion may be written as [22]: where ρ and U are the reference air density and horizontal flow speed; the aerodynamic centre x AC (where the circulatory load acts) and the control point x CP (where the non-penetration boundary condition for the vertical flow speed V(y, t) is imposed) lay at the first-and third-quarter chords, respectively; x MC refers to the mid-chord (where the non-circulatory load acts); C l/α is the derivative of the sectional lift coefficient C l with respect to the nominal angle of attack α.
Within tuned strip theory [21], the scaling factor κ introduces steady three-dimensional downwash effects due to the tip vortices [33] and takes part in the derivative of the wing lift coefficient C L/α = κC l/α [34], while the complex lift-deficiency function C T (k) is defined in the reduced-frequency k domain and embeds the lag due to two-dimensional inflow effects from the travelling flat wake [35], namely: with being the aerofoil thickness ratio, η the wing aspect ratio, e(η) Jones' edge-velocity factor (i.e., the ratio between wing semi-perimeter and span [36]) and H 2 n (k) Hankel's functions of the second type and n-th order. Note that the proposed approximate correction for the wing lift coefficient is analogous to applying Oswald's efficiency factor [37], whereas κ ≡ 1 for standard strip theory.
In summary, unsteady aerodynamics account for both circulatory lag and non-circulatory inertia effects [38], whereas quasi-unsteady aerodynamics neglect circulatory lag terms [39] within a quasi-static approach (i.e., C T ≈ 1 with k << 1). Quasi-steady and steady aerodynamics then disregard non-circulatory inertia terms too, within a static approach (i.e., C T = 1 and k ≈ 0) where the control point coincides with the elastic axis in the former case (in order to avoid unrealistic aerodynamic damping [23]) and all time derivatives are eventually discarded in the latter case. Note that the structural motion represents a time-dependent boundary condition for the incompressible airflow in all cases, but a synchronous proportional variation of the circulatory airload occurs within quasi-static approaches only [40].

Modal Approach and Stability Analysis
Given the natural vibration modes of the wing, the relative equations of motion can be transformed into ordinary differential equations (ODEs) by employing Ritz's method [41], where generalised coordinates {χ(t)} multiply the natural vibration mode shapes φ. Recasting the whole aero-structural model in vector-matrix form, the governing modal equations then read: and enforce a monolithic FSI, with {F A (t)} the aerodynamic load; as anticipated, [C S ] = [0] is hereby assumed for convenience without loss of generality. Due to linearity, the parametric aeroelastic stability of the wing is then monitored from the root locus of the characteristic equation for the flutter determinant as the reference airspeed increases; in particular, the system becomes metastable when the real part of at least one eigenvalue λ vanishes and leaves the response undamped, namely: with λ = iω in the case of dynamic flutter where a coupled resonant harmonic motion excites the eigenmode {u} (Hopf's bifurcation occurring when a couple of complex conjugate eigenvalues cross the imaginary axis), whereas λ = 0 in the case of static divergence (occurring when the aeroelastic stiffness matrix becomes singular). Note that the eigenvalues and eigenvectors are suitably assumed to form a complete distinct set, as per most of the practical cases [42]; of course, a sufficient number of natural vibration modes shall be employed in order to grant proper convergence of the modal analysis [43].

Sensitivity Analysis
The aeroelastic stability boundary acts as a typical constraint in aircraft MDO problems, where parametric sensitivities of both flutter and divergence speeds with respect to variations in the design variables are required for gradient-based optimisation algorithms [44,45]. In general, such derivatives may numerically be obtained by performing many simulations with small alterations of each and every design variable; nevertheless, the aircraft's conceptual and preliminary design stages would require extensive robust computations for the trends of both objective function and constraints while changing the design parameters [46,47]. Effective analytical solutions (with explicit expressions in a few special cases) have hence been derived by differentiating the governing equations, taking advantage of the self-adjoint nature of complex eigenproblems and then separating real and imaginary parts [48,49].
Especially when the sensitivity of the critical mode shape is also required, the continuation method may efficiently be employed to calculate the critical point and obtain all desired derivatives at the same time [50]. The eigenproblem arising in flutter and divergence analysis may indeed be considered as a system of nonlinear equations with a normalisation condition for the eigenvectors; the solution is then numerically sought by the Newton-Raphson method, where the aeroelastic equations are differentiated with respect to a design parameter of interest p and the system is linearised at each iteration [51]. Within a single unified procedure, the eigenvalues' and (right) eigenvectors' derivatives may then automatically be determined along with the eigensolution itself at once, with no information about the transposed (left) problem being required; higher-order derivatives may be obtained by further differentiation. In particular, with U = U c (p), ω = ω c (p) and {u} = {u c (p)} in the critical flight condition, the nonlinear eigenproblem for the flutter mode reads [9]: and is differentiated to give the eigensolution sensitivity with respect to the design parameter as: where the submatrices and subvectors involve the derivatives of the aeroelastic model, namely: in the degenerate case of static divergence with Im {u c } = 0 and ω c = 0, all quantities are inherently real and the matrix of algebraic equations above specialises without the second and fourth rows and columns, respectively. Following previous studies [27], the derivatives may finally be normalised with respect to the reference values of the involved quantities, thereby comparing all sensitivities within a coherent representation.

Lower-Fidelity Model
The Typical Section abstraction [52] is hereby employed as the analytical lower-fidelity model, as it is conceptually illuminating, inherently robust and computationally efficient within its limitations. It resembles the section of the uniform wing and its fundamental aeroelastic behaviour for the bending-torsion coupling mechanism; matching the wing inertia and natural frequencies, such an abstraction, may be regarded as a condensed ROM. The structural arrangement is a classic mass-spring system for the wing section, where the coupled pitch and plunge motion is restrained by vertical and angular linear springs with equivalent stiffnesses k h = mω 2 h and k θ = µω 2 θ (see Appendix A), respectively; it is representative of the inertial and elastic properties per unit length of the wing at about 75% of its span, where the structural arrangement becomes progressively flexible and the aerodynamic load is still high. The model is hereby extended to include all bending modes having natural frequency below that of the first torsion mode, as the instability may involve any combination of them in general.
Considering the first two uncoupled bending modes and the first uncoupled torsion mode of Euler-Bernoulli's beam model (with γ 1 ≈ 1.875, γ 2 ≈ 4.694 and ν 1 ≈ 1.571) [12], the aeroelastic equations of motion for the wing's Typical Section in (time-varying) steady incompressible flow are: where the cross-projections of the non-orthogonal modes scale the off-diagonal coupling terms: in particular, all the latter are constant and read f ≈ 0.959, g ≈ 0.274 and r ≈ 39.275 for homogeneous aero-structural properties. Note that assuming a static aerodynamic approach is consistent with the rigorous application of the scaling factor κ, which is hereby based on steady lifting-line theory; no damping is hence provided to the aeroelastic response. It is also worth stressing that the equations for the first and second bending/plunge modes are not coupled directly, but the equation for the torsion/pitch mode couples them indirectly.

Aero-Structural Parametric Derivatives
The sensitivity of the aeroelastic matrices with respect to aero-structural parameters can readily be obtained in explicit analytical form, where the chain rule applies [28]. In particular, the derivatives with respect to the material density ρ m and elastic modulus E are given by: whereas those with respect to reference flow speed U and perturbation frequency ω read: and with respect to the wing semispan l it is:

Higher-Fidelity Model
The high-fidelity model is based on a FEM of a slender beam solved with the commercial code Nastran [53] for the structural analysis; the latter is coupled with either the doublet lattice method (DLM) available in the same software [54] or an in-house panel code based on Morino's boundary element method (BEM) for the steady and unsteady aerodynamic analysis [55,56], thereby enforcing a numerical FSI. Both Nastran-based structural FEM and in-house aerodynamic BEM are generated with automatised routines, in order to ease parametric studies; this capability is exploited to perform numerical convergence studies (see Appendix B) and compare the analytical flutter derivatives with their numerical counterparts obtained via finite difference.

Structural Model
Following previous works [21], the wing structure is modelled with CBEAM elements and accounts for the distance between inertial and elastic axes; the node lying at the wing root is then clamped (see Figure 1). Further FEM nodes are placed at the leading and trailing edges of the wing and connected to the beam nodes with the rigid elements RBE2, in order to support mapping with the aerodynamic grid (in terms of both structural deformation and aerodynamic load: the latter modifies the former, which changes the airflow boundary condition in turn) according to the closely-spaced rigid diaphragm (CSRD) assumption [57]. For the natural vibration analysis, shear deformation is neglected in order to obtain the Euler-Bernoulli beam model, with PBEAM defining the properties (i.e., inertia and stiffness) of the beam element and SPC1 defining the single-point constraint for the clamped root. Using Nastran's SOL103 to obtain the structural eigenvalues and eigenvectors, the vibration analysis is performed while selecting Lanczos' method [58] as available in EIGRL and normalising the modes to unit values of the generalised mass.

Aerodynamic Model
Standard flutter prediction methods in the industrial environment are based on aerodynamic panel codes (mostly DLM [59,60]) that idealise wing and empennage as lifting surfaces; if considered, the aircraft fuselage is treated with dedicated elements for non-lifting bodies. The lower computational cost compensates for the lower fidelity of the model; however, the reliability may be increased by correcting the aerodynamic influence coefficients with higher-fidelity data (which are generally not available at the preliminary stage of an MDO process though). Considering small perturbations of an unsteady irrotational flow, the BEM proposed by Morino [55,56] is hereby adopted and consists of an integral representation of the velocity potential at any point of the computational domain in terms of the values on the surface surrounding the aircraft body and wake; the principle of superposition applies. The linear equation of acoustic waves propagation governs the unsteady flow for the case of isentropic compressible fluid and accounts for the finite speed of sound [61], whereas Laplace's equation is resumed for the case of steady flow by taking advantage of Prandtl-Glauert's transformation [62]. Green's function representing a unit-impulsive point source, the theoretical formulation of the boundary value problem for the perturbation potential is then based on Green's formula [33], with (Neumann-type, generally instationary) non-penetration boundary conditions on the aircraft surface and wake as well as (Dirichlet-type, generally stationary) unperturbed asymptotic conditions far from the latter. Bernoulli's theorem is linearised to calculate the pressure coefficient and Kutta's condition is imposed at the trailing edge of lifting surfaces, where the wake detaches and is shed back with the reference airspeed, trailing circulation variations without sustaining any pressure load [37]; its trajectory then represents a streamline at any time and may generally become part of a nonlinear iterative solution where roll-up occurs due to downwash effects, if it is not prescribed a priori (as is in fact typical for most practical applications [53,63], especially when characterised by moderately unsteady flow). Thus, Morino's BEM is able to deal with arbitrary 3D geometries in a unified manner, reducing the effort of the abstraction process and easing the integration with complex structural models. Although the original theory is valid for arbitrary motion in time or frequency domain, only harmonic oscillations in the frequency domain are here considered as suitable for flutter analysis.
The current implementation of Morino's BEM is described in previous works [64] and is equivalent to an appropriate combination of doublets (for lifting bodies and wake) and point sources/sinks (for thickness effects and non-lifting bodies). The 3D geometry is approximated with flat quadrilateral panels that follow the local wing surface and the aerodynamic potential is assumed constant over them, with an analytic expression for their mutual induction. Note that similar flow conditions can be modelled by either vortex or doublet distributions and a quadrilateral doublet element is equivalent to a vortex ring placed at the panel edges [65]; higher-order methods have also been formulated in order to improve accuracy and computational efficiency for complex geometries and configurations, but they are more prone to becoming ill-posed and typically not necessary as long as enough quadrilateral panels refine the aerodynamic grid for lower-order methods to capture high pressure gradients while still staying away from ill-conditioning [66]. Through the wake being assumed as flat, rigid and parallel to the free-stream velocity in order to perform rigorous comparisons with Nastran's DLM and Theodorsen's theory without loss of generality (in the absence of aircraft fuselage and empennages [67]), a linear system is obtained: where frequency-dependent T σ and T ϑ are the matrices of aerodynamic influence coefficients (AIC), whereas σ k and ϑ k are the aerodynamic potential and normal wash associated with the k-th mode shape of the structural displacement [68], having mapped the latter onto the aerodynamic grid with an effective implementation of the infinite plate spline (IPS) method [69]. The elements Q hk of the generalised aerodynamic forces (GAF) matrix Q(s) are then calculated from the work done by the aerodynamic pressure C p due to the k-th mode on the displacement defined by the h-th mode φ h for a prescribed range of reduced frequencies [40], namely: where n i and A i are the normal vector and area of the i-th aerodynamic panel, respectively. The NACA0002 aerofoil [26] is suitably adopted to obtain a baseline 3D representation of Loring's wing surface (see Figure 2); this particular choice of thin symmetric aerofoil is consistent with the original work with no loss of generality, as the parametric sensitivity of the aeroelastic stability boundary with respect to aerofoil thickness was already found to be poor for small perturbations of subsonic potential flow (i.e., in the absence of strong shock waves and/or significant flow separation) [70].

Aeroelastic Model
After obtaining the natural vibration modes, the GAF are computed with Morino's BEM, and the continuation method [50] is used to solve the aeroelastic stability eigenproblem and trace the root locus. In order to avoid its costly re-computation during the stability analysis, the GAF matrix is approximated by a rational expression in which the non-linear dependency on the complex reduced frequency s appears explicitly. Here the matrix fraction approach (MFA) [71] is used and preferred to the rational function approximation (RFA) [72], since the former exhibits higher accuracy than the latter for the same number of poles [73]; yet, both methods provide the analytical continuation of the GAF for complex reduced frequency, increasing the accuracy of damped solutions with respect to the p − k method. According to MFA, the GAF matrix is expressed as fraction of matrix polynomials: and the accuracy of the approximation increases with increasing the number of poles M, which is equal to the size of the state-space system [74]; in particular, three poles have been used in this proof-of-concept work. The approximation matrices D i and N i are both obtained by solving a least-square problem in which the distance from the GAF samples Q(ik) is minimised.
The continuation method hereby adopted provides a straightforward and efficient tracking technique without using any correlation function, such as the modal assurance criterion (MAC) [75]. For this reason, the method is rather insensitive to the number of poles used for the aerodynamic finite-state approximation and is able to distinguish between actual and artificial aerodynamic states by construction. In the continuation method, the aeroelastic equations are differentiated with respect to the free-stream airspeed, thereby resulting in a system of ODEs which are solved with a predictor-corrector integration schema starting from an initial airspeed value for which a true solution is known (e.g., the natural vibration frequencies for U = 0).
For validation purposes, flutter analysis is also carried out with Nastran's DLM directly: a flat plate aligned with the free-stream is used as lifting surface and the aerodynamic panels are defined with CAERO1 and PAERO1 cards; the distribution and number of such panels follow established guidelines and convergence studies (see Appendix B). The mode shapes are mapped onto the aerodynamic grid with the IPS method [76] using the SPLINE1 card and the GAF matrix is then generated for sixteen reduced frequencies specified in the MKAERO1 card, suitably ranging from k = 0 to k = 1 (with logarithmic-like spacing) as per literature studies and common practice [53,63]; a cubic spline is then exploited to interpolate the GAF therein. For the dynamic aeroelastic stability analysis, the p − k method as available in the FLUTTER card is used: the aerodynamic damping is approximated as the imaginary part of the GAF matrix computed for (undamped) harmonic motion, limiting results' accuracy to the case of lightly damped aeroelastic solutions [12]. The p − k approximation of the aeroelastic eigenproblem [77] is then solved for all the free-stream speed values defined in the FLFACT card, according to the reference length and air density specified in the AERO card. After damping ratios are obtained in the given speed range for all modes, the (undamped) flutter point is determined by linear interpolation.

Aero-Structural Parametric Derivatives
To calculate the sensitivity of the flutter boundary with respect to any aero-structural parameter, the derivatives of the aeroelastic matrices are also necessary in the first place. In particular, the derivatives of structural eigenvalues and eigenvectors with respect to structural parameters, such as wing material properties (e.g., density or elastic modulus) and element properties (e.g., moments of inertia or distance between elastic and inertia axis), are obtained with Nastran SOL200. The design variable label and initial value are defined with the DESVAR card, which is connected to the relative bulk data entry by DVMREL1 or DVPREL1 cards for material or element properties, respectively; derivatives are hence computed for the structural responses defined by DRESP1 cards, which are the first n structural eigenvalues and eigenvectors.
The GAF matrix is differentiated with respect to both reduced frequency and design variables: according to MFA, the finite-state approximation allows the analytic differentiation with respect to the complex reduced frequency s, whereas the method presented in previous works [73] was used for the derivatives with respect to the design variables p. In particular, by defining the functions P σ k , p and R σ k , p as: a sub-differentiation process is set up, and depending on the number of design variables and mode shapes to be considered [11], the derivatives may be obtained by either a direct approach: or an adjoint approach: In computing the partial derivatives, it is important to note that the imaginary component is introduced only by the normal wash ϑ and the wake's influence coefficients inside matrix T σ ; therefore, the partial derivatives of the steady (real) contribution to the GAF are computed via complex step [78] and the partial derivatives of the unsteady (complex) contribution are then analytically assembled.

Results and Discussion
Loring's uniform slender thin wing [25] is hereby considered, as both experimental wind-tunnel data and numerical results assuming two-dimensional incompressible potential flow are available and all results can hence be explained from both physical and mathematical standpoints; moreover, this fundamental benchmark embeds the full complexity of the aeroelastic problem without introducing detrimental uncertainties. The wing's chord is 0.305 m and the semi-span is 2.057 m, giving an aspect ratio 13.5. The inertial axis lays at 42.3% of the chord, with mass 8.05 kg/m and mass moment of inertia 0.0471 kg·m. The wing root is clamped at the elastic axis, with flexural stiffness 677.3 N·m 2 and torsional stiffness 1018.9 N·m 2 placed at 30% of the chord.
The coupled natural vibration frequencies of the first bending and torsion modes were observed at 1.29 and 18.1 Hz, respectively; that of the second bending mode was detected in between at 7.7 Hz. With a fluid density 1.11 kg/m 3 , the flutter speed and frequency were experienced at 90.3 m/s and 10.2 Hz, which give a reduced frequency around 0.11; the subsonic flow may then be considered as moderately unsteady. A Euler-Bernoulli beam model calculated the first three coupled vibration frequencies as 1.29, 7.65 and 17.98 Hz; once coupled with two-dimensional unsteady aerodynamics for a flat plate using standard strip theory, flutter was predicted at 90.7 m/s and 9.2 Hz with good accuracy.
In order to elaborate on the literature results and visualise the flutter mechanism, the aeroelastic stability analysis of Loring's wings is first performed with the same assumptions as in the original publication and the approach here described in the Problem Formulation; the p − k method (with Theodorsen's exact function [35]) has consistently been used and results have been cross-verified against the state-space representation (with a common two-term RFA of Theodorsen's function) [21]. By employing the first two bending modes and the first torsion mode, the respective coupled vibration frequencies are calculated as 1.21, 7.59 and 17.91 Hz; once still coupled with two-dimensional unsteady aerodynamics for a flat plate using standard strip theory, flutter is consistently predicted at 91.15 m/s and 9.2 Hz (which is indeed an excellent approximation, as confirmed by a modal convergence study up to the first three bending and torsion modes). The evolution of the aeroelastic system's eigenvalues is presented in Figure 3 and confirms static divergence to arise well beyond dynamic flutter (i.e., U d >> U f ). Due to the remarkable agreement between measurements and simulations, the assumptions of slender beam structure and two-dimensional incompressible potential flow in the latter are hence deemed justified. Note that the instability mechanism involves second bending and first torsion modes directly, but the first bending mode is indirectly essential for their coupling to occur before that between first bending and first torsion modes (see Appendix A).

Aeroelastic Analyses
The aeroelastic stability of Loring's wing is investigated and compared using the proposed lower and higher-fidelity models, while still assuming incompressible flow. The analyses encompass parametric variations of the aero-structural properties, quantifying their effects on the aeroelastic stability boundary and critically assessing the differences in the multifidelity results from both theoretical and computational perspectives, for possible practical applications in airplane design and optimisation adopting hybrid strategies. Figure 4 shows the aeroelastic stability analysis from the higher-fidelity model, focusing on the same flutter mechanism as found earlier: the first torsion mode becoming unstable and extracting energy from the airflow through the coupled second bending mode, flutter is found at 94.77 m/s and 10.04 Hz in excellent agreement with the experimental results and exhibits a slightly higher vibration frequency than the theoretical one. The aeroelastic stability analysis performed with Nastran and its embedded DLM for the flat wing is also presented and confirms the higher-fidelity results based on Morino's BEM with NACA0002 aerofoil, flutter being found at 94.44 m/s and 10.38 Hz. Note that both sets of code share the same structural model, and the coupled natural vibration frequencies in the void are 1.21, 7.55 and 21.03 Hz for the first three bending modes, whereas it is 17.88 Hz for the first torsion mode; these first four modes have been used in all analyses (see Appendix B). As expected from the wing aspect ratio being relatively large, three-dimensional downwash effects are actually moderate and the (beneficial) steady ones on the (attenuated) airload distribution are roughly compensated by the (detrimental) unsteady ones on the (accelerated) airload evolution [21]; thus, the flutter speed is just slightly higher than that predicted assuming unsteady two-dimensional flow. Still focusing on the flutter mechanism, Figure 5 then shows the aeroelastic stability analysis from the lower-fidelity model: despite its inherent simplicity, the latter predicts flutter quite accurately at 92.1 m/s and 9.09 Hz, thereby proving the effectiveness of such an idealisation (analogous stability diagrams have indeed been obtained by retaining only steady aerodynamic terms in the higher-fidelity model, for cross-verification purposes). The coupled natural vibration frequencies in the void are 1.21 and 7.59 Hz for the two bending modes, and 17.91 Hz for the torsion mode, in excellent agreement with the higher-fidelity ones. The aeroelastic response is then metastable until the first instability occurs, as no damping is provided by either the elastic structure or static aerodynamics; to this respect, it is worth mentioning that the lift-derivative correction for the steady three-dimensional downwash (which tends to increase the flutter speed, by providing lower airload) incidentally compensates the lack of aerodynamic lag (which tends to decrease the flutter speed, by neglecting wake inflow). Note that, although the flutter reduced-frequency is rather high for steady aerodynamics to be rigorously applicable, the latter was still adopted in order to exacerbate the difference between lower and higher-fidelity models and stress the multilevel approach for more general conclusions.

Sensitivity Study
Following previous studies on the aeroelastic stability boundary of slender cantilever wings and its sensitivity [27,28], linear FSI models give nonlinear flutter trends as relevant aero-structural parameters are perturbed. The reciprocal positions of aerodynamic, elastic and inertial axes being often constrained by the available space inside the wing as well as the chosen structural layout and systems arrangement, the flutter point's sensitivity to the wing's structural properties is individually explored by changing its material density ρ m (which scales m and µ) and elastic modulus E (which scales EI and GJ), whereas varying the wing's semi-span l (which modifies the lift coefficient derivative C l/α ) alters its geometry and affects both structural and aerodynamic properties at the same time.
From Figure 6, it can be seen that changing the wing's material density alters all natural frequencies through the inertial properties and hence the flutter frequency, but has marginal/no influence on the flutter speed; all symbols give the actual nonlinear percentage variation of the flutter point, whereas all lines with corresponding colour draw its normalised linear prediction. In particular, a 1% increment in the material density causes about a 0.5% decrement in the flutter frequency; the negligible variation of the flutter speed is confirmed by previous works [27] as being due to the (small) beneficial effect of an increase in sectional mass being compensated by an (almost) equal and opposite detrimental effect of an increase in mass moment of inertia. This outcome may indeed help while minimising the aeroplane weight, as reducing the material density does not significantly affect the present flutter boundary due to wing bending-torsion instability (especially in the absence of lumped masses). However, it shall be recalled that varying the wing inertia might cause other types of aero-structural instabilities to arise at the aircraft level due to modal coalescence, and monitoring the variation of the flutter frequency is then important to preventing potential resonance. It is also worth mentioning that the small higher-fidelity effect on the flutter speed is due to the unsteady airload being frequency dependent through the lift-deficiency function and the related aerodynamic lag. As per Figure 7, changing the wing's material elastic modulus alters all natural frequencies through the stiffness properties and hence the flutter speed and frequency; in particular, a 1% increment in the elastic modulus induces about a 0.5% increment in all the latter, leaving the flutter reduced frequency practically unchanged. This outcome is quantitatively confirmed by previous works [27] as being due to the detrimental effect of an increase in flexural stiffness being much lower than the beneficial effect of an increase in torsional stiffness. The striking agreement between lower and higher-fidelity results is mostly driven by the respective structural models being equivalent (as seen from the close agreement between the natural vibration frequencies), while the different aerodynamics play a minor role (as the airflow is moderately unsteady). The quasi-linear trend of the percentage variations reveals a rather large confidence interval for the normalised sensitivity; however, the dimensional counterpart of the latter to be used by optimisation routines follows a highly nonlinear pattern (note that the explicit lower-fidelity expression given in Appendix A for the torsional static divergence speed provides a straightforward theoretical check). Finally, Figure 8 shows that changing the wing semispan has a large effect on the stability boundary through the stiffness properties as well as the aerodynamic loads and hence on the flutter speed and frequency. In particular, a 1% increment in the wing semispan induces about a 0.8% decrement in the critical speeds and about a 1.6% decrement in the flutter frequency; note that analogous trends for the variations and orders of magnitude for the sensitivities were obtained in previous studies on similar slender wings [28], as a qualitative means of validation. Lower and higherfidelity results are still in very good agreement and differences are mainly due to unsteady downwash effect, as the aspect ratio of the wing changes considerably with the span and so does the sectional airload. For the sake of completeness, it is worth mentioning that variations and sensitivities of the divergence speed (see Appendix A) with respect to the same aero-structural parameters considered above were consistently found to have trends and orders of magnitude very close to those pertaining the flutter speed. In particular, it is observed that changing the wing's geometry has a much larger impact on the flutter boundary than changing the material properties (especially density); this is particularly true for the flutter frequency, the variations of which exhibit a more significant nonlinear behaviour and local curvature.
Due to the sufficient mutual separation of the natural vibration modes, it is also worth stressing that the instability mechanism did not change across the parametric variations, and the derivatives of the flutter stability boundary have accurately been calculated by the lower-fidelity model at a (marginal) fraction of the computational complexity and costs. In particular, semi-analytical solutions drastically reduced the latter, being almost instantaneous and demanding minimal pre-and post-processing (if any at all) while granting an enhanced theoretical understanding. Thus, the healthy combination of lower and higher-fidelity models enables efficient multidisciplinary exploration of a large design variable space for innovative aeroplane concepts and configurations; further numerical savings may still be obtained by exploiting reliable surrogate models for the higher-fidelity solutions, with an effective synthesis of the underlying complexity [20].

Conclusions
Within the context of aircraft multidisciplinary design and optimisation, a comparative sensitivity study for the bending-torsion flutter instability of flexible aircraft wings in subsonic flow has been presented. Analytical models and numerical tools with different complexities and fidelities have been used, in view of possible practical applications exploiting multilevel approaches within the conceptual and preliminary MDO phases. Parametric studies have been performed where the effects of varying the wing's aero-structural properties on the aeroelastic stability boundary have been quantified and critically assessed from both theoretical and computational perspectives. When the natural vibration modes of the wing are well separated from all other natural vibration modes of the aircraft and the aeroelastic instability mechanism does not change in nature, an efficient hybrid strategy is then recommended where the flutter analysis is performed using higher-fidelity approaches, whereas the sensitivity analysis of the flutter boundary is performed using lower-fidelity approaches, thereby improving theoretical understanding and reducing computational costs while retaining accuracy. Future works are encouraged to investigate additional effects (e.g., wing sweep and flow compressibility) and increased complexity (e.g., presence of a control surface) or perform the full MDO of a flexible aircraft wing, effectively exploiting the proposed multifidelity strategy.

Appendix A. Aeroelastic Stability of the Typical Section with Steady Aerodynamics
Providing full control on the results with respect to the specific assumptions of the methods and tools employed for the analysis, theoretical formulations grant a clear and complete overview of the problem which is essential for any engineering application. Although inherently limited in their general capabilities, analytical solutions then offer a wealth of qualitative information and quantitative details as well as fundamental insights and rigorous validation for both educational and practical purposes.
Considering a single mode for the wing bending and torsion, respectively, the aeroelastic equations for the Typical Section with steady aerodynamics read [79]: where f ≡ 1 for two-dimensional aerofoils [12], with φ h ≡ 1 and φ θ ≡ 1. By neglecting the aerodynamic load, the coupled bending and torsion natural vibration frequencies are explicitly obtained as: where the uncoupled counterparts are resumed whenever inertial and elastic axes coincide (i.e., with x CG = 0). Otherwise, the characteristic equation for the flutter determinant provides with the metastable boundary and gives the flutter frequency as: and setting the inner radical discriminant equal to zero then gives the flutter speed as: where the sign before the inner radical shall provide with the smallest positive real value; finally, note that setting both ω f and P 0 equal to zero provides with the torsional static divergence speed. In particular, when the aerodynamic centre is ahead of the elastic axis (i.e., with x AC < 0), the latter explicitly reads: regardless the bending stiffness; when the elastic axis is ahead of the inertial axis (i.e., with x CG > 0), the flutter condition is given as: where all coefficients are analytical functions of the aero-structural parameters, namely: The formulas above were implemented and compared with literature results for the stability boundary of a Typical Section assuming standard strip theory (SST) with steady aerodynamics [12,79]: exact agreement was always found and provided rigorous validation, also with respect to the relative root locus. The same formulas have then been used to investigate the stability boundary of Loring's wing Typical Section when either the first or the second bending mode is coupled with the first torsion mode, employing tuned strip theory (TST) with steady aerodynamics [21]. Without accounting for the modal cross-projection (i.e., f ≡ 1 for a "pitch & plunge" apparatus), flutter is calculated at 106.5 m/s and 4.32 Hz in the former case whereas at 73.9 m/s and 11.28 Hz in the latter case; when the modal cross-projection is considered (i.e., f < 1 for a slender beam), flutter is calculated at 109.7 m/s and 4.28 Hz in the former case whereas at 139.2 m/s and 9.60 Hz in the latter case. The static divergence speed pertains the torsion mode only and is found at 210.2 m/s, regardless the bending modes and their cross-projections.
A direct comparison with the low-fidelity results hence reveals that the interaction between bending and torsion modes as well as their cross-projections are essential to reproduce the correct behaviour of the flutter mechanism, where the coalescence between first bending mode and torsion mode drives that between second bending mode and the latter (which have closer natural vibration frequencies) to occur earlier.

Appendix B. Higher-Fidelity Model Results Convergence Study
The number of elements for the structural and aerodynamic models was defined according to rigorous convergence studies. Following Nastran's best practice [80] with k max = 1 the maximum reduced frequency employed in the aeroelastic analysis, fifteen DLM panels have uniformly been distributed along the wing chord. When adopting Morino's BEM, thirty aerodynamic panels were symmetrically placed along both upper and lower aerofoil surfaces according to the convergence study in Figure A1 (left), with a suitable refinement around the leading edge in order to capture high pressure gradients. The number of aerodynamic panels in the span-wise direction was determined according to another convergence study shown in Figure A1 (right) and all higher-fidelity results presented in this work were obtained with fifteen panels strips uniformly distributed along the wing span, since sufficient to grant a relative error below 1% for both flutter speed and flutter frequency. Following convergence studies in previous works [37,64], the wake extends for a hundred chords behind the wing and is modelled with a hundred rows of shed panels, the length of which is cubically increased with increasing their distance from the trailing edge. Figure A2 (left) shows the indicial lift-deficiency function Φ(τ) (which represents the equivalent of the complex lift-deficiency function C T in the reduced-time domain τ) from a unit step in angle of attack, whereas Figure A2 (right) depicts the normalised distribution of the steady circulation κ(y) (which represents the ratio between the sectional circulation in three-and two-dimensional flow within the framework of modified strip theory, MST [13]); note that the unsteady lift development approaches quasi-steady theory (QST) asymptotically. These results also justify the assumption of two-dimensional flow for the lower-fidelity model, within the framework of TST [21] (to this respect, it is worth mentioning the higher-fidelity lift coefficient derivative of the wing is obtained as C L/α = 5.25, in excellent agreement with the lower-fidelity analytical estimation C L/α = 5.21).   Figure A3 (left) then shows the evolution of the natural vibration frequencies percentage error with varying the number of FEM nodes of the structural beam; the value obtained with the most dense grid is used as reference to compute the plotted error. The first four natural vibration frequencies exhibit a good convergence behaviour, with the third one (i.e., the first torsional mode) almost insensitive to the number of nodes. All higher-fidelity results presented in this work were hence obtained with 20 nodes uniformly distributed along the wing span, since sufficient to grant a relative error below 2%. Finally, Figure A3 (right) shows that employing the first four natural vibration modes granted convergence of both flutter speed and flutter frequency.