Coupling Mid-Fidelity Aerodynamics and Multibody Dynamics for the Aeroelastic Analysis of Rotary-Wing Vehicles

Abstract

A mid-fidelity aerodynamic solver based on the vortex particle method for wake modeling, DUST, is coupled through the partitioned multi-physics coupling library preCICE to a multibody dynamics code, MBDyn, to improve the accuracy of aeroelastic numerical analysis performed on rotary-wing vehicles. In this paper, the coupled tool is firstly validated by solving simple fixed-wing and rotary-wing problems from the open literature. The transient roll maneuver of a complete tiltrotor aircraft is then simulated, to show the capability of the coupled solver to analyze the aeroelasticity of complex rotorcraft configurations. Simulation results show the importance of the accurate representation of rotary wing aerodynamics provided by the vortex particle method for loads evaluation, aeroelastic stability assessment, and analysis of transient maneuvers of aircraft configurations characterized by complex interactional aerodynamics. The limited computational effort required by the mid-fidelity aerodynamic approach represents an effective trade-off in obtaining fast and accurate solutions that can be used for the preliminary design of novel rotary-wing vehicle configurations.

1. Introduction

The design of complex rotary-wing vehicles, such as tiltrotors, represents a challenge for engineers and scientists. Being able to perform complete aeroelastic simulations considering the coupling of rotary-wing aerodynamics with structural dynamics of the vehicle is essential for the development of novel aircraft configurations [1]. Sophisticated numerical tools have been developed to support the analysis of rather different operating flight conditions typical of VTOL vehicles with increasing level of detail. Indeed, a substantial effort was spent in past years to develop sophisticated structural dynamics codes (CSD) that were effectively used for rotorcraft applications (e.g., References [2,3,4]). In detail, structural dynamics of rotary-wing vehicles was typically investigated using the multibody approach [5,6], which takes into account the nonlinear dynamics of the interconnected rigid and flexible bodies representing the aircraft components during transients. The multibody approach was also used to investigate aeroelastic phenomena, especially in airplane mode flight where whirl-flutter instabilities may occur [7]. A particular effort in this research field was spent at Politecnico di Milano, where, starting in the 1990s, a free general-purpose multibody software called MBDyn (http://www.mbdyn.org/) was developed, with the aim of gaining autonomous modeling capabilities of generic problems related to the dynamics of complex aeroelastic systems, specifically rotorcraft and tiltrotors [8]. A tiltrotor flight mission is characterized by vertical take off and landing, typical of helicopters, and by a cruise flight condition typical of fixed-wing airplanes. Transition between the two flight conditions occurs through the so-called conversion maneuver, in which the proprotor-nacelle system is tilted. Thus, an accurate model for the evaluation of unsteady loads occurring during transients must be considered for a correct analysis of the complex dynamics characterizing of these vehicles. However, multibody solvers, such as MBDyn, are typically equipped with simplified aerodynamics models based on Blade Element/Momentum Theory (BE/MT). This type of approach does not consider aerodynamic interference between rotors and the actual geometry of lifting surfaces, possibly leading to a misrepresentation of the aerodynamic loads and loss of information related to periodic actions.

The accurate computation of unsteady aerodynamic loads for rotary-wing vehicles was the main goal of several research groups that developed, in recent decades, different high fidelity Computational Fluid Dynamics (CFD) codes employing the use of high-fidelity Navier–Stokes equations solvers. To name a few examples, ONERA, University of Glasgow, DLR and Airbus Helicopters Deutschland and Politecnico di Milano developed, respectively, elsA [9], HBM3 [10], FLOWer [11], and ROSITA [12], high-fidelity CFD codes based on the block-structured grid, finite volume, and Chimera approach for the simulation of rotating bodies. These tools were purposely developed in Europe for rotorcraft and tiltrotor application studies [13,14,15]. Similarly, considerable research effort was dedicated in the USA to the numerical study of rotorcraft, particularly to tiltrotor aerodynamics, as shown for instance by the works by Meakin [16], Potsdam and Strawn [17], and Wissink et al. [18], where hover configurations were simulated using different implementations of the Navier–Stokes equations, and by the recent works of Lim, Tran et al. [19,20,21,22] that investigated the aerodynamic interaction between rotors and wing of the XV-15 tiltrotor aircraft.

Coupling of CSD codes with high fidelity CFD solvers was successfully investigated and implemented in the last two decades years for aeroelastic simulation of rotorcraft applications [23,24,25,26,27,28]. The coupled CSD/CFD numerical approach was successfully validated against experimental results, e.g., for the flutter calculations of a vertical tail model [29] and for the analysis of rotor blade structural loads of a complete helicopter model tested in a transonic wind tunnel [30]. Nevertheless, despite continuous advances in the field of high performance computing, coupled simulations of CSD and time-accurate URANS simulations of complete rotorcraft configurations still require a robust computational effort, not suitable for the preliminary design stage of novel VTOL aircraft configurations as tiltrotors, which requires a great number of simulations to reproduce the different flight conditions that characterize their mission.

In recent years, considerable effort was dedicated by several research groups to the development of mid-fidelity aerodynamic solvers based on the use of the vortex particle method (VPM) [31,32] for wake modeling. This numerical methodology showed a quite accurate representation of the aerodynamic interactions among several bodies typical of complex rotorcraft configurations and limited computational time with respect to URANS CFD simulations. To cite few examples, Lu et al. [33] developed an optimization methodology for helicopter design based on an viscous VPM model combined with an unsteady panel hybrid method. Alvarez and Ning developed a VPM-based code [34] for the investigation of multi-rotor configurations. Tan et al. [35] used a vortex-based approach coupled with a viscous boundary model to study rotor-to-rotor interactional problems occurring during shipboard operations [35]. Recently, Politecnico di Milano developed a novel, flexible mid-fidelity computational tool, called DUST (www.dust-project.org) aimed at representing a fast and reliable asset for the simulation of the aerodynamics of complex rotorcraft configurations, such as the electrical Vertical Take Off and Landing (eVTOL) aircraft. DUST is an open source code, released under MIT license, integrating different aerodynamic models for solid bodies, as thick surface panels, thin vortex lattices elements, and lifting lines elements. Moreover, a VPM method was implemented for wake modeling providing a stable Lagrangian description of free-vorticity flow field, which is suitable for numerical simulations of configurations characterized by strong aerodynamic interactions. DUST was thoroughly validated by comparison with both high-fidelity CFD simulations results and experimental data over complex rotorcraft configurations, such as eVTOLs [36] and tiltrotors [37]. Consequently, this novel open source tool is reaching maturity for the simulation of the aerodynamics of complex rotorcraft configurations accurately reproducing the interaction between rotors and wings.

The combination of a multibody dynamics solver with a mid-fidelity aerodynamic tool aims at representing an ideal trade-off between speed of execution and accuracy of the solution, aimed at the preliminary design of novel rotary-wing aircraft configurations. A novel open access aeroelastic tool was built by coupling MBDyn with DUST. The coupling of the two codes relies on the partitioned multi-physics coupling library preCICE [38], a very useful and robust tool for managing the communication between different solvers. After a brief description of the multibody dynamics and mid-fidelity aerodynamics solvers, this paper describes the details of the methodology used for the coupling. An interesting novelty proposed by this tool is the capability of modeling the deflection of a control surface, representing an essential aspect in the simulation of aircraft maneuvers. The coupled code is validated considering two test cases [39,40]. The first validation is provided by analyzing Goland’s wing, a numerical test case that was widely used in literature as a benchmark for flutter predictions. The second test case involves a rotary-wing application, using experimental data available for the XV-15 tiltrotor in hover. Then, the paper shows the capability of the tool to simulate a more complex application, namely the coupled simulation of a complete tiltrotor during a transient maneuver, a novel application for a coupled CSD/mid-fidelity aerodynamic tool in the rotorcraft field. Indeed, the control of these vehicles is often obtained by concurrently operating several control surfaces and actuation systems associated with the rotors, mixed by a complex Flight Control System (FCS) that takes into account the different flight conditions [41,42]. Thus, the design of the control surfaces and the selection of the actuators requires the ability to evaluate the unsteady aeroelastic loads acting on the aircraft during the maneuvers, to improve the vehicle response, increase effectiveness and efficiency, and reduce weight and complexity of the control system. Consequently, the coupled simulation of the complete tiltrotor model, representative of the Bell XV-15, was aimed at indicating the role of an accurate reproduction of rotary-wing aerodynamics and aeroelasticity, resulting from VPM and including the effects of the mutual interactions between the components of the complete aircraft, to correctly reproduce the dynamics of a roll maneuver. The results represent a key-feature for sizing the wing movable surfaces during the preliminary design of a novel vehicle. In fact, this work was performed in the framework of the FORMOSA and ATTILA projects, the former aimed at the design of the novel wing movable surface system of the NextGen Civil Tiltrotor aircraft (NGCTR), and the latter aimed at designing and testing a wind-tunnel model for the experimental verification of the whirl-flutter stability boundaries of the NGCTR. Both projects are developed within the framework of the Clean Sky 2-H2020 Program.

2. Description of the Coupled Software

2.1. Multibody Software MBDyn

MBDyn solves constrained dynamics of rigid and flexible systems by solving the Newton-Euler equations of motion of rigid and flexible bodies connected by kinematic constraints [8]. The problem is cast in first-order form as:

$$\begin{array}{cc}\mathbf{M}(\mathbf{x},\mathbf{t})\dot{\mathbf{x}}& =\mathbf{p}\end{array}$$

(1)

$$\begin{array}{cc}\dot{\mathbf{p}}& ={\mathit{\varphi }}_{/\mathbf{x}}^T\lambda +{\mathbf{f}}_i(\dot{\mathbf{x}},\mathbf{x},t)+{\mathbf{f}}_e(\dot{\mathbf{x}},\mathbf{x},t)\end{array}$$

(2)

$$\begin{array}{cc}\mathit{\varphi }\left(\mathbf{x}\right)& =\mathbf{0},\end{array}$$

(3)

where $\mathbf{x}$ is the vector of the kinematic unknowns, $\mathbf{p}$ is the vector of the momentum unknowns, $\lambda$ is the vector of the Lagrangian multipliers, $\mathbf{M}$ is a configuration (and possibly time) dependent inertia matrix, ${\mathbf{f}}_i,{\mathbf{f}}_e$ are arbitrary internal and external forces, $\mathit{\varphi }\left(\mathbf{x}\right)$ are the algebraic constraint equations (only holonomic constraints are shown for conciseness), and ${\mathit{\varphi }}_{/\mathbf{x}}^T$ is the Jacobian matrix of the holonomic constraints with respect to the kinematic unknowns. Each structural node instantiates six force and moment balance equations (2), while only nodes with associated inertia properties instantiate the equations that define the momenta (1). The details associated with the handling of finite rotations are not made explicit; the interested reader is referred to Reference [8]. Additional states, associated with scalar fields (for example hydraulic pressure, temperature, voltage) and, thus, the associated differential and/or algebraic balance equations, can be taken into account through specialized sets of nodes.

The nodes that describe the kinematics of the structural problem can be connected either by elastic/viscoelastic internal forces (namely lumped structural components [43], beams [44,45], shells [46], Component Mode Synthesis (CMS) elements [47]) expressed by ${\mathbf{f}}_i$, with a variety of viscoelastic constitutive laws, or by the kinematic constraints of Equation (3). Simple aerodynamics can be modeled by built-in elements that exploit the 2D strip theory model by look-up tables of the aerodynamic coefficients and classical rotor inflow models based on the momentum theory.

Coupling with an external solver is performed through external force element. This type of element communicates with the external solver by passing the kinematics of the model and receiving the corresponding loads. The coupling can be loose (i.e., explicit scheme) or tight (i.e., implicit scheme). For further details on the possible kind of coupling with an external solver, the reader is referred to Reference [48].

The external forces can be formulated directly in the absolute frame, or referred to a reference node. In the former case, operations are straightforward; in the latter, the kinematics are first reported in the reference frame of the reference node and then sent to the peer along with the motion of the reference node. The latter returns nodal forces and moments oriented according to the reference frame of the reference node.

2.2. Mid-Fidelity Aerodynamic Software DUST

DUST is an open source software designed by using object-oriented paradigms of the latest FORTRAN standards. DUST is a mid-fidelity aerodynamic software released under open source MIT license. The code was implemented using object-oriented paradigms of FORTRAN. The aerodynamic problem, based on incompressible flow hypothesis, is formulated by using the Helmholtz’s decomposition of the velocity field that is given by the irrotational (${\mathbf{u}}_{\phi }$) and solenoidal ${\mathbf{u}}_{\psi }$ contributions of velocity as $\mathbf{u}={\mathbf{u}}_{\phi }+{\mathbf{u}}_{\psi }$. Different classical potential-based elements are implemented in the code, such as lifting line elements [49,50], typically used to model aerodynamics of slender lifting bodies as blades, surface panels [51], and vortex lattice elements (Reference [52], Section 10.4.3), typically used for aerodynamic modeling, respectively, of thick and thin solid bodies. Even if incompressible potential flow is considered in the code assumptions, compressibility effects are considered in the computations. In particular, a Prandtl-Glauert correction [53] is considered for steady aerodynamic loads calculation with surface panels and vortex lattice elements. On the other hand, lifting line elements intrinsically include compressibility and viscous effects introduced by Mach-depending tabulated sectional aerodynamic coefficients, obtained by experiments or by 2D RANS numerical simulations.

DUST simulations only require surface meshes, as the core of the mid-fidelity aerodynamic code is the vortex particles method (VPM) [31,32], a Lagrangian grid-free approach that is used to model the free vorticity of wakes and that does not require a volume mesh of the flow surrounding the object of investigation. Vortex particles reproduction of the flow enabled to obtain a robust representation of interacting wakes issued by lifting surfaces and bodies, as typically occurs in rotary-wing vehicles applications. In particular, in order to reduce the computational cost related to reproduce vortex particle interactions, a Cartesian Fast Multipole Method (FMM) [54], providing a background hierarchical decomposition of the domain into clusters of cells is implemented.

The code calculates aerodynamic loads using a different method according to the element considered to model the problem. For lifting line elements, the aerodynamic loads are calculated starting from the tabulated sectional aerodynamic coefficients used to solve their intensity. For vortex lattice elements, aerodynamic loads are calculated using the unsteady formulation of the Kutta-Joukowsky theorem. For surface panel elements, aerodynamic loads are calculated using the unsteady formulation of the Bernoulli theorem, considering the vorticity of the flow.

For a detailed insight on the mathematical formulation of the DUST code, including the implemented governing equations, the reader is referred to a recent work published by some of the present article authors [36].

2.3. Description of MBDyn-DUST Coupling

Communications between DUST and MBDyn are managed by preCICE (Precise Code Interaction Coupling Environment), a coupling library for partitioned multi-physics simulations, originally developed for fluid-structure interaction and conjugate heat transfer simulations. preCICE offers methods for transient equations coupling, communication means, and data mapping schemes. It is written in C++ and offers additional bindings for C, Fortran, MATLAB, and Python. preCICE (https://github.com/precice/) is a free software released under the LGPL3 license. While MBDyn uses its own C++, C, Fortran, and Python Application Programming Interface (API) to communicate with external software without any further modification to the C++ source code, no API was already available in DUST. Thus, new Fortran modules collecting all the classes, subroutines, and functions required by the adapter for preCICE were implemented to support coupling with DUST. The optional coupling with external codes was managed through pre-processor directives.

A new adapter was implemented for supporting the communication of all the kinematic variables (position, orientation, velocity, and angular velocity) plus forces and moments acting on the nodes of a MBDyn model exposed through an external structural force element. Figure 1 shows a scheme of the communication and information exchange, managed through the adapters for the two solvers.

The interface between structural and aerodynamic grids is obtained as a weighted average of the distance between the nodes of the two grids and it is used for motion interpolation and consistent force and moment reduction. Figure 2 shows the q nodes of the structural grid, namely $Q_q$. The centers and the vertices of each aerodynamic mesh are, respectively, $P_e$ and $P_p$, where e and p are the corresponding indices.

The kinematic variables, ${\varphi }_p$, of a point p positioned on the aerodynamic surface of a DUST component is evaluated as the following weighted-average,

$${\varphi }_p=\sum _q{w}_{pq}{\varphi }_q,$$

(4)

where ${\varphi }_q$ is the same kinematic variable associated with the qth structural node of the MBDyn model.

Weights $w_{pq}$ could be any set of non-negative real numbers, satisfying the normalization conditions

$$\sum _q{w}_{pq}=1\forall p,$$

(5)

since they define the weighted average of the variables associated with the structural nodes q on the aerodynamic nodes p. These coefficients could be proportional to some negative power, defined as a user input, of a norm of the vectors $(P_p-Q_q)$. As an example, using the local coordinates in the reference configuration ${\mathit{r}}_{pq}$, the norm of these vectors can be defined as:

$$\parallel (P_p-Q_q){\parallel }^2:={\mathit{r}}_{pq}^T\mathbf{W}{\mathit{r}}_{pq},$$

(6)

where $\mathbf{W}$ is a positive (semi-)definite matrix, providing an anisotropy degree of freedom to the user in defining the (semi-)norm [55]. Threshold values and maximum number of influencing weights are two criteria, defined as user inputs, used to restrict the average only to the significant structural nodes for each aerodynamic point.

2.3.1. Kinematic and Load Variables

The position of a point P in the global reference frame g of the aerodynamic surface is evaluated as

$${(P_p-O)}^g=\sum _{Q\in I_P}{w}_{pq}\left\{{(Q_q-O)}^g+{\mathbf{R}}_Q^{r\to g}(P_p-Q_q)\right\},$$

(7)

where $Q\in I_p$ indicates the subset of structural points $Q_q$ that belong to the $I_p$ aerodynamic component, ${(Q_q-O)}^g$ is the distance from the origin of the $Q_q$ structural point, and ${\mathbf{R}}_Q^{r\to g}(P_p-Q_q)$ rotates in the global coordinates the distance between the aerodynamic point and the structural one. Its angular velocity ${\mathbf{\omega }}_P$ and velocity ${\mathit{v}}_P$, respectively, are

$$\begin{array}{cccc}\hfill {\mathbf{\omega }}_P& =\sum _{Q\in I_P}{w}_{pq}{\mathbf{\omega }}_Q,\hfill & \hfill {\mathit{v}}_P& =\sum _{Q\in I_P}{w}_{pq}\left\{{\mathit{v}}_Q+{\mathbf{\omega }}_Q\times (P-Q)\right\}.\hfill \end{array}$$

(8)

The aerodynamic forces and moments are evaluated at the evaluation points $P_e$ located at the centers of each panel and then transferred to the structural nodes using the summation of forces and transport of moments as follows:

$$\begin{array}{cccc}\hfill {\mathit{f}}_Q& =\sum _{e\in J_Q}{w}_{qe}{\mathit{f}}_e,\hfill & \hfill {\mathit{m}}_Q& =\sum _{e\in J_Q}{w}_{qe}\left\{{\mathit{m}}_e+(P_e-Q_q)\times {\mathit{f}}_e\right\},\hfill \end{array}$$

(9)

where $e\in J_Q$ indicates the subset of evaluation points that belong to each sub-component $J_Q$, and the weights $w_{qe}$ are calculated using Equations (6) and (7), by computing the distance between each structural node and evaluation point.

2.3.2. Implementation of Software Coupling

Figure 3 shows the flow of information during the coupled simulation between the two solvers for an implicit tight serial scheme. First, the object precice of class t_precice is declared for handling a coupled simulation through the preCICE library. This object is used both for managing data communication, and for updating coupled components of the aerodynamic model. Then, the instance of DUST participating in the coupled simulation is created, reading the XML preCICE configuration file. After some preliminary operations, the mesh used to couple the codes is defined, and the fields involved in the communication are initialized. Initialization of the Fast Multiple kernels, wake, and linear system follows. Before the time loop starts, communication is first established between the coupled codes. The time loop starts with the update of DUST’s explicit aerodynamic elements (lifting lines and actuator disks). A checkpoint of the exchanged fields is stored, to be reloaded during sub-iterations of preCICE’s implicit coupling. DUST then receives the kinematic variables of the structural nodes from the external software, MBDyn, and updates the surfaces of the coupled components and the near-field wake elements. The linear system is then updated and solved, calculating the strengths of the vortexes of the surface panels and vortex lattice elements.

The solution of the non-linear lifting line problem follows. The circulation $\mathbf{\Gamma }$ is computed by using analytical expression from the Kutta-Joukowski theorem, where lift coefficient is evaluated from tabulated sectional data. Once the intensity of the surface singularities has been evaluated, surface pressure distribution and elementary forces and moments are retrieved using the unsteady Kutta-Joukowski theorem for the vortex lattice and lifting lines elements, as well as the unsteady Bernoulli theorem for the 3D-panels.

Aerodynamic forces and moments are reduced to the nodes of the interface between the aerodynamic and structural meshes and sent to MBDyn. A convergence check on the kinematics variables follows. If convergence is not reached, the checkpoint fields are reloaded, and a new sub-iteration begins. If convergence is attained, the time step is finalized, saving the status and updating the wake and the geometry of the uncoupled components for the next time step.

2.4. Hinged Surfaces Modeling

Modeling the deflection of a control surface is essential for the simulation of aircraft maneuvers. In MBDyn, the deflection of a control surface can be easily modeled as a rigid/deformable body properly constrained to the fixed part of the vehicle to allow only its relative rotation about the hinge axis. In DUST, the possibility to include a control surface in the aerodynamic mesh has been only recently introduced. In the following, the description of the implemented model for hinged surfaces in DUST is introduced with a two-dimensional example first, and then extended to three-dimensional deformable components.

As outlined by the scheme in Figure 4 (left), in a two-dimensional problem the control surface can be defined in the local reference frame of the component, by means of the hinge axis position H, the chord-wise direction $\xi$, and a blending region $[-u,u]$ introduced to avoid irregularities in the mesh as the surface rotates by an angle $\theta$. Moreover, in the 2D modeling, the rotation axis $\widehat{\mathit{h}}$ is assumed to be orthogonal to the plane of the airfoil.

As shown in Figure 4, a orthonormal reference frame for the hinge is defined with origin in H and axes $\widehat{\mathbf{\xi }}$, $\widehat{\mathbf{\eta }}=\widehat{\mathit{h}}\times \widehat{\mathbf{\xi }}$. The position of a point with respect to this reference frame is

$$\mathit{r}=\xi \widehat{\mathbf{\xi }}+\eta \widehat{\mathbf{\eta }}+h\widehat{\mathit{h}}.$$

(10)

Three regions are defined using the coordinates based on this reference frame:

• $\xi \le -u$: no influence of the control surface rotation;
• $\xi \ge u$: rigid rotation about the hinge:

  $$\Delta \mathit{r}=sin\theta \widehat{\mathit{h}}\times \mathit{r}+(1-cos\theta )\widehat{\mathit{h}}\times \widehat{\mathit{h}}\times \mathit{r};$$

  (11)
• $-u\le \xi \le u$: blending region to avoid irregularities, defined as an arc of a circle whose center is located at point C and whose radius is:

  $$\overline{C{C}^{\prime }}={\displaystyle \frac{\overline{C^{\prime }H}}{tan{\displaystyle \frac{\theta }{2}}}}\mathrm{where}\overline{C^{\prime }H}=u.$$

  (12)

In a three-dimensional problem, the reference configuration of a control surface for a generic swept wing is defined in the wind axis reference frame of the component, as shown in Figure 5.

The aerodynamic sections that are involved in the deflection of the control surface are, thus, the ones that satisfy the condition $y\left(A\right)<y\left(P\right)<y\left(B\right)$, where $y\left(P\right)$ is the ordinate of the $P_i-th$ aerodynamic mesh point expressed in the wind reference system. As in the 2D case, one can define the three regions for each stripe identified at the previous point. The y coordinate of the origin of the sectional reference frame is determined by linear interpolation between points A and B.

When the movable surface is coupled with a structural component in MBDyn, the orientation of the hinge in DUST nodes comes from the orientation of the nodes of the MBDyn model. The rotation axis is defined as $\widehat{\mathbf{h}}={\displaystyle \frac{(B-A)}{\parallel B-A\parallel }}$. Each point of the movable surface is linked to the hinge nodes, obtaining its kinematic variables as a weighted average of the motion induced by the rotation of the hinge nodes. Weights $w_{ph}$ are evaluated using only the h components of the vectors connecting the control surface points to the hinge nodes.

The weights $w_{ph}$ are then combined to the interpolations weight w described in Section 2.3 to allow the deformation of the structure, as shown in Figure 6. The weights are calculated for each point of the aerodynamic mesh, whereas the weights $w_{ph}$ are specifically used to impose the rigid deflection of the movable surface.

3. Validation of the Coupled MBDyn-DUST Software

3.1. Goland’s Wing

This section presents the application of the coupled software to Goland’s wing [56], which is widely used in the literature to test and validate aeroelastic codes. This low aspect ratio wedged wing ($AR\approx 3.33$) is also interesting to highlight the impact on flutter calculations of 2D and 3D aerodynamic models. Figure 7 shows the layout of the problem. EA indicates the elastic axis; CG indicates the center of gravity axis. $U_{\infty }$ is the free-stream velocity. All the relevant geometrical and structural properties are reported in

Table 1. Goland’s wing properties and flight condition [57].

Wing properties
Half span	6.096 $\mathrm{m}$	Inertia/unit span	8.64 $\mathrm{k}\mathrm{g}\mathrm{m}$
Chord, c	1.8288 $\mathrm{m}$	Mass/unit span	35.71 $\mathrm{k}\mathrm{g}/\mathrm{m}$
Elastic axis	0.33$·$c	Torsional stiffness	0.99$\times {10}^6$ $\mathrm{N}{\mathrm{m}}^2$
Center of gravity	0.43$·$c	Bending stiffness	9.77$\times {10}^6$ $\mathrm{N}{\mathrm{m}}^2$
Flight condition
Air density	1.020 $\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$	AoA Perturbation	0.05 deg

. They have been obtained from Reference [57].

Two aerodynamic meshes were considered. In the first case, the wing was modeled as a flat plate using vortex lattice (VL) elements, while, in the second case, it was modeled with surface panels (SP) reproducing its geometrical shape and thickness. The results of a convergence analysis on the computed aerodynamic loads indicated the need to use 30 elements in the span-wise direction to obtain a correct spatial discretization for both models. The VL flat mesh requires 30 uniform divisions in the chord-wise direction, while the SP model requires 30 divisions for the lower and upper side of the wing, with half-cosine refined distribution at leading edge. The structural model built in MBDyn consists of four beams based on a $C^0$ beam discretization based on the finite volume concept proposed in Reference [44]. A satisfactory number of beam elements was obtained by satisfying a convergence requirement on the first four modes of the wing, as indicated by comparison of the frequencies of the first four normal vibration modes of the wing computed by MBDyn with respect to Goland’s work [56] and NASTRAN shown in

Table 2. Comparison of the first four natural frequencies computed for the Goland’s wing [56].

	Goland [56] Hz	NASTRAN Hz	MBDyn Hz
1st Bending	7.66	7.66	7.66
1st Torsion	15.24	15.24	15.21
2nd Bending	38.80	38.59	38.54
2nd Torsion	55.33	54.84	54.79

.

In the present study, tight coupled time-marching simulations were performed using a time discretization of 0.001 s and considering for the evolution of the wake particles a fourteen chords long box behind the wing, resulting in a developed wake made by about three thousand vortex particles. The computation time to perform a time-marched coupled simulation with a total duration of one second, using a workstation equipped with an Intel^® Core^TM i9-9980XE processor running on a base frequency of $3.00$ GHz, with 18 physical cores and 2 threads for each core, was about 19 min for VL and about 29 min for SP.

To study the flutter instability, a non-zero angle of attack of 0.05 deg was introduced as perturbation, as in Reference [58]. The frequency and damping of the response were identified from the time history of the wing-tip deflection using the matrix pencil estimation (MPE) method [59]. Figure 8 presents the z-displacement of the last structural external node across the flutter onset computed with the SP model. The red line corresponds to a stable damped response, whereas the blue line shows the condition of incipient flutter, in which a constant amplitude free oscillation is reached. The green line shows the unstable response at a speed greater than the flutter speed. A corresponding representation of the deformed mesh associated with the bending-torsion instability and the related distribution of the pressure coefficient is shown in Figure 8 (left). These results indicate the ability of the coupled simulation to correctly capture fixed-wing flutter.

Table 3. Comparison of flutter speed and frequency computed for Goland’s wing.

Author	Model	${\mathit{V}}_{\mathit{f}}$ ms${}^{-1}$	${\mathit{f}}_{\mathit{f}}$ Hz
Goland [56]	Analytical	137.2	11.25
Patil et al. [57]	Intrinsic beam + strip theory	135.6	11.17
Wang et al. [60]	ZAERO	174.3	-
Wang et al. [60]	Intrinsic beam + UVLM	163.8	-
Murua et al. [58]	SHARP, Displacement beam + UVLM	165.0	10.98
Present Work	MBDyn’s built-in strip theory	135.1	11.07
Present Work	DUST (VL)-MBDyn	168.2	10.84
Present Work	DUST (SP)-MBDyn	174.2	11.06

reports a comparison of flutter speed and frequency computed by several independent authors. Results computed with the coupled code are in quite good agreement with those obtained by similar codes using 3D aerodynamic models [57,58,60]. In detail, the discrepancy with the results obtained with the same MBDyn structural model, but using its built-in aerodynamics based on two-dimensional unsteady strip theory, indicates the superior capability of the coupled code in the investigation of aeroelastic problems. The comparison of the results from the literature reported in Table 3, obtained with 2D and 3D models, confirms the need of a three-dimensional aerodynamic model for a correct and realistic flutter analysis of low aspect ratio wings.

To evaluate the capabilities of the coupled code using the different employed aerodynamic models, Figure 9 presents the frequency and damping of the first beam torsional mode of the wing as functions of the free-stream velocity.

The numerical results of the coupled simulations using a panel mesh (SP) show a slightly higher aerodynamic damping with respect to those obtained using VL. A predicted flutter speed increase of approximately $3.7\%$ is observed (Figure 9 right). Nevertheless, the quite limited differences in the results of the two models indicate that, for simple configurations, without complex aerodynamic interaction between bodies, a vortex lattice mesh is more convenient than surface panels, as the computational cost reduces significantly.

3.2. XV-15 Proprotor in Hover

A rotary-wing test case was also considered to validate the MBDyn-DUST coupling. It consists of the XV-15 proprotor equipped with metal blades, a three-blade stiff-in-plane rotor with gimballed hub. For the work presented here, the gimbal universal joint has been modeled as ideal homokinetic joint, neglecting the 2/rev components caused by rotor flapping [61]. This test was selected thanks to the availability in the literature of the rotor geometrical and structural data [62,63] and of experimental data in hover condition. The experimental data were obtained during the test campaign that took place at the Outdoor Aerodynamic Research Facility (OARF), which is described in the work by Felker et al. [64].

The layout of a control chain representative of that of the XV-15 proprotor is shown in Figure 10. It was modeled using information from Reference [63]. The flowchart in Figure 11 details the blade pitch control system implemented in MBDyn. The role of each component is:

• Pylon: this node represents the actual connection between the pylon extremity and the rotor; when the isolated rotor is analyzed, this node is grounded.
• Control: through this node the controls (cyclics and collective) are prescribed. To decouple longitudinal and lateral cyclic inputs, the node is defined in a reference system that is rotated about the shaft by the angle ${\psi }_{sp}={tan}^{-1}(x_{sp}/y_{sp})$, where $x_{sp}$ and $y_{sp}$ are the in-plane components of the point where the pitch link is connected to the swashplate.
• Fixed Swashplate: this node represents the non-rotating portion of the swashplate. Its in-plane displacements and rotation about the shaft axis are constrained to the Control node. To account for the flexibility of the control chain, it is connected to the Control node by three equally radially spaced elastic rods.
• Rotating Swashplate: this node is connected to the Fixed Swashplate by a revolute hinge, and to the Mast node (defined later) by a joint that makes it rotate along with the shaft itself.
• Collective Head: this node is connected to the Control node by a deformable spring along the shaft axis, to account for the flexibility of the collective control path; it is positioned in the rotating reference system.
• Head Rocket Arm: this node is connected to the Collective Head node by a revolute hinge, and to the rotating swashplate by means of a rod (cyclic tube). Then, through the pitch link, the cyclic and collective commands are transmitted to the blade.
• Engine: this node represents the torsional inertia associated with the engine.
• Mast: this node drives the Hub (defined later) and the Rotating Swashplate nodes. It is connected to the Pylon node by means of a revolute hinge, whose angular velocity is prescribed.
• Hub: this node is constrained to the Mast node by a spherical hinge and a gimbal joint [65]; the combination of these two joints represents an ideal constant velocity joint.
• Yoke: this component represents the rigid body that connects the hub with the blade. The blade-to-yoke connection is obtained using a single load path through a revolute hinge.
• Blade: this node condenses both the inertial and aerodynamic properties of each proprotor blade.

All rotor data were taken from the original CAMRAD II model presented in Reference [62], here reported for clarity in

Table 4. XV-15 main rotor data [62].

Rotor Data
Blades	3
Solidity	0.0891
Radius	3.81	$\mathrm{m}$
Precone	2.5	deg
Chord	0.3556	$\mathrm{m}$
Nominal speed	589	RPM

.

Table 5. XV-15 Blade airfoil distribution [64].

Airfoil Data
Profile	start	end
NACA 64-935	0.09	0.13
NACA 64-528	0.13	0.34
NACA 64-118	0.34	0.655
NACA 64-(1.5)12	0.655	0.9
NACA 64-208	0.9	1

presents the blade span-wise airfoil distribution, from Reference [64].

For validation purposes, simulations using MBDyn alone, DUST alone, and the coupled solvers have been performed. In detail, in the MBDyn alone simulations the blade aerodynamics was modeled with the Blade Element/Momentum Theory using a Glauert model [66] for the induced velocity. A tip loss correction was introduced to scale the value of the normal aerodynamic force with a prescribed span-wise weight.

The aerodynamics of each blade was modeled in DUST using lifting line elements, naturally encompassing both compressibility and viscous effects, as described in detail in Reference [36]. This aerodynamic model provides accurate results on high aspect ratio bodies, as blades are, while being computationally very efficient, as shown in Reference [67].

The DUST aerodynamic model of the full-scale XV-15 proprotor was built considering the airfoil geometry, chord, twist, sweep, and dihedral distributions reported in Reference [68]. Lifting line elements were used for the aerodynamic model of the three blades. The tabulated aerodynamic performance coefficients of the blade airfoils were calculated using XFOIL [69] in the range of angle of attack before stall. The two-dimensional aerodynamic load curves are extended to the $[-180°:180°]$ range of angles of attack using the methods described in References [2,70]. While the MBDyn model applies a non-smooth transition between adjacent airfoil sections, DUST interpolates the aerodynamic properties of the airfoil sections used to build the blade model. A convergence analysis on computed rotor thrust and torque indicated the need to use 40 lifting line elements for the spatial discretization of each blade and a time discretization of 100 time steps to simulate a complete rotor revolution. For the evolution of the wake particles, a 5 rotor radius long box was considered in the simulations, resulting in a developed wake of around 52 thousand vortex particles. The computational time required to perform a DUST alone simulation over 5 complete rotor revolutions was about 324 s, while tightly coupled simulations using the same aerodynamic DUST model took around 395 s using the 18-core workstation described in Section 3.1.

Figure 12a shows the solidity-scaled thrust coefficient $C_T/\sigma$ versus torque coefficient $C_Q/\sigma$ curves computed for the XV-15 rotor in hover with MBDyn alone, DUST alone, and the coupled MBDyn-DUST solver, compared with experimental data reported in Reference [64]. It is worth noticing that:

(i)

the curves from the MBDyn-DUST coupled simulation agree fairly well with experimental data, particularly at high values of thrust coefficient, thus indicating the suitability of the coupled tool to improve accuracy in the evaluation of rotor performance with a quite limited computational effort;

(ii)

on the contrary, the curves from the MBDyn alone and DUST alone simulations both appreciably depart from experimental data and the coupled simulation results, although somehow capturing the observed trends;

(iii)

the curves from the MBDyn alone and DUST alone simulations appear to be in good agreement with each other, especially at high values of thrust coefficient;

(iv)

however, this agreement actually tells two different stories: MBDyn alone uses a poor aerodynamic model, but correctly describes the kinematics of the problem, whereas DUST alone uses a better aerodynamic model but with prescribed kinematics; it appears that, for this specific problem, a more accurate structural model and a more accurate aerodynamic model alone give a comparable contribution to model fidelity; clearly, this observation can hardly be generalized;

(v)

it also appears that, in the MBDyn alone analysis, the presence of an inflow model and of standard, well understood and accepted empirical corrections, can substantially improve the quality of the prediction, compared to the significant error that was observed in Goland’s wing flutter analysis.

4. Simulation of the Roll Maneuver on the Complete XV-15 Tiltrotor

A full-span aeroelastic model reproducing the Bell XV-15 research aircraft equipped with rigid metal blades [62] is considered for the coupled roll maneuver simulations. The purpose of this analysis is to show the ability of the code to model complex problems, particularly those involving the deflection of movable control surfaces, and to give an indication of the possible outcome of the proposed novel coupled numerical tool.

4.1. Numerical Model

The numerical model of the XV-15 tiltrotor is built considering the full scale dimensions and components of the aircraft. The model includes the fuselage, the horizontal and vertical tail-planes, the wing equipped with control surfaces (flap and flaperon) and the two proprotors with the corresponding nacelles. The main geometric characteristics of the entire vehicle, including the airfoils, are reported in

Table 6. Geometric characteristics of the XV-15 [41].

	Wing	Horizontal Tail	Vertical Tail
Airfoil	NACA 64A223	NACA 64015	NACA 0009
Span	9.8 m	3.91 m	2.34 m
Mean aerodynamic chord	1.60 m	1.20 m	1.13 m
Sweep ($c/4$)	$-{6.5}^{\circ }$	$0^{\circ }$	${31.6}^{\circ }$
Dihedral	${2.0}^{\circ }$	$0^{\circ }$	-
Incidence	${3.0}^{\circ }$	$0^{\circ }$	$0^{\circ }$
	Flap	Flaperon
Span along hinge line	1.30 m	2.40 m
Chord/Wing chord	$0.25$	$0.25$
Maximum deflection	${75}^{\circ }$	${47}^{\circ }$

.

The aerodynamic model uses different types of aerodynamic elements: lifting lines for the blades of the proprotor, using the same spatial discretization of Section 3.2, and surface panels for all the other aerodynamic components: 2514 elements are used for the fuselage, 1620 for the tail planes, 835 for each nacelle, and 3500 for the wing. The aerodynamic mesh of the complete XV-15 tiltrotor is shown in Figure 13.

The dynamic model includes:

• the wing, modeled as a rigid body, including flaps and flaperons;
• the fuselage and the empennages, modeled as rigid bodies, including the rudders and the elevator control surfaces;
• the pylon/nacelle system, mounted to the wing-tip; it can be tilted with respect to the wing to reproduce tiltrotor in airplane mode (APMODE), helicopter mode (HEMODE), or in any intermediate attitude;
• the rotor with the exact kinematics of the blade pitching mechanism, as described in Section 3.2.

The inertial structural data used for the multibody model were taken from the CAMRAD II and NASTRAN models presented in Reference [62].

Table 7. Center of mass, mass, and inertia tensor of the complete aircraft [62].

$x_{CG}$	$y_{CG}$	$z_{CG}$	M
7.214 $\mathrm{m}$	0 $\mathrm{m}$	−0.425 $\mathrm{m}$	5896.8 $\mathrm{k}\mathrm{g}$
$I_{xx}$	$I_{yy}$	$I_{zz}$	$I_{xy}$	$I_{xz}$	$I_{yz}$
67,512.6 $\mathrm{k}\mathrm{g}{\mathrm{m}}^2$	335,546 $\mathrm{k}\mathrm{g}{\mathrm{m}}^2$	398,235 $\mathrm{k}\mathrm{g}{\mathrm{m}}^2$	4.4 $\mathrm{k}\mathrm{g}{\mathrm{m}}^2$	23,136.5 $\mathrm{k}\mathrm{g}{\mathrm{m}}^2$	−0.5 $\mathrm{k}\mathrm{g}{\mathrm{m}}^2$

reports the inertia data of the complete aircraft considering a reference system having the x-axis pointing towards the tail, the y-axis aligned with the right wing, and the z-axis pointing upward according to the right-hand rule.

From the structural point of view, the flaperon is modeled as a rigid body connected to two nodes representing the hinges. The first node is constrained through a spherical hinge, while the second is constrained by an inline joint. The rotation of the flaperon is obtained by imposing a prescribed rotation on the second constraint. The combination of these three joints makes the constraint statically determined. The line connecting the two nodes identifies the axis of rotation of the movable surface; the aerodynamic mesh is deformed according to the methodology illustrated in Section 2.4.

To evaluate the aerodynamic effects of the tail-planes and their interaction with the wake of the rotors on roll performance, coupled simulations were performed for three different aircraft configurations, according to

Table 8. Aircraft configurations tested for roll maneuver coupled simulations.

Configuration	Airframe	Tail	Rotor
Configuration I (C.I)	✓	✗	✗
Configuration II (C.II)	✓	✓	✗
Configuration III (C.III)	✓	✓	✓

and Figure 14. First, the airframe only was considered; then, the tail-planes were added, and, finally, the rotors, increasing the complexity and the completeness of the model. The structural model and the relative mass properties were the same for the three simulated configurations.

A trimmed flight condition reported by Ferguson [71] was considered for the roll maneuver coupled simulations. The flight condition and trim parameters are reported in

Table 9. Flight condition parameters used for roll maneuver simulation [71].

Flight Condition
Air density	1.225	$\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$
Speed	72.022	$\mathrm{m}/\mathrm{s}$
Trim parameters
Mode	Airplane
Pitch angle	6.944	deg
Rotor speed	517	RPM
Rotor collective	29.5015	deg
Elevator	−1.2398	deg

.

During the roll maneuver simulations, the flaperons are deflected according to a step function from ${\delta }_a=0°$ to ${\delta }_a=\pm 20°$. The control surfaces start moving after 0.5 s from the beginning of the simulation, to make sure the maneuver starts after any initial aerodynamic transient vanish. At the same time, the roll degree of freedom of the entire model is unlocked. In the simulations, the aircraft rolls about the longitudinal axis, positive starboard (right) wing up. Yaw rotation is about the vertical body axis, positive nose left, while pitch rotation is about the axis normal to longitudinal plane of symmetry, positive nose up.

Considering a tightly coupled simulation using 145 time steps over a complete blade rotor revolution, the computational time to perform a one second long simulation for the full vehicle C.III was about 8 h using the previously mentioned 18-core workstation. In detail, the full developed wake for C.III consisted of approximately 84 thousand particles, considering a box domain twice the length of the aircraft in the x-direction. No convergence studies were performed for the full vehicle coupled simulation. Indeed, the choice of the selected spatial and time discretizations for the full vehicle was dictated by best practices (see Reference [37]) and by the will to limit the computational effort of the simulations. Indeed, one of the goals of the present activity is to show the capabilities of the novel coupled numerical tool for the preliminary design of innovative rotary-wing aircraft configurations that requires a wide number of numerical simulations to reproduce the different attitudes of their flight mission.

4.2. Results and Discussion

The present section reports the main results of the coupled simulations for the complete vehicle, starting from the discussion of the response to the roll maneuver in terms of flight mechanics performance, here intended in terms of bank angle and roll rate for a prescribed aileron deflection.

Figure 15a shows the comparison of the bank angle $\varphi$ evolution during the simulated roll maneuver for the three aircraft configurations tested. In particular, the figure clearly shows that the aerodynamic effects of the tailplanes and rotor wakes changes the slope of the bank angle curve and influences the roll maneuver performance.

Specifically, roll performance appears to degrade when additional contributions to the aerodynamics are considered.

Table 10. Percent difference of time to bank 45° for aircraft C.I and C.II with respect to C.III.

Configuration	$\mathbf{\Delta }$%
C.I	−4.46 %
C.II	−3.33 %

shows the negative percentage differences between the time computed to bank by a 45° angle for aircraft C.I and C.II with respect to the complete aircraft equipped with rotors, i.e., C.III.

The performance reduction due to the tail surfaces amounts to about 1%; it is explained by the change in angle of attack due to the roll rate, which produces a negative roll moment contribution. The performance reduction due to the proprotor aerodynamics, confirmed by the comparison of roll rate evolution ($\dot{\varphi }$) presented in Figure 15b, is related to a backward tilting of the rotor induced by the component of reference velocity associated with roll rate in the rotor disk plane. It is also observed that the aircraft roll causes an opposite variation of the rotors thrust, as shown in Figure 16.

This imbalance between the right and left rotor thrust generates a non-negligible contribution to the yaw moment reaction at the revolute joint that grounds the aircraft. Since the thrust of the left rotor increases and that of the right one decreases, the sign of the yaw moment perturbation is negative.

Figure 17a shows the measured yaw moment reaction at the ground joint. Since the position of this constraint is fixed during the simulation, the purpose of these results is not to establish the actual behavior of the aircraft but to estimate the impact of the different vehicle parts on this quantity.

The analysis shows that tailplanes generate a proverse yaw contribution compared to the simplest C.I. On the contrary, the introduction of the aerodynamic contribution of the proprotors generates an adverse yaw moment. This effect is related to the opposite variation of the rotor thrusts highlighted in Figure 16.

Considering dynamic oscillations of the aircraft during the roll maneuver, a quite complex response is observed for C.III, the full aircraft configuration with rotors. In particular, the spectrogram of the yaw moment time history computed for C.III shown in Figure 17b clearly identifies the correspondence of these oscillations with the multiples of the rotor $N/rev$. Furthermore, it is evident that almost all harmonics increase in power as the simulation advances over time. The development of these aerodynamic states is due to the capability of the aerodynamic solver DUST to capture aerodynamic interactions.

Thus, in order to deeply understand the effect of possible aerodynamic interactions on the maneuver, the effect of the two proprotors on the normal force acting on the wing is investigated. Figure 18 shows the convention adopted for the azimuthal angle of the blade $\psi$.

Figure 19 shows the azimuthal histories of the normal force computed on the left and right side of the wing for the three investigated configurations over a complete rotor revolution at the same bank angle $\varphi =30°$. As one may expected, for aircraft configurations I and II, there are no oscillations due to the aerodynamic loads, as the proprotor’s wake is absent. A slight difference between results obtained for configurations C.I and C.II can be noted on the right wing only, as the introduction of the tail provides an increase of normal aerodynamic load during the maneuver. The blue line curves, corresponding to the wing normal force for C.III, show the three peaks related to the passage of the blades in front of the wing. This comparison shows the capability of the mid-fidelity aerodynamic model to reproduce the detailed effects of proprotor wake interaction on the wing, providing an increase of the aerodynamic load on both sides. This effect is mainly due to the acceleration of the flow impinging the wing caused by the proprotor inflow, which provides a local increase of the wing loading.

Figure 20 shows the span-wise distribution of the normal force acting on the wing for the three simulated aircraft configurations at bank angle $\varphi =30°$. Results for C.III are plotted at three different blade azimuth angles $\psi$. The span-wise wing load comparison confirms the increase of the normal force due to interaction of the proprotor wake. The effect on the aerodynamic loading due to flaperons deflection is apparent, as well. Indeed, the aerodynamic model accurately captures the reduction of wing loading on the left, due to the control surface negative deflection, and the increase of loading on the right due to the control surface positive deflection. Focusing on C.III simulation results, the span-wise load variation slightly changes with blade azimuth angle $\psi$, due to the interaction between the tip vortices released by the proprotor blades and the wing. Indeed, the normal force acting on the wing presents a peak at $\psi =0$, as also shown in the time history of Figure 19.

A more detailed insight into the interactional effects of the proprotors on the wing can be deduced from Figure 21. The figure presents a visualization of the flow field behind the proprotors through Q-criterion iso-surfaces and the distribution of the pressure coefficient $C_P$ over the aircraft surfaces for the three tested configurations. In detail, the flow field representation for C.III clearly shows the interference of the proprotor helical system with the wing. This interaction mainly provides an up-wash due to the blade tip vortex encountering the inboard region of the wing [37] and a consequent local increase of pressure with respect to C.I and C.II, as shown by the larger and more intense pressure region computed over the wing surface for the aircraft configuration equipped with proprotors (see Figure 21c).

Figure 22 compares the pressure coefficient $C_P$ distribution evaluated at the mid-span section of the flaperon. The $C_P$ curves show that interaction with the proprotor wake is responsible for an increase of the suction peak on the upper surface at the leading-edge region of the airfoil. Otherwise, the pressure coefficient distribution is quite similar for all configurations in the aft portion of the airfoil, which corresponds to the deflected flaperon. A higher aerodynamic loading of the complete aircraft wing results in C.III with respect to C.I and C.II due to aerodynamic interaction with the proprotor’s wake.

Finally, the effect of the aerodynamic interference of the wing on the proprotors is investigated. Figure 23 shows the contours of the blade lift coefficient ($C_L$) computed over a complete rotor revolution for different bank angles $\varphi$. The polar plots show that the mid-fidelity aerodynamic model captures quite well the variation of the thrust force distribution in the region corresponding to the passage of the blade in front of the wing. Considering the blade lift coefficient distributions on the proprotors at the same bank angle, a quite similar distribution is observed between the two rotors, except for a small difference in the innermost area of the disk, which is related to the different position of the aileron (positive deflection on the right, negative deflection on the left); see Figure 23a,b referring to the moment in which the ailerons are deflected, and the bank angle is null. As the maneuver progresses and the bank angle increases, the lift distribution on the two proprotors presents more differences. In particular, focusing the attention on the right rotor, it is evident that increasing the bank angle between $\varphi =0°$ (see Figure 23a) and $\varphi =45°$ (see Figure 23g), the lift coefficient positive variation in the outer region of the blade decreases in the range of blade azimuthal angle between $\psi =210°$ and $\psi =270°$, characterized by the passage of the blade in front of the wing. The left rotor behaves in the opposite way showing, as the bank angle increase, a positive variation of the lift coefficient in the outer region of the blade for the same azimuthal blade angle range (see Figure 23b–h). In detail, the magnitude of this lift variation within the maneuver is greater for the left proprotor. This opposite trend of the proprotors local loads reflects the physics of the unbalance of the two proprotors thrust discussed in Figure 16 and the consequent growth of an adverse contribution to the yaw moment.

5. Conclusions

A new aeroelastic numerical tool was presented, obtained by coupling a mid-fidelity aerodynamic solver with a multibody dynamics software. It showed the advantages and limits of the proposed approach for the aeroelastic simulation of rotary-wing vehicles. The resulting solver was thoroughly described, highlighting the key-aspects of the mathematical formulation of both solvers, as well as the implementation details of their coupling. Furthermore, the capability to simulate movable surfaces deflection implemented in the aerodynamic solver was highlighted, due to their importance in simulating maneuvers and control studies that represent essential aspects in the design of innovative air vehicles.

The combination of a multibody structural model with a mid-fidelity numerical representation of the aerodynamics enabled to accurately capture aerodynamic interactional effects that characterize rotary-wing aircraft with limited computational cost compared to conventional, higher-fidelity CSD/CFD tools. The coupled code was validated by comparing its results with numerical and experimental data available in the open literature for both fixed- and rotary-wing problems. In detail, validation tests showed the importance of a more accurate description of the aerodynamics to reproduce the dynamic behavior of a wedged wing in flutter condition, and to improve the performance evaluation of a proprotor in hover.

Finally, a complex operating condition was studied by simulating a complete model of the XV-15 tiltrotor performing a roll maneuver. The discussion of the simulation results highlighted the flexibility of the proposed tool in simulating complex rotary-wing aircraft configurations and its capability to capture the fine details of flow physics and the effects on vehicle performance of the aerodynamic interaction between wakes and bodies that are typical of rotary-wing aircraft.

The outcome of this work opens a novel scenario for the scientific and industrial community. Indeed, the significantly lower computational effort required by the coupled simulations with respect to conventional, higher-fidelity CSD tools coupled with URANS solvers, as well as the quite high level of accuracy that could be obtained by these simulations, indicate that the use a mid-fidelity aerodynamic solver for coupled aerodynamics-multibody dynamics simulations represents a valuable instrument for performing the large number of aeroelastic analyses required in the preliminary design phase of innovative rotary-wing vehicles, not only in the rotorcraft field but also in other domains, for example, wind energy and turbomachinery applications.