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

**Alberto Savino, Alessandro Cocco, Alex Zanotti, Matteo Tugnoli, Pierangelo Masarati and Vincenzo Muscarello**

Dipartimento di Scienze e Tecnologie Aerospaziali, Politecnico di Milano, Via La Masa 34, 20156 Milan, Italy; alberto.savino@polimi.it (A.S.); alessandro.cocco@polimi.it (A.C.); matteo.tugnoli@polimi.it (M.T.); pierangelo.masarati@polimi.it (P.M.); vincenzo.muscarello@polimi.it (V.M.)

**\*** Correspondence: alex.zanotti@polimi.it

**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.

**Keywords:** aeroelasticity; fluid-structure interaction; rotary-wing aerodynamics; multibody dynamics; tiltrotor; computational fluid dynamics; vortex particle method

## **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–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

**Citation:** Savino, A.; Cocco, A.; Zanotti, A.; Tugnoli, M.; Masarati, P.; Muscarello, V. Coupling Mid-Fidelity Aerodynamics and Multibody Dynamics for the Aeroelastic Analysis of Rotary-Wing Vehicles. *Energies* **2021**, *14*, 6979. https:// doi.org/10.3390/en14216979

Academic Editor: Chunhua Liu

Received: 20 September 2021 Accepted: 19 October 2021 Published: 25 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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 highfidelity 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–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–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–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:

$$\mathbf{M(x,t)\dot{x} = p} \tag{1}$$

$$
\dot{\mathbf{p}} = \boldsymbol{\Phi}\_{/\mathbf{x}}^T \boldsymbol{\lambda} + \mathbf{f}\_{\mathrm{i}}(\dot{\mathbf{x}}, \mathbf{x}, t) + \mathbf{f}\_{\mathrm{e}}(\dot{\mathbf{x}}, \mathbf{x}, t) \tag{2}
$$

$$\boldsymbol{\phi}(\mathbf{x}) = \mathbf{0},\tag{3}$$

