Numerical Whirl–Flutter analysis of a tiltrotor semi-span wind tunnel model

This work presents the modeling and preliminary Whirl–Flutter stability results achieved within the Advanced Testbed for TILtrotor Aeroelastics (ATTILA) CleanSky2 project. The project addresses the design, manufacturing, and testing of a semi-span wind-tunnel model of the Next Generation Civil TiltRotor. The preliminary multibody models developed in support of the wind-tunnel testbed design are described, illustrating the modeling technique of each subcomponent of the model, namely the wing, the rotor, the blades, and the yoke. The methodologies used to analyze the stability of systems subjected to periodic aerodynamic excitation when the problem is modeled using full-featured multibody solvers are presented in support of Whirl–Flutter identification during wind-tunnel testing.


Introduction
Tilting rotor aircraft, or tiltrotor, with their outstanding capability of taking off and landing vertically like helicopters and, at the same time, of achieving high speeds-up to twice that of helicopters-in forward flight when operating like conventional turboprop aircraft, represent one of the few successful examples of advanced vertical lift configurations.
After a long development phase that encompassed few experimental aircraft that successfully made it to flightnoticeably the Transcendental Model 1-G (a tilting rotor concept), which flew about 100 h without ever completing a full conversion [26] and the Bell XV-3 (a tilting rotor configuration) in the 1950s [12], and the Bell XV-15 (a tilting rotor/nacelle configuration) in the 1980s [26], unlike other less fortunate technology demonstrators-the concept finally proved its soundness with the Bell-Boeing V-22 (a tilting rotor/nacelle configuration) [46], a military aircraft that after an almost 25 years long development phase, has been operated by the US Marine Corps with an increasing reward for the last 15 years, and is entering service for carrier onboard delivery (COD) in the US Navy. Remaining in the military arena, the Bell V-280 Valor (a tilting rotor configuration) after a test campaign started in 2017 is currently under evaluation by the US Army for its Future Long-Range Assault Aircraft (FLRAA) program [13]. However, with the then Bell-Agusta and now Leonardo AW609 (a tilting rotor/ nacelle configuration) [36] about to become operational after long and thorough development [15], the tiltrotor design appears to be mature enough to enter also the civil air transport market [1].
Nevertheless, tiltrotor design remains a rather challenging engineering task, considering the various operating conditions and multipurpose missions that are expected to be accomplished by this complex type of aircraft. In particular, the problem of assessing Whirl-Flutter stability limits is at the same time fundamental and challenging. Whirl-Flutter is an aeroelastic stability phenomenon that affects both turboprop and tiltrotor aircraft [6]. When a rotor mounted on a flexible structure rotates, the normal vibration modes associated with the supporting structure may interact with the precession motion of the rotor.
When its motion is perturbed, each point on the rotor axis of rotation draws paths about its reference position. This motion changes the way each rotor blade is affected by the incoming airspeed, correspondingly altering the overall aerodynamic loads. At the verge of whirl flutter, when this phenomenon is triggered, perturbations result in a periodic orbit. The resulting forces can lead to the divergence of the system response, in what represents a truly aeroelastic instability [45].
Nowadays, although the key aspects of the phenomenon in tiltrotor aircraft are understood, the capability to predict it is still limited, despite the availability of sophisticated aeroelastic analysis tools. The difficulty lies in its dependence on several factors, including the geometrical design, the structural properties, the dynamics of the actuators, etc. which can all contribute to the Whirl-Flutter characteristics in ways that are not always intuitive. The problem may be exacerbated when unconventional configurations are considered, as is the case of propeller-and rotor-driven advanced air mobility concepts [22].
Aeroelastic testing, especially when it concerns stability and flutter, can be extremely dangerous. For this reason, wind-tunnel testing is needed to understand phenomena and increase the level of trust in predictions in support of their mandatory verification in flight. As a consequence, the possibility to validate numerical predictions using experimental data is of paramount importance from an industry standpoint.
Wind-tunnel testing of tiltrotor configurations, and its use for validation and correlation with numerical predictions, has been widely used for tiltrotor development. Due to the need to aeroelastically scale also the wing portion of the model, of extreme miniaturization of the rotor-consider for example the need to reproduce in detail rather complex blade root attachments and pitch control mechanisms-aeroelastic wind-tunnel models of tiltrotor aircraft can be even more complex than helicopter ones.
Several wind-tunnel models have been developed to explore the concept of Whirl-Flutter in prop-and tiltrotor configurations. However, specific wind-tunnel testing of configurations that were intended to be tested in flight only dates back to the early development of the V-22 when a 1/5scale semi-span model was designed and manufactured by Bell and used for wind-tunnel tests and correlation with simulation capabilities available at that time [42]. Eventually, this model evolved in the Wing and Rotor Aeroelastic Test System (WRATS), which has been in use for several years to explore various aspects of tiltrotor aeroelasticity, including aeroelastic tailoring of the wing [10,35], several parametric studies (e.g., in [40]), active control for flutter suppression [17,25], an interesting four-blade, articulated soft-in-plane configuration for hub load reduction [33,34], and several numerical studies of stiff-and soft-inplane configurations successfully correlated with experiments [29,30,48,49].
Another effort towards the wind-tunnel testing of the aeroelastic stability of tiltrotors took place in Europe, within the EUROFAR program [21]. Subsequent efforts aimed at the development of an original European tiltrotor concept, ERICA [32]. This concept is characterized by the tilting of the outer portions of the wing, affected by the downwash of the rotor. These efforts resulted among the others in the project ADYN [5,24], aimed at studying the aeroelasticity and the noise characteristics of the concept.
At about the same time, yet another effort is ongoing in Europe, within the CleanSky2 initiative, within the project ATTILA, to design, manufacture and test a tiltrotor windtunnel model in support of the Next-Generation Civil Tilt Rotor-Technology Demonstrator (NGCTR-TD), a fixedengine tilting rotor configuration. This paper illustrates the preliminary modeling efforts of the wind-tunnel testbed, as initially presented in [8,9]. Preliminary studies were also conducted in [18]. Numerical models of the reference testbed have been developed using FLIGHTLAB 1 and MBDyn 2 [28], respectively, by collaborating research teams at the Royal Netherlands Aerospace Centre (NLR) and Politecnico di Milano. To the authors' knowledge, a direct correlation of Whirl-Flutter results from the two solvers has not been published elsewhere.
These numerical models are developed in support of the design of the physical wind-tunnel testbed. Pending experimental characterization of the ATTILA testbed dynamics and Whirl-Flutter stability properties, code-to-code verification results are here presented, compared with data obtained from a reference CAMRAD II model of the full-scale aircraft. In addition to component-level models verification, preliminary Whirl-Flutter predictions are presented. In this context, the Matrix Pencil Estimation (MPE) and Periodic Operational Modal Analysis (POMA) eigensolution identification methods are used and compared for the first time, to the authors' knowledge.