where **x** is the vector of the kinematic unknowns, **p** is the vector of the momentum unknowns, *λ* is the vector of the Lagrangian multipliers, **M** is a configuration (and possibly time) dependent inertia matrix, **f***<sup>i</sup>* ,**f***<sup>e</sup>* are arbitrary internal and external forces, *φ*(**x**) are the algebraic constraint equations (only holonomic constraints are shown for conciseness), and *φ<sup>T</sup>* /*x* 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 **f***<sup>i</sup>* , 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 (**u***ϕ*) and solenoidal **u***<sup>ψ</sup>* contributions of velocity as **u** = **u***<sup>ϕ</sup>* + **u***ψ*. 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

<sup>198</sup> indices.

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. <sup>175</sup> For a detailed insight on the mathematical formulation of the DUST code, including the <sup>176</sup> implemented governing equations, the reader is referred to a recent work published by some of <sup>177</sup> the present article authors [36].

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]. <sup>178</sup> *2.3. Description of MBDyn-DUST coupling*

## *2.3. Description of MBDyn-DUST Coupling* <sup>179</sup> Communications between DUST and MBDyn are managed by preCICE (Precise Code Interaction

<sup>174</sup> of the Bernoulli theorem, considering the vorticity of the flow.

Version October 8, 2021 submitted to *Energies* 5 of 29

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. <sup>180</sup> Coupling Environment), a coupling library for partitioned multi-physics simulations, originally <sup>181</sup> developed for fluid-structure interaction and conjugate heat transfer simulations. preCICE offers <sup>182</sup> methods for transient equations coupling, communication means, and data mapping schemes. It is <sup>183</sup> written in C++ and offers additional bindings for C, Fortran, Matlab, and Python. preCICE (https: <sup>184</sup> //github.com/precice/) is a free software released under the LGPL3 license. <sup>185</sup> While MBDyn uses its own C++, C, Fortran and Python Application Programming Interface (API) <sup>186</sup> to communicate with external software without any further modification to the C++ source code, no <sup>187</sup> API was already available in DUST. Thus, new Fortran modules collecting all the classes, subroutines <sup>188</sup> and functions required by the adapter for preCICE were implemented to support coupling with DUST. <sup>189</sup> 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. <sup>190</sup> A new adapter was implemented for support the communication of all the kinematic variables <sup>191</sup> (position, orientation, velocity and angular velocity) plus forces and moments acting on the nodes of a <sup>192</sup> MBDyn model exposed through an external structural force element. Figure 1 shows a scheme of the <sup>193</sup> communication and information exchange, managed through the adapters for the two solvers.

**Figure 1.** Scheme of the communication managed through adapters for MBDyn and DUST **Figure 1.** Scheme of the communication managed through adapters for MBDyn and DUST.

<sup>194</sup> The interface between structural and aerodynamic grids is obtained as a weighted average of the <sup>195</sup> distance between the nodes of the two grids and is used for motion interpolation and consistent force <sup>196</sup> and moment reduction. Figure 2 shows the *q* nodes of the structural grid, namely *Qq*. The centers and <sup>197</sup> the vertices of each aerodynamic mesh are respectively *P<sup>e</sup>* and *Pp*, where *e* and *p* are the corresponding 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 *Qq*. The centers and the vertices of each aerodynamic mesh are, respectively, *P<sup>e</sup>* and *Pp*, where *e* and *p* are the corresponding indices.

> The kinematic variables, *φp*, of a point *p* positioned on the aerodynamic surface of a DUST component is evaluated as the following weighted-average,

$$
\phi\_p = \sum\_q w\_{pq} \phi\_{q\prime} \tag{4}
$$

where *φ<sup>q</sup>* is the same kinematic variable associated with the *q*th structural node of the MBDyn model.

Weights *wpq* could be any set of non-negative real numbers, satisfying the normalization conditions

$$\sum\_{q} w\_{pq} = 1 \qquad \forall p\_{\prime} \tag{5}$$

**Figure 2.** Scheme for motion interpolation (left) and force and moment transfer (right). Structural points are represented by red dots. Nodes of the aerodynamic mesh and panels centers are represented 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<sup>p</sup>* − *Qq*). As an example, using

with plain dot and crosses, respectively;

the local coordinates in the reference configuration *rpq*, the norm of these vectors can be defined as:

$$\|(P\_p - Q\_q)\|^2 := r\_{pq}^T \mathbf{W} r\_{pq} \tag{6}$$

where **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

$$\mathbb{E}(P\_p - O)^{\mathcal{g}} = \sum\_{Q \in I\_P} w\_{pq} \left\{ (Q\_q - O)^{\mathcal{g}} + \mathbf{R}\_Q^{r \to \mathcal{g}} (P\_p - Q\_q) \right\},\tag{7}$$

where *Q* ∈ *I<sup>p</sup>* indicates the subset of structural points *Q<sup>q</sup>* that belong to the *I<sup>p</sup>* aerodynamic component, (*Q<sup>q</sup>* − *O*) *g* is the distance from the origin of the *Q<sup>q</sup>* structural point, and **R** *r*→*g Q* (*P<sup>p</sup>* − *Qq*) rotates in the global coordinates the distance between the aerodynamic point and the structural one. Its angular velocity *ω<sup>P</sup>* and velocity *vP*, respectively, are

$$
\omega\_P = \sum\_{Q \in I p} w\_{pq} \,\omega\_{Q \cdot \prime} \quad \text{ $\mathbf{v}\_P = \sum\_{Q \in I p} w\_{pq} \{ \mathbf{v}\_Q + \omega\_Q \times (P - Q) \}$ }.\tag{8}
$$

The aerodynamic forces and moments are evaluated at the evaluation points *P<sup>e</sup>* 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:

$$f\_Q = \sum\_{\varepsilon \in I\_Q} w\_{q\varepsilon} f\_{\varepsilon} \quad \mathfrak{m}\_Q = \sum\_{\varepsilon \in I\_Q} w\_{q\varepsilon} \{ \mathfrak{m}\_{\varepsilon} + (P\_{\varepsilon} - Q\_q) \times f\_{\varepsilon} \} \tag{9}$$

where *e* ∈ *J<sup>Q</sup>* indicates the subset of evaluation points that belong to each sub-component *JQ*, and the weights *wqe* are calculated using Equations (6) and (7), by computing the distance between each structural node and evaluation point.

**Figure 2.** Scheme for motion interpolation (**left**) and force and moment transfer (**right**). Structural points are represented by red dots. Nodes of the aerodynamic mesh and panels centers are represented with plain dot and crosses, respectively.

## 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.

[39].

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. <sup>216</sup> follows. Before the time loop starts, communication is first established between the coupled codes. <sup>217</sup> The time loop starts with the update of DUST's explicit aerodynamic elements (lifting lines and <sup>218</sup> actuator disks). A checkpoint of the exchanged fields is stored, to be reloaded during sub-iterations of <sup>219</sup> preCICE's implicit coupling. DUST then receives the kinematic variables of the structural nodes from <sup>220</sup> the external software, MBDyn, and updates the surfaces of the coupled components and the near-field <sup>221</sup> wake elements. The linear system is then updated and solved, calculating the strengths of the vortexes <sup>222</sup> of the surface panels and vortex lattice elements.

Version October 8, 2021 submitted to *Energies* 7 of 29

<sup>215</sup> in the communication are initialized. Initialization of the Fast Multiple kernels, wake and linear system

**Figure 3.** Flowchart of the implicit communication managed by preCICE between DUST and MBDyn **Figure 3.** Flowchart of the implicit communication managed by preCICE between DUST and MBDyn [39].

<sup>223</sup> The solution of the non-linear lifting line problem follows. The circulation **Γ** is computed by <sup>224</sup> using analytical expression from the Kutta-Joukowski theorem where lift coefficient is evaluated <sup>225</sup> from tabulated sectional data. Once the intensity of the surface singularities has been evaluated, <sup>226</sup> surface pressure distribution and elementary forces and moments are retrieved using the unsteady The solution of the non-linear lifting line problem follows. The circulation **Γ** 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.

<sup>227</sup> Kutta-Joukowski theorem for the vortex lattice and lifting lines elements, and the unsteady Bernoulli <sup>228</sup> 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

<sup>230</sup> aerodynamic and structural meshes and sent to MBDyn. A convergence check on the kinematics

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 *ξ*, and a blending region [−*u*, *u*] introduced to avoid irregularities in the mesh as the surface rotates by an angle *θ*. Moreover, in the 2D modeling, the rotation axis *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 *<sup>H</sup>* and axes **<sup>ˆ</sup>***ξ*, *<sup>η</sup>***<sup>ˆ</sup>** <sup>=</sup> *<sup>h</sup>***<sup>ˆ</sup>** <sup>×</sup> **<sup>ˆ</sup>***ξ*. The position of a point with respect to this reference frame is

$$
\mathbf{r} = \boldsymbol{\xi}\,\mathbf{\hat{f}} + \boldsymbol{\eta}\,\boldsymbol{\eta} + h\,\mathbf{\hat{f}}.\tag{10}
$$

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


$$
\Delta \mathbf{r} = \sin \theta \hat{\mathbf{h}} \times \mathbf{r} + (1 - \cos \theta) \hat{\mathbf{h}} \times \hat{\mathbf{h}} \times \mathbf{r};\tag{11}
$$

3. −*u* ≤ *ξ* ≤ *u*: blending region to avoid irregularities, defined as an arc of a circle whose center is located at point *C* and whose radius is:

**Figure 4.** Scheme of the two-dimensional hinged surface configuration.

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.

**Figure 5.** Hinge reference system for a swept wing.

**Figure 5.** Hinge reference system for a swept wing. The aerodynamic sections that are involved in the deflection of the control surface are, thus, the ones that satisfy the condition *y*(*A*) < *y*(*P*) < *y*(*B*), where *y*(*P*) is the ordinate of the *P<sup>i</sup>* − *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 aerodynamic sections that are involved in the deflection of the control surface are, thus, the ones that satisfy the condition *y*(*A*) < *y*(*P*) < *y*(*B*), where *y*(*P*) is the ordinate of the *P<sup>i</sup>* − *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.

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 **hˆ** = (*B* − *A*) k*B* − *A*k . Each point of the movable surface is 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 **hˆ** = (*B* − *A*) k*B* − *A*k . 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 *wph* are evaluated using only the *h* components of the vectors connecting the control surface points to the hinge nodes.

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 *wph* are evaluated using only the *h* components of the vectors connecting the control surface points to the hinge nodes. The weights *wph* are then combined to the interpolations weight *w* described in The weights *wph* 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 *wph* are specifically used to impose the rigid deflection of the movable surface.

**Figure 6.** Control surface deflection for a deformable wing.

## **3. Validation of the Coupled MBDyn-DUST Software**

## *3.1. Goland's Wing*

*3.1. Goland's Wing*

**Figure 6.** Control surface deflection for a deformable wing. **3. Validation of the Coupled MBDyn-DUST Software** 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* ≈ 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*<sup>∞</sup> is the free-stream

This section presents the application of the coupled software to Goland's wing [56],

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*<sup>∞</sup> is the free-stream

velocity. All the relevant geometrical and structural properties are reported in Table 1. They have been obtained from Reference [57].

**Table 1.** Goland's wing properties and flight condition [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* <sup>0</sup> 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.

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


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® CoreTM 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.

**Figure 8.** Time history of the Goland's wing-tip deflection evaluated with surface panels aerodynamic mesh at three different wind speeds.

Table 3 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.

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


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.

**Figure 9.** Frequency and damping versus velocity for Goland's wing. Coupled simulations results (VL and SP mesh) and MBDyn results with 2D Strip Theory aerodynamics.

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:


Version October 8, 2021 submitted to *Energies* 13 of 29

**Figure 10.** Layout of the XV-15 proprotor control chain.

**Figure 11.** Flowchart of the MBDyn model of the XV-15 proprotor, particularly showing the individual blade pitch control system components and their connections for the dual control path. **Figure 11.** Flowchart of the MBDyn model of the XV-15 proprotor, particularly showing the individual blade pitch control system components and their connections for the dual control path.


All rotor data were taken from the original CAMRAD II model presented in Reference [62], here reported for clarity in Table 4. Table 5 presents the blade span-wise airfoil distribution, from Reference [64].

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


**Table 5.** XV-15 Blade airfoil distribution [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 *CT*/*σ* versus torque coefficient *CQ*/*σ* 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:


**Figure 12.** Comparison of results obtained for the XV-15 full-scale single proprotor in hover: (**a**) torque coefficient *CQ*/*σ* versus thrust coefficient *CT*/*σ*; (**b**) figure of merit (FM) versus thrust coefficient *CT*/*σ*.

## **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.


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

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.

**Figure 13.** Aerodynamic mesh of the XV-15 tiltrotor model.

The dynamic model includes:


The inertial structural data used for the multibody model were taken from the CAM-RAD II and NASTRAN models presented in Reference [62]. Table 7 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.


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

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 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.

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


**Figure 14.** Geometry of aircraft configurations tested for roll maneuver coupled simulations.

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.

During the roll maneuver simulations, the flaperons are deflected according to a step function from *δ<sup>a</sup>* = 0° to *δ<sup>a</sup>* = ±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.


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

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 *φ* 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.

**Figure 15.** Comparison of: (**a**) the computed bank angle evolution and (**b**) roll rate evolution for the tested configurations.

Specifically, roll performance appears to degrade when additional contributions to the aerodynamics are considered. Table 10 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.

**Table 10.** Percent difference of time to bank 45° for aircraft C.I and C.II with respect to 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 (*φ*˙) 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.

**Figure 16.** Rotor thrust with respect to bank angle during the roll maneuver.

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.

**Figure 17.** Yaw moment during the roll maneuver. (**a**) Yawing moment reaction evolution during roll. (**b**) Spectrogram of yaw moment time history for C.III.

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 <sup>501</sup> aerodynamic interactions.

504

Version October 8, 2021 submitted to *Energies* 20 of 29