FLIGHTLAB
FLIGHTLAB is a state-of-the-art multibody, componentbased, selective fidelity modeling and analysis software package. It supports modeling and simulation of rotorcraft, fixed-wing aircraft, compound aircraft, helicopters, multi-copters, drones, flying cars and experimental aircraft configurations. Rotorcraft and other aircraft models can be developed to fit their application with the desired level of fidelity. The numerical models can be used for engineering analysis, real time simulation, or both. The development system is also used to generate run-time models for real time applications. The key capabilities of the software include: -Multiple bodies, multibody dynamics -Nonlinear unsteady aerodynamics -Flight dynamics and real time simulation (including piloted, full flight simulators) -Flight performance, stability, controllability, and handling qualities -Aeroelastic stability, vibration, and loads -Aircraft systems analysis and hardware-in-the-loop (HIL) simulation -Coupling with external programs, including CFD and Matlab/Simulink The Whirl-Flutter analysis is usually performed through direct linearization of the problem about a steady trim solution, using central difference perturbation of the firstand second-order degrees of freedom. This technique permits to identify the wing-pylon, drive system and rotor modes, providing valuable information about the underlying dynamics. A nonlinear transient analysis resulting from the application of appropriate excitations, replicating the expected wind-tunnel test procedure, has also been implemented to verify the results of the linearized stability analysis, investigate various experimental excitation strategies, and determine the required excitation magnitude and the associated structural loads.