<sup>499</sup> that almost all harmonics increase in power as the simulation advances over time. The development <sup>500</sup> of these aerodynamic states is due to the capability of the aerodynamic solver DUST to capture

<sup>502</sup> Thus, in order to deeply understand the effect of possible aerodynamic interactions on the

Figure 18 shows the convention adopted for the azimuthal angle of the blade *ψ*.

is investigated. Figure 18 shows the convention adopted for the azimuthal angle of the blade *ψ*.

$$\textbf{Figure 18.}\text{ Convention adopted for the rotor blade azimuthal angle }\psi\text{.}$$

**Figure 18.** Convention adopted for the rotor blade azimuthal angle *ψ*. <sup>505</sup> Figure 19(b) shows the azimuthal histories of the normal force computed on the left and right <sup>506</sup> side of the wing for the three investigated configurations over a complete rotor revolution at the same <sup>507</sup> bank angle *φ* = 30°. As one may expected, for aircraft configurations I and II there are no oscillations <sup>508</sup> due to the aerodynamic loads, as the proprotors wake is absent. A slight difference between results <sup>509</sup> obtained for configurations C.I and C.II can be noted on the right wing only, as the introduction of <sup>510</sup> the tail provides an increase of normal aerodynamic load during the maneuver. The blue line curves 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 *φ* = 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 19.** Azimuthal history of wing normal force during a rotor revolution at bank angle *φ* = 30°, (**a**) left side and (**b**) right side.

Figure 20 shows the span-wise distribution of the normal force acting on the wing for the three simulated aircraft configurations at bank angle *φ* = 30°. Results for C.III are plotted at three different blade azimuth angles *ψ*. 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 *ψ*, 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 *ψ* = 0°, as also shown in the time history of Figure 19.

**Figure 20.** Span-wise normal force distribution on the wing at bank angle *φ* = 30°.

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<sup>P</sup>* 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 21.** Coupled simulation results for the three configurations at bank angle *φ* = 30°: pressure coefficient *C<sup>P</sup>* contours over the aircraft surface for the three configurations and iso-surface of Q-Criterion (*Q* = 1200) for flow field description for C.III.

Figure 22 compares the pressure coefficient *C<sup>P</sup>* distribution evaluated at the mid-span section of the flaperon. The *C<sup>P</sup>* 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.

**Figure 22.** *C<sup>P</sup>* distribution at flaperon mid-span for the three configurations at bank angle *φ* = 30°. (**a**) *C<sup>P</sup>* left wing; (**b**) *C<sup>P</sup>* right wing.

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 (*CL*) computed over a complete rotor revolution for different bank angles *φ*. 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 *φ* = 0° (see Figure 23a) and *φ* = 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 *ψ* = 210° and *ψ* = 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.

**Figure 23.** Lift coefficient disks at various bank angles *φ*, front view as shown in Figure 18.

## **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.

**Author Contributions:** Conceptualization, A.S., A.C.; methodology, A.S., A.C., A.Z.; software, A.S., A.C., M.T., P.M.; validation, A.S., A.C.; formal analysis, A.S., A.C.; investigation, A.S., A.C.; resources, A.S., A.C.; data curation, A.S., A.C.; writing—original draft preparation, A.Z., A.S., A.C.; writing review and editing, A.Z., A.S., A.C., P.M., V.M.; supervision, A.Z., V.M., P.M.; project administration, A.Z., V.M., P.M.; funding acquisition, P.M., V.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received funding from the Clean Sky 2—H2020 Framework Program, under grant agreement N. 885971 (the FORMOSA project) and under grant agreement N. 863418 (the ATTILA project).

**Conflicts of Interest:** The authors declare no conflict of interest.

## **Abbreviations**

The following nomenclature and abbreviations are used in this manuscript:



## **References**