MBDyn
MBDyn is a general-purpose multibody solver, developed at the Aerospace Science and Technology department of Politecnico di Milano and distributed as free software.
MBDyn automatically writes and solves the equations of motion of a system of entities possessing degrees of freedom (nodes) connected through algebraic constraints and subjected to internal and external loads. Constraint equations are explicitly accounted for, following a redundant coordinate set approach. Thus, the resulting system of Differential-Algebraic Equations (DAE) takes the form where x is the vector of the kinematic unknowns, p that of the momentum unknowns, that of the algebraic Lagrangian multipliers, is a configuration-and time-dependent inertia matrix, f i , f e are arbitrary internal and external forces, (x) is the vector of the (usually nonlinear) algebraic equations that express kinematic (holonomic) constraints, and ∕ is the Jacobian matrix of the constraints with respect to the kinematic unknowns. Each node instantiates the corresponding balance equations (1b), while only nodes with associated inertia properties instantiate the related momenta definitions (1a). Additional states associated with scalar fields (namely, hydraulic pressure, temperature, electric potential) and thus the corresponding balance equations, can be taken into account through dedicated sets of nodes.
Elements are responsible for the contributions to the balance equations through (visco-)elastic, internal forces i , possibly state-dependent external force fields e (e.g., aerodynamic forces), and reaction forces f c = T ∕x , introduced using the Lagrange multipliers and the Jacobian matrix of the algebraic constraint equations in Eq. (1c).
The DAE system can be integrated using several different A/L stable integration methods, among which is an original multistep method with tunable algorithmic dissipation, specifically designed for the class of problems usually solved with MBDyn.
The stability analysis of a nonlinear multibody model, formulated in an absolute reference frame, is critical since the problem may be time-variant. In principle, since in the configuration under analysis, the flow is essentially axial, this specific problem could be reduced to time-invariant through the so-called multiblade coordinates, with the possibility to determine a steady solution, when asymptotically stable, for its linearization and analysis using methods for linear, time-invariant (LTI) problems. However, owing to the flexibility of the wing and the opportunity to trim the aircraft at non-zero angles of attack, at the trim point the flow might not be perfectly axial; furthermore, aerodynamic interferences between the rotor and the wing may break the symmetry of the axial flow, introducing 1/

3
rev periodicity in the motion of each blade, and N b /rev periodicity in the overall response of the system.
It is thus convenient to use methods specifically designed to address problems of this kind. A typical approach, based on Floquet's theory as introduced in the stability analysis of rotorcraft by Peters & Hohenemser in 1971 [38], evaluates the stability of linear, time-periodic (LTP) problems. It requires the integration of the state transition matrix of the problem over a revolution, which is called the monodromy matrix; this is often impractical when using multibody solvers, which formulate the problem as a large number of nonlinear DAEs, as discussed in [4]. A generalization of Floquet's approach to nonlinear, aperiodic problems, possibly subjected to random excitation, has been recently proposed in [52], based on Lyapunov's theory of characteristic exponents (LCE); however, the applicability of this method to systems of DAE is still debated [27], and hardly practical from a computational viewpoint for large systems. An interesting approach, based on the proper orthogonal decomposition (POD) of snapshots of large sets of coordinates of the problem, taken during free-response transients at each rotor revolution has been proposed in [43] and successfully used with MBDyn. The equivalents of Floquet's characteristic coefficients are obtained from the eigenanalysis of a matrix that corresponds to Floquet's monodromy matrix, resulting in a least-squares fitting of the snapshots.
All the previously mentioned methods can, to a different extent, evaluate the "contraction" characteristic of the free response of the system (e.g., the real part of the eigenvalues, when defined, or the characteristic coefficients of Floquet's theory, or the characteristic exponents of Lyapunov's theory), but some are unable to determine its periodicity, when the free response is oscillatory [4,43,52], while others suffer from limitations in the interpretation of the imaginary part of the eigenvalues of the monodromy matrix [38], although attempts have been made to overcome them (see for example [39]).
In this work, two somewhat complementary approaches are followed: -Matrix Pencil Estimation (MPE), used to retrieve the main wing modes; -Periodic Operational Modal analysis (POMA), used to take into account the periodicity of the system to extract not only the wing's normal modes, which represent fundamental harmonics of the system's free response, but also the corresponding superharmonics, namely the fundamental harmonic plus or minus integer multiples of the rotor angular velocity, resulting from the interaction with the rotor itself.
The workflow of the two methods is sketched in Fig. 1, from the MBDyn simulation results to frequency and damping, identified by both methods, and participation factors and mode shapes, only reconstructed by POMA. The details of each process are discussed in the following sections.

Matrix pencil estimation
The Matrix Pencil Estimation method (MPE) was designed by Hua et al. [20] to estimate parameters of exponentially damped or undamped sinusoids in the presence of noise. A multiple-input algorithm was subsequently proposed by Favale et al. [41], who applied the methodology also to more complex problems, such as tiltrotor Whirl-Flutter analysis. The method can estimate the modal parameters of a system from its free responses to external input. Starting from a multibody simulation of the complete semi-span tiltrotor wind-tunnel model, excited using a sinusoidal input at the swashplate, all the available strain values in the beam elements and accelerations of the nodes are extracted. A selection process is then performed to ensure reliable identification and select the most meaningful data. After simple data pre-processing, consisting of filtering and re-sampling, the identification algorithm is applied to compute the system poles. A further postprocessing step involving the creation of stabilization diagrams, as suggested in [41], is finally executed to evaluate the quality of the identification process and retrieve the desired results, such as frequency and damping ratio.
The MPE method can estimate the modal parameters of a system starting from its free response to external finite input. A general response signal y(t) will be composed of two contributions: the system free decay x(t) and an unknown noise n(t) . In discrete-time, the signal can be written as: where N is the total number of samples obtained using a sampling frequency f = 1∕dt . According to the above hypothesis, the response signal of the system in discrete time can be approximated as a sum of complex exponentials: where M is the order of the system and z i is equal to: The free-response of the system can be described by its modal parameters: the complex amplitude R i and poles s i , with frequencies i and damping ratios i . The aim is then to find these parameters from the recorded response of the system. Starting from the noiseless signal x(k) a rectangular Hankel matrix with dimensions (N − L) × (L + 1) can be created: where L is called the pencil parameter, which should be smaller than (N − 1) . Two matrices will be then created, excluding the first and the last column of X: The matrices 1 and 2 can be re-written as: where the matrices R , Z 0 , Z 1 , and Z 2 are functions of the complex amplitudes R i and system's poles z i : It can be demonstrated that M is the rank of the matrix pencil if M ≤ L ≤ N − M and the diagonal elements of 0 are the eigenvalues of the matrix pair 2 and 1 . Finding the system's poles reduces to solving the following eigenvalue problem:

3
The same procedure described here holds for signals disturbed by noise. Hua et al. [20] suggested that to obtain the best results the parameter L should be chosen as high as possible but limited between N/3 and N/2. After creating the block Hankel matrix of the collected responses, , a fundamental step is to perform its singular value decomposition. Then, a new set of filtered matrices can be obtained by truncating the matrices at the first dominant values. From their analysis, it is possible to determine, in first approximation, the order of the system M: when using very high signal-tonoise ratio data, the first M singular values can be of orders of magnitude greater than the others. However, when very noisy signals are considered, the singular values relative to effective system contribution do not significantly stand out among all the others. In this case, after having normalized the singular values, a possible approach is to compute their first derivative and identify the order corresponding to its maximum value. For simplicity, the previous discussion was performed considering a single input response. However, the algorithm can be easily extended to multiple input channels considering a matrix of responses instead of a column array.

Periodic operational modal analysis (POMA)
Rotorcraft modal identification is typically performed in nonrotating coordinates by applying a multiblade coordinate transformation (MBC). However, when dealing with real-world problems, several practical issues may arise due to, e.g., small anisotropy of the rotor blades, or slightly different axial position of the sensors on each blade. Furthermore, when performing a periodic stability analysis, richer insight into the system behavior can be retrieved: for each harmonic contribution captured with the MBC transformation approach, a (theoretically) infinite number of sub-and super-harmonics is obtained, each with its participation factor. The algorithm called Periodic Operational Modal Analysis (POMA), originally proposed by Wereley and Hall [54], is used in this work. It is based on the harmonic transfer function concept, namely the periodic counterpart of the frequency response function for LTI systems. Wereley defined a new fundamental signal space for periodic systems, containing the socalled exponentially modulated periodic signals (EMP). These have been defined as the complex Fourier series of a periodic signal of frequency p , modulated by a complex exponential signal: where s n = s + nj p , and u n are the Fourier coefficients of u(t). The harmonic transfer function can be defined as: The output can be expressed as: where Ω = 2 ∕T rev is the characteristic frequency of the system (i.e. the frequency of rotation of the rotor in the case of a tiltrotor), and Y( ) is the Fourier transform of y(t).
In theory, to consider all the harmonics that characterize a linear time-periodic (LTP) system, ( ) should be an infinite matrix. However, for most applications, a satisfactory approximation can be obtained with (often quite) a limited number of harmonics. Wereley also showed the expression of the harmonic transfer function in terms of modal parameters of the state transition matrix: where is the direct transmission term. The terms ̄ r,l and ̄ r,l are the Fourier coefficients obtained from the expansion of the following matrices: where r (t) and r (t) T are defined as: and matrix (t) can be computed as the product of (t) and a matrix containing the eigenvectors of the Floquet [14] factor .
The power spectrum of the output yy ( ) can be expressed as a function of the harmonic transfer function: Note that in this case, yy ( ) is the power spectrum of the output signals of the system which have been previously exponentially modulated to be able to apply the relation of Eq. (15).
Substituting Eq. (18) into (21) and setting to zero:  (20) where ( ) is a function of the input spectrum and the input characteristics of the system. Equation (23) shows that the autospectrum of the output can be approximated by a sum of modal contributions if ( ) is reasonably flat at least for the band of interest of a specific mode. A further simplification can be introduced with the assumption of suitably separated modes reducing Eq. (22) to: The expression in Eq. (24) can be compared to the one for the power spectrum of an LTI system: Equations (25) and (24) show that the peak related to any super-harmonic of a given mode can be viewed as the peak of a mode of a generic LTI system. Standard identification techniques can then be used to compute frequencies, damping factors, and modal shapes from the measured spectra.
Under the assumption that the input spectrum is relatively flat in the range of frequencies of the modes of interest, the ATTILA wing-pylon model in MBDyn was excited by superimposing a white noise disturbance to all the components of the wind speed. The resulting strains and displacements were pre-processed to retrieve zeromean signals, which were then exponentially modulated. Finally, a Stochastic Subspace Identification (SSI) [7,44] algorithm was applied to identify the frequency and damping from the strains to remove the rotor periodicity, and the mode shapes from the displacements. Using displacement signals allowed to extract mode shapes of the system and better visualize the system's deformation. This result is not easy to achieve for a complex multibody system. However, the procedure can become increasingly computationally expensive when a huge number of degrees of freedom are considered. Additionally, to retrieve accurate and reasonable mode shapes, a selection of the input data for the algorithm is not always feasible, increasing the odds of encountering clusters of modes, also considering the possible presence of noisy data in the identification process. When performing the analysis using strain measurements is easier to isolate modes' contributions, minimizing the possibility of overlapping modes by carefully selecting data among all available measurement channels. For the frequency and damping identification of wing, blade and (23) ( ) r,s,l,k =̄ r,l uu ( )̄ H s,k (24) yoke strain measurements were employed. Note that for the extraction of mode shapes and other modal parameters, both rotating and non-rotating measurements were used in the same identification process. As a consequence, the computed rotor-related frequencies are expressed in a rotating reference frame. In short, the steps of the identification algorithm are: -record the response y(t) of the system to a broadband input; -exponentially modulate the response y(t) to generate the signals ŷ n (t) = y(t)e −in p t . The value of n is in the range where n h is the maximum number of Fourier coefficients, equal to the number of harmonics, used to approximate the periodicity of the system; -compute the autospectrum of ŷ(t) using a standard approach (e.g., Welch's method); -apply classical identification techniques to extract the system modal parameters. In this work, a covariancedriven stochastic subspace identification (cov-SSI) algorithm has been applied.

Model description
The entire model can be divided into two main subcomponents: the wing-pylon assembly and the rotor. A render of the ATTILA multibody model is showed in Fig. 2.

Wing-pylon model
The wing-pylon model has been developed to reproduce the fundamental frequencies and mode shapes of an equivalent finite element stick model tuned to match the full-scale aircraft dynamics at the rotor hub.

FLIGHTLAB
The FLIGHTLAB model contains 16 elastic beam segments for the wing and 4 beam segments for the nacelle. The wing airloads are modeled using an enhanced lifting line model with a Peters-He [37] finite-state dynamic inflow model. Each structural beam segment is connected to a quasisteady airloads component, which uses 2D (AoA, Mach/ Reynolds) look-up table data to calculate the airloads.
The connection between the wing and the nacelle is modeled by three torsional spring-damper components, collocated and connected in series.
The shaft tilting hinge is modeled by a gimbal hinge, with a pitch-yaw stiffness that can be changed at run-time to switch between the downstop ON and OFF configurations (see Figs. 3 and 4).

MBDyn
The wing model consists of 3 finite volume three-node beam elements [3,16] for the stiffness part, and one body element for each node to model the inertial component. The nacelle part is divided into a tilting and a non-tilting part, both modeled as rigid bodies.
The parts are connected with deformable joints, which represent the flexibility of downstop and wing-pylon attachment. The aerodynamic loads are introduced through MBDyn's aerodynamic beam elements, based on simple strip theory, each linked to the corresponding structural element.

Rotor model
The ATTILA proprotor is a three-bladed stiff-in-plane rotor with a gimballed hub. It consists of the control chain, three blades, and the yoke.

FLIGHTLAB
In FLIGHTLAB, the blade is modeled using elastic beam components. The beam axis of each finite element is defined by the locus of shear centers (the elastic axis) of the physical blade. Appropriate sweep and droop rotations are applied to approximate the position of the elastic axis relative to the feathering axis.
The aerodynamics of the blade are modeled using look-up tables with a correction for unsteady circulatory effects. The transition between the airfoils is non-smooth. Provisionally, the induced velocity is modeled as a uniform inflow with the Glauert distribution.
The blade is connected to the yoke at two locations, at the inner and outer bearings. This creates a dual load path. The FLIGHTLAB solver is single load path-based, where calculations are performed sequentially from the blade tip towards the rotor hub. To facilitate the modeling of the dual load path resulting from the combination of blade and yoke, the root end of the blade is modeled as a separate beam, inverted and connected at the physical root via two-parent springs representing the inner bearing. In this setup, the root of the torque tube acts as a second "blade tip" as far as the solver is concerned.
The main "blade" load path consists of the yoke segments and the blade segments outboard from the outer bearing. The torque tube load path consists of the blade segments between the inner and outer bearings. The outer bearing is modeled as a series of three hinges with zero spring stiffness, allowing rotation around all three axes but constraining translation. Due to the single load path nature of the FLIGHTLAB solver, the inner bearing cannot be modeled as a hinge-slide combination. Instead, two perpendicular rigid flap/lag offsets are located at the inner bearing, perpendicular to the yoke. The inboard end of the torque tube is then connected to the rigid offsets through two two-parent translational linear spring-dampers that only constrain translation. The springs have been assigned high stiffness to approximate a rigid constraint. Both the springs and the rigid offsets are identical in length. The length of the spring is arbitrary but large enough to minimize the spring restraint in the axial direction of the yoke in the presence of relative displacements/rotations.
As mentioned, there are two offset/spring combinations: one in the flapwise and one in the edgewise direction. In this way, translation of the inboard end of the torque tube is constrained to the yoke at the inner bearing location, in both the flapwise and edgewise directions. Due to the length of  the springs and rigid offsets, the angle between them will be small, resulting in very small off-axis forces. As such, translation in the spanwise direction is nearly unobstructed, whereas rotation is free in all axes (see Fig. 5).
The control chain is modeled as a conventional swashplate arrangement. The rotor hub is rigidly connected to the tip of the pylon shaft when analyzing the half-wing model and rigidly connected to the inertial system when analyzing the isolated rotor. A bearing component drives the rotor rotation over the gimbal, which is connected to the yoke via an underslung offset.
The swashplate is located above the hub, connected via a rigid offset. The non-rotating swashplate node is connected to the pylon through a translational linear spring-damper, representing the collective spring stiffness of the control chain. A collocated gimbal hinge with spring-damper effects models the cyclic pitch spring stiffness.
The non-rotating swashplate node translates along the shaft axis, following the collective input through a controlled slider. Azimuthal rotation of the rotating swashplate node is achieved by a controlled hinge, which is slaved to the rotational speed of the hub. The swashplate mass is fixed to the rotating swashplate node and is required to avoid a singular matrix during linearization. Each pitch link is connected to the rotating swashplate and offset from the shaft through azimuthal rotation and rigid translation. Cyclic control inputs are introduced at the root end of the pitch links through a controlled slider. A two-parent linear spring-damper represents the pitch link stiffness and is connected to the pitch horn on the blade. The pitch link spring does not constrain rotation.

MBDyn
The blade and the yoke are modeled in MBDyn using the three-node beam element, similar to the modeling of the wing (Fig. 6).
To describe the orthotropy of blade and yoke, the stiffness matrix of the beam sections incorporates the offsets and relative rotations between the feathering axis and the neutral and elastic axes.
The aerodynamic model is constructed using MBDyn's aerodynamic beam. Each aerodynamic panel can incorporate aerodynamic twist variation and airfoil transitions. Along the blade, six airfoils are placed with a non-smooth transition. The rotor aerodynamics are completed by a uniform inflow model, based on momentum theory with Glauert's distribution and empirical corrections for tip loss and other standard airflow features. The resulting model represents a low-fidelity description of the rotor aerodynamics; however, in the tiltrotor aeromechanics literature, it is considered adequate for Whirl-Flutter. Future investigation will also consider mid-fidelity aerodynamics, involving the Vortex-Particle Method (VPM) presented in [47].
The blade is connected to the yoke in two locations, at the inner and outer bearing. In MBDyn, these two bearings are modeled with ideal rigid constraints: for the inner one, both flapwise and chordwise displacements are constrained, whereas for the outer, all three components of translation are constrained.
The control chain has a traditional helicopter-like configuration: it is formed by seven MBDyn nodes joined following the scheme of Fig. 7: -Pylon: this node physically connects the pylon extremity and the rotor; when the isolated rotor is analyzed, this node is clamped. -Airframe: this node receives the commanded collective and cyclic pitch controls. To decouple the two cyclic inputs, the node is positioned on a reference system that is rotated by the angle sp = tan −1 (x sp ∕y sp ) , where x sp and y sp are the locations of the pitch link attachment to the swashplate. -Fixed Swashplate: this node's in-plane displacement components and axial rotation are rigidly constrained to the airframe. To take into account the flexibility of the control chain, a collective spring, and two cyclic springs are positioned in between the airframe node and the fixed swashplate. -Rotating Swashplate: this node is connected to the fixed swashplate by a revolute hinge, and its axial rotation is constrained to the mast. -Engine: this node is connected to the mast by a torsional spring to reproduce the drive-train dynamics. -Mast: this node transmits the rotation to the hub and the rotating swashplate. It is connected to the pylon node by a revolute hinge. -Hub: This node is constrained to the mast node by a spherical hinge and a gimbal rotation constraint: these two combined elements create an ideal constant velocity joint.

Wing-pylon
To validate the wing-pylon model, a modal analysis has been performed comparing the frequencies and mode shapes of both the MBDyn and FLIGHTLAB models with a target NASTRAN stick model. In the NASTRAN model, the wing is modeled with 12 CBAR elements (two-node beam elements), whereas the nacelle is modeled with 4 CBAR elements. Lumped stiffness elements are used to join the two subcomponents. Compared with the NASTRAN model, the MBDyn model of the wing is discretized with a smaller number of beam elements (3, with 6 nodes vs. 12); instead, the FLIGHTLAB model is a one-to-one translation of the NASTRAN model.
At this stage, the comparison of the mode shapes is based on the Modal Assurance Criterion (MAC) defined by Eq. (26). In this equation, FEM i represents the ith mode shape calculated using the original finite element stick model, whereas the term MBD j represents the jth mode shape obtained from the two multibody codes.
The MAC matrix is presented in Fig. 8. The closer the matrix is to identity, the closer the results are. The differences between the results obtained from MBDyn and FLIGHT-LAB and those from NASTRAN appear to be negligible. Close inspection of the 6-DOF mode shape at the location of the rotor hub confirms the rather satisfactory correlation.

Rotor
To validate the dynamic behavior of the rotor, its frequencies in vacuo (i.e., neglecting aerodynamic loads) at different collective pitch settings have been evaluated and compared to the fan plots obtained with an equivalent full-scale CAM-RAD II rotor model. Moreover, the rotating blade mode shapes have been compared for different collective angles and with and without the presence of the drive train system. Remarkably good agreement has been obtained between the three models, despite their relative complexity.

Trim procedure
The trim targets depend on the configuration being investigated but are identical for both FLIGHTLAB and MBDyn.
Power-on trim is based on prescribing the desired rpm and finding the collective and cyclic pitch controls that achieve a target thrust and zero cyclic gimbal flap angle, increasing the airstream speed until the maximum torque is reached. From that point on, while the airstream speed is increased further, a constant torque trim is maintained to represent a steady powered descent.
During power-off trim, the rotor is de-clutched from the engine (i.e., no torque is transmitted), and the rotor speed is held constant by acting on the collective pitch while the airstream speed is increased.
In FLIGHTLAB, trim is achieved by a gradient-based iteration process. In power-off trim, the requested rotor speed is frozen at the target value while the swashplate controls are adjusted to achieve an average rotor azimuthal acceleration equal to zero. Post-trim, the rotor speed degree of freedom is released and the collective is fine-tuned to achieve the desired steady power-off rotor speed.
In MBDyn, power-on trim is achieved by starting from an initial guess of the collective angle and then applying the desired torque at the engine side. The desired rotor speed is maintained using a PI controller that sets the required swashplate collective as a function of the errors on the rotor angular velocity and its integral. The power-off condition is achieved by simply setting the applied torque to zero.

Whirl-Flutter prediction
Whirl-Flutter results mainly concern the low-frequency modes that characterize the aeroelasticity of the model in the configuration that will be tested in the wind tunnel, that is the wing flatwise and chordwise bending modes (respectively, called "wing bending" and "wing chord' in the following), the "wing torsion" mode, and the "pylon-yaw" mode. The latter mode is characterized by substantial participation of a rotation of the pylon-nacelle body relative to the wing about an axis loosely normal to the wing surface. This latter mode is particularly sensitive to the locking of the downstop mechanism. The focus is on evaluating the damping associated with those modes as the wind tunnel speed is increased and its sensitivity to several parameters of the problem, with particular attention on the prediction capabilities and accuracy of the algorithms considered in this work. Figure 10 shows the Whirl-Flutter stability predictions for the power-on configuration with the downstop engaged. The figure compares the modes obtained in FLIGHTLAB through a direct linearization and those obtained by applying the Periodic Operational Modal Analysis to MBDyn's results.

POMA mode shapes
The stability analysis of a complex non-linear system modeled using the multibody approach cannot be easily addressed using analytical methods based on a linearized expression of the governing equations. This problem has been previously addressed, in analogy to what is done here, extracting the characteristics of the system from numerical experiments. In [31], the stability properties of a complete tiltrotor multibody model were computed using the previously mentioned POD-based approach [43].
The method used in this work, however, has the benefit of extracting a prescribed number of periodic mode shapes starting from the main harmonic contribution in the collected data. The results will be presented independently for each identified wing mode. For each of them, the visualized mode shape and the participation factors computed performing the analysis using strain measurements are compared. This additional check allowed to verify the agreement of the results of the two analyses. Furthermore, after having performed this comparison, it was easier to identify the modes looking at the patterns of the extracted participation factors j n that are calculated as: where j n is the j − th eigenvector associated to the n − th harmonics.
As an example, the wing bending mode is considered. The channels used in the identification process are the wing root out of plane bending moment and additional (27)  ones related to the rotor, the gimbal angle, and the blades out of plane bending moments, to also capture the rotor contribution in the superharmonics. Figure 14 shows the complex mode indicator function for the above listed channels for the powered condition at minimum speed. The peaks associated with the wing bending mode and its relative superharmonics are clearly visible and distinguishable, assuring a high-quality identification. The wing bending frequency, b , is smaller than the rotor angular velocity, so the b − Ω superharmonic is reflected in the positive half-plane of the complex mode indicator function, resulting in the lower superharmonic in the Figure. Of particular interest is the response of the system at the identified superharmonics frequencies. Clearly, in the non-rotating frame, the wing subsystem is almost timeinvariant, while the measurements in the rotating frame, related to the rotor, exhibit a more periodic behavior dominating the two superharmonics. This is observed in the mode shapes identified using displacement measurements and by computing the participation factors of the identified eigenvectors using strain data.
As described in Fig. 15, it is evident that the wing contribution in the response is limited to the main harmonic; the participation of the wing bending moment to the eigenvectors is negligible in both superharmonics. Nevertheless, the rotor contribution is predominant and is mainly associated with a variation of the rotor gimbal angle. The bending of the wing and the flapping of the rotor dominate the overall response of this mode. The latter is induced by the transverse wind component the rotor experiences during wing out-ofplane motion.
Observing the mode shapes, depicted in Fig. 16, the behavior noticed from the computed participation factor is more easily visualized. In addition, the +Ω superharmonic shows, not only a major contribution of gimbal mode, but also of in-plane blade bending.

Discussion
In Fig. 10, the results from the two codes show similar trends for the wing bending mode (red) and the wing chord mode (blue). The FLIGHTLAB model predicts a lower damping value when considering the torsional mode (black) and pylon-yaw mode (cyan).
The firsts three modes (wing bending, wing chord, and wing torsion) show a good match in Fig. 11. A nonnegligible difference is found for the pylon yaw mode. The damping curve predicted by the MPE method changes slope at a lower speed compared to the one obtained with the POMA method. This discrepancy may be due to the different types of excitation used to trigger the time responses analyzed by the two methods. Since the MBDyn model is nonlinear, what is interpreted as a modal response may differ when excited by turbulence rather than a deterministic excitation applied to the blade pitch through the swashplate.
The comparison of the Whirl-Flutter stability predictions for the power-off configuration presented in Fig. 12 and obtained with MBDyn, identified using the POMA method, and FLIGHTLAB exhibits an overall good agreement. In this case, the wing bending mode identified in FLIGHTLAB shows a slope slightly steeper than that identified from MBDyn's results, whereas the damping predicted for the pylon yaw mode is higher in the MBDyn analysis. The wing chord and torsion modes, instead, match almost perfectly.
The damping factors shown in Fig. 13, identified from MBDyn's results with the two proposed algorithms, show an almost perfect agreement for the wing beam bending and the wing chord modes. When identified through the MPE method, the wing torsion mode does not show a change in slope. As in the power-on configuration, the pylon-yaw mode shows different trends between the two identification methods.

Conclusions
Two multibody models of the Next-Generation Civil TiltRotor-Technology Demonstrator (NGCTR-TD) wind-tunnel testbed, independently developed using FLIGHTLAB and MBDyn, are presented. Both models are in good agreement with a reference CAMRAD II rotor model and NASTRAN wing-pylon stick model. After thorough validation of the aeromechanics predictions, aeroelastic characteristics of the experimental set-up have been assessed utilizing three methods to identify its aeroelastic modes. The first method is based on state linearization as adopted in FLIGHTLAB. The Matrix Pencil Method and Periodic Operational Modal Analysis have been applied to time histories simulated using MBDyn. The results show good agreement for the wing beam bending, chord and torsion modes, with residual discrepancies for pylon-yaw. They are attributed to residual differences in the underlying dynamics of the models rather than to the identification methodologies. The last method can easily extract the shapes of both wing-pylon and rotor modes. Their visualization grants the possibility of a deeper understanding of the system behavior and instability mechanisms. The analyses confirm the absence of Whirl-Flutter instability in the airstream speed range of interest. The satisfactory agreement between the different approaches considered in this work gives us sufficient confidence in view of the continuation of the research towards the experimental verification of the aeroelastic stability of the system. Future work will address more accurate predictions using mid-fidelity aerodynamic models, model updating after preparatory and ground vibration testing of the wind-tunnel model, and correlation with test results when available.
Funding Open access funding provided by Politecnico di Milano within the CRUI-CARE Agreement. This research received funding from the Clean Sky 2 -H2020 Framework Program, under grant agreement N. 863418 (the ATTILA project).