A non-linear unsteady vortex-lattice method for rotorcraft applications

Abstract

The present work aims to extend the capabilities of DUST, a mid-fidelity aerodynamic solver developed at Politecnico di Milano, for the aerodynamic simulation of rotorcraft applications. With this aim, a numerical element was implemented in the solver obtained by a coupling between the potential unsteady vortex lattice method and viscous aerodynamic data of aerofoil sections available from two-dimensional high-fidelity computational fluid dynamics (CFD) simulations or experimental wind-tunnel tests. The paper describes the mathematical formulation of the method as well as a validation of the implementation performed by comparison with both high-fidelity CFD simulation results and experimental data obtained over aerodynamics and aeroelastic fixed-wing benchmarks. Then, the method was used for the evaluation of the aerodynamic performance of two rotorcraft test cases, i.e. the full-scale proprotor of the XV-15 tiltrotor operating in different flight conditions and two propellers in tandem with overlapping disks. Simulation results comparison with high-fidelity CFD and data from wind tunnel tests highlighted the potentialities and advantages of the implemented approach to be used for the design and investigation of rotorcraft configurations characterised by consistent viscosity effects.

Nomenclature

${C_L}$

lift coefficient, $ = L/\!\left( {0.5{ho} U_\infty ^2S} \right)$

${C_P}$

power coefficient, $ = P/\!\left( {{ho} {n^3}{D^5}} \right)$

${C_Q}$

torque coefficient $ = Q/\!\left( {{ho} {n^2}{D^5}} \right)$

${C_T}$

thrust coefficient $ = T/\!\left( {{ho} {n^2}{D^4}} \right)$

$D$

propeller diameter [m]

$f$

torsional mode frequency [Hz]

$J$

advance ratio $ = {V_\infty }/\!\left( {nD} \right)$

${L_x}$

longitudinal distance between the propeller axis [m]

${L_y}$

lateral distance between the propeller axis [m]

LL

lifting line

$n$

rotational speed [rad/s]

NL-VL

non-linear vortex lattice

$P$

propeller power [W]

$Q$

rotor torque [Nm]

$R$

propeller blade radius [m]

SP

surface panel

$T$

propeller thrust [N]

$u$

freestream velocity component [m/s]

${U_\infty }$

wind tunnel freestream velocity [m/s]

VL

vortex lattice

Greek symbol

${{\eta}}$

rotor efficiency

${{{ho} }}$

air density [kg/m $^3$ ]

${{\xi}}$

torsional mode damping

1.0 Introduction

The study of the complex aerodynamics characterising the novel concept of electric take-off and landing (eVTOL) aircraft for advanced air mobility (AAM) represents a challenge for computational fluid dynamics (CFD) tools. These kinds of flying machines experience a wide variety of flight conditions in which the rotor works with a high disk loading in helicopter mode, while a lower disk loading is experienced in aircraft mode, where the rotor behaves as a propeller. Currently, to obtain appreciable performance for both these flight conditions, the optimal design of eVTOL rotors leads to the use of highly twisted blades. Consequently, blades’ aerofoil sections experience a wide range of angles of attack that could overcome the stall condition, particularly in hovering. Moreover, the architectures of these innovative vehicles are characterised by multi-rotor configurations that lead to complex interactional mechanisms between propeller wakes in the different stages of their flight envelope. Thus, aerodynamic simulations of innovative rotorcraft configurations as tiltrotors are performed using solvers based on a finite volume implementation of the Reynolds Averaged Navier-Stokes (RANS) equations [1, 2]. Nevertheless, time-accurate RANS simulations of a complete vertical take-off and landing (VTOL) aircraft are still time-consuming and require a huge amount of computational resources. Thus, high-fidelity CFD tools are usually devoted to a limited number of detailed analyses, but they are still not suitable to be used in the design process of a new VTOL aircraft. In the preliminary design phase of innovative rotorcraft concepts such as eVTOLs, a key aspect is to obtain a fast and accurate evaluation of the loads acting on rotors and lifting surfaces. Indeed, as previously cited, the flight mission of these novel aircraft is characterised by different configurations, from hover to conversion phase and aircraft mode. Thus, a huge number of simulations are required during their design phase to cover the entire aircraft flight envelope. Consequently, in recent years, the interest of industrial and scientific communities concerning the use of mid-fidelity aerodynamic solvers, based on the vortex-particle method (VPM) for wake modeling [3, 4], has grown in the field of rotorcraft simulations. This interest was finalised to obtain fast and accurate numerical tools to be used in the preliminary design phase of novel rotorcraft architectures. Due to the lower computational costs required by mid-fidelity numerical methods with respect to high-fidelity CFD simulations based on Navier-Stokes solvers, these tools represent a quite interesting solution to be used for the development of innovative rotorcraft configurations as well as to investigate interactional aerodynamics mechanisms typical of these vehicles.

Just to cite a few examples of mid-fidelity aerodynamics solvers, German Aerospace Center (DLR) developed Unsteady Panel Method (UPM) [5] an unsteady panel and free-wake code originally intended for aeroacoustic simulations of helicopters but recently applied on arbitrary complex configurations as compound rotorcraft [6]. Moreover, the National Technical University of Athens (NTUA) developed GENeral Unsteady Vortex Particle (GENUVP), a software based on a panel method coupled with a VPM solver used for both aerodynamic and aeroacoustic simulations of rotorcraft [7]. The University of Glasgow developed an unsteady aerodynamic analysis method based on VPM recently applied for the investigation of the complex wake of coaxial rotors in the work by Tan et al. [8], while a vortex-based approach coupled with a viscous boundary model was used by the same authors to investigate complex rotorcraft-to-rotorcraft interference problems that occur during shipboard operations [9]. Brigham Young University developed FlowUnsteady, an open-source variable-fidelity framework for unsteady aerodynamics and aeroacoustics based on the reformulated vortex particle method (rVPM), widely employed for the investigation of rotor-on-rotor aerodynamic interactions and complete eVTOLs configurations [10].

The mid-fidelity open-source software DUST has been developed by Politecnico di Milano since 2017 for the simulation of the interactional aerodynamics of rotorcraft and unconventional aircraft configurations [11]. The code is released under the open-source Massachusetts Institute of Technology (MIT) license. The capabilities of the code have been quite extended in recent years, and DUST has been also coupled to the open-source multibody solver MBDyn [12], also developed at Politecnico di Milano, enabling to perform aeroelastic analysis of complete rotorcraft configurations. The solver has reached a mature state and was robustly validated against experimental and high-fidelity CFD codes by simulating complex rotorcraft configurations as complete eVTOL and tiltrotor aircraft [13, 14], thus is currently widely used by both academic and industrial communities. The mathematical formulation of DUST relies on an integral boundary element formulation of the aerodynamic problem and vortex-particle modelling [3, 4] of the wakes. This choice naturally fits the Helmholtz decomposition of the velocity field from a mathematical point of view and avoids the numerical instabilities occurring with connected models of the wake. A model can be composed of several components, connected to user-defined reference frames, whose position and motion can be defined in a hierarchical way. Different aerodynamic elements allow for different levels of fidelity in the model, ranging from lifting line elements to zero-thickness lifting surfaces and surface panels for thick solid bodies.

The present work describes a new numerical element aimed to extend the core capabilities of the methods already implemented in DUST for the simulation of rotorcraft applications characterised by important viscous effects. In particular, the numerical element is based on a coupling between the potential unsteady vortex lattice method and viscous aerodynamic data of aerofoil sections available from two-dimensional high-fidelity CFD simulations (RANS) or experimental data. This approach was successfully used in recent literature by several authors with applications on both fixed and rotary wings. To cite a few examples, Dos Santos and Marques [15] propose a pre- and poststall correction based on the Kirchhoff flow approach to non-linear lift applied to the vortex lattice method. The proposed methodology was applied to different wing geometries, including changes in the taper ratio, and showed a good matching with experimental and high-fidelity CFD data for the lift coefficient evaluation. Further application to fixed wings is provided by Parenteau et al. [16] that proposed coupling of classical vortex lattice method with viscous corrections based on sectional RANS computations to evaluate the aerodynamic performance of aircraft wings in transonic cruise and take-off/landing. Mukherjee and Gopalarathnam [17] proposed a novel scheme suitable for implementation in lifting-line and vortex-lattice methods based on an iterative decambering approach for the prediction of poststall characteristics of wings using known section data. This new scheme, differing from how the residual for the Newton iteration is computed, was found to be more robust at achieving convergence and was validated by comparison with experimental data over fixed wing application in stall behaviour. Now considering rotary wing applications, Lee and Lee [18] used a non-linear vortex lattice method for the evaluation of the aerodynamic performance and wake structures of a wind turbine. In particular, this method showed the capability to handle the non-linear stall and post-stall behaviours of blades by coupling vortex lattice elements with experimental look-up tables of blade sectional aerodynamic coefficients.

In the present paper, the description of the mathematical formulation of the non-linear vortex lattice method (NL-VL) implemented in DUST is firstly given, followed by a validation of the implementation based on comparison with high-fidelity CFD and experimental data obtained over simple fixed-wing aerodynamic and aeroelastic benchmarks described in the literature. Then, the paper describes the method’s performance for calculating the aerodynamic performance of two different rotorcraft applications. The first consisted of the study of the full-scale proprotor of the XV-15 tiltrotor operating in hover condition, forward flight in helicopter mode, and aircraft mode. The second is the study of the aerodynamic interaction between two propellers in tandem with overlapping rotor disks, a typical configuration occurring in eVTOL aircraft. Thus, the main goal of the paper is to outline the advantages and limits of such an approach with respect to high-fidelity CFD and already implemented methods in DUST, such as lifting lines and surface panels for the calculation of aerodynamic performance and flow physics of complex rotary-wing applications typical of tiltrotors and eVTOLs.

2.0 Non-linear vortex lattice method

As previously stated, the DUST solver implemented different classical aerodynamic elements, i.e. lifting lines, vortex lattices and surface panels. The lifting line (LL) element implemented in DUST is a 1-D model of thin, slender lifting bodies, whose sectional aerodynamic coefficients of lift, drag and pitching moment, i.e. ${C_l}$ , ${C_d}$ , ${C_m}$ , are provided as a function of the local angle-of-attack $\alpha $ , the local Reynolds number $Re$ and the local Mach number $M$ . This allows taking into account the aerofoil camber and thickness and both viscous and compressibility effects on aerodynamic load computation. However, LL suffers numerical instabilities when used to model low-aspect-ratio wings. These instabilities may occur particularly for aeroelastic problems, where a tight coupling between the aerodynamic and the structural solver is used for the analysis. Moreover, this 1-D element is particularly accurate for the simulations of slender bodies as blades but can lose accuracy for the investigation of low aspect ratio wings typical of tiltrotor and eVTOLs architectures. The surface panels (SP) are formulated as a Morino-like problem [19] in implicit form. With this element, the real shape of the component can be described by applying the non-penetration boundary condition in the physical position. However, this element is not able to capture the nonlinear effect on the aerodynamic loads.

The vortex lattice (VL) method provides the aerodynamic elements for the discrete representation of the mean surface of thin lifting bodies, modeled as a sheet of vortex rings of intensity ${\rm{\Gamma }}$ equivalent to a piecewise-uniform surface doublet distribution. The compressibility effects are taken into consideration by applying a Prandtl-Glauert correction on the loads, while the element is not able to capture the non-linear behaviour of the aerodynamic loads. On the other hand, this element relies on a fully implicit numerical scheme thus providing a high degree of numerical stability. More details about the mathematical formulation of these elements can be found in Ref. [11].

In order to overcome the previously mentioned limitations of these numerical elements, a new non-linear vortex lattice element was developed and implemented in DUST by joining the capabilities of lifting lines and vortex lattices. The description of the mathematical formulation of this new element is provided in the following section.

2.1 Mathematical formulation

In the first step, the potential solution of the linear vortex lattice is calculated, by solving the linear system that imposes the non-penetration boundary condition on the vortex lattice panel as follows:

(1) \begin{align}{\hat{\textbf n}} \cdot {{\textbf{u}}_\phi } = {\hat{\textbf n}} \cdot \left( {{{\textbf{u}}_b} - {{\textbf{U}}_\infty } - {{\textbf{u}}_\psi }} \right),\end{align}

where ${{\textbf{u}}_\phi }$ is the potential velocity, ${{\textbf{u}}_b}$ is the body velocity, ${{\textbf{U}}_\infty }$ is the free-stream velocity ${{\textbf{u}}_\psi }$ is the rotational perturbation velocity. In the vortex lattice method, a surface is discretised as a sheet of vortex wing with intensity ${{\rm{\Gamma }}_{{i_v}}}$ , equivalent to a piecewise-uniform surface doublet distribution with the same intensity ${\mu _{{i_v}}}$ . The boundary condition is written for each panel collocation point and the problem can be rewritten in a linear system form as

(2) \begin{align}{\textbf{A}}{\rm{\Gamma }} = {\rm{RHS}},\end{align}

where each term ${a_{ij}}$ of the matrix ${\textbf{A}}$ is defined as the velocity component induced by unit strength ${{\rm{\Gamma }}_i}$ , normal to the surface at collocation point $i$ . Then, a fixed-point iterative problem is solved by imposing a convergence between the lift coefficient derived from the aerodynamic 2D look-up table and the one calculated from the stripe intensity. The total velocity $\boldsymbol{u}_{{v_{nl}}}^i$ acting on the control point (centre) at time $t$ and iteration $i$ of each stripe can be calculated as follows:

(3) \begin{align}\boldsymbol{u}_{{v_{nl}}}^i = {\boldsymbol{U}_\infty } - {\boldsymbol{u}_{{v_{nl}},b}} + {\boldsymbol{u}_{tot\backslash {v_{nl}}}} + \boldsymbol{u}_{{v_{nl}}}^{i - 1}\end{align}

where ${\boldsymbol{U}_\infty }$ is the free stream velocity, ${\boldsymbol{u}_{{v_{nl}},b}}$ is the body velocity, ${\boldsymbol{u}_{tot\backslash {v_{nl}}}}$ is the velocity induced by all elements (potential wake, vortical particles, and lifting body except the non-linear vortex lattice elements), $\boldsymbol{u}_{{v_{nl}}}^{i - 1}$ is the induced velocity of the stripe on the other stripes which is calculated as the velocity induced by a surface doublet whose vertices are the leading and trailing edge nodes and the intensity is obtained as follows:

(4) \begin{align}{\rm{\Gamma }}_{{v_{nl}}}^{i - 1} = d{\rm{\Gamma }}_{{v_{nl}}}^1 + \mathop {\mathop \sum \limits_{k = 2} }\limits^N d{\rm{\Gamma }}_{{v_{nl}}}^k - d{\rm{\Gamma }}_{{v_{nl}}}^{k - 1}\end{align}

where $d{\rm{\Gamma }}_{{v_{nl}}}^k$ is the intensity panel $k$ -th of the stripe calculated at iteration $i - 1$ , and $N$ is the number of chordwise elements. Then, the velocity is projected on the plane formed by the stripe normal $\hat{\boldsymbol{n}}$ and chord $\hat{\boldsymbol{t}}$ :

(5) \begin{align}\boldsymbol{u}_{2D}^i = \left( {\boldsymbol{u}_{{v_{nl}}}^i \cdot \hat {\boldsymbol{n}}} \right)\hat{\boldsymbol{n}} + \left( {\boldsymbol{u}_{{v_{nl}}}^i \cdot \hat{\boldsymbol{t}}} \right)\hat{\boldsymbol{t}}\end{align}

Angle-of-attack $\alpha $ is calculated as:

(6) \begin{align}{\alpha ^i} = {\tilde \alpha ^i} - \alpha _{2D}^i\end{align}

where ${\tilde \alpha ^i}$ is calculated as:

(7) \begin{align}{\tilde \alpha ^i} = {\rm{atan}}\!\left( {\frac{{\boldsymbol{u}_{2D}^i \cdot \hat{\boldsymbol{n}}}}{{\boldsymbol{u}_{2D}^i \cdot \hat{\boldsymbol{t}}}}} \right)\end{align}

This angle must be corrected by taking out the velocity induced by the stripe itself. According to Ref. [20], this correction is equal to:

(8) \begin{align}\alpha _{2D}^i = \frac{{2 \cdot {\rm{\Gamma }}_l^{i - 1}}}{{{C_{L/\alpha }}{\rm{\;}}c\left| {\boldsymbol{u}_{2D}^i} \right|}}\end{align}

where ${C_{L/\alpha }}$ is the slope of the lift curve in the linear range obtained from the tabulated data, and $c$ is the stripe chord. For each iteration, the vortex lattice linear system is solved by modifying only the right-hand side and taking as a result the updated circulations. To numerically stabilise the method two kinds of relaxation are available:

• • Constant relaxation, where the updated right-hand side of the system is

  (9) \begin{align}{\rm{RH}}{{\rm{S}}_i} = {\rm{RH}}{{\rm{S}}_{i - 1}} + \alpha \cdot {\textbf{r}}\end{align}
  with $\alpha $ being the constant relaxation factor and ${\textbf{r}}$ is the residual vector defined as the difference between the lift coefficient solved from the linear system $C_L^{{\rm{inv}}}$ and the one obtained for the aerodynamic 2D look-up table at iteration $i$ $C_L^{{\rm{visc}}}$ :
  (10) \begin{align}{r_s} = \left( {C_L^{{\rm{inv}}} - C_L^{{\rm{visc}}}} \right) \cdot {\rm{sign}}\!\left( {C_L^{{\rm{inv}}} - C_L^{{\rm{visc}}}} \right)\end{align}

• • Aitken acceleration [21], where the updated right-hand side of the system is

  (11) \begin{align}{\alpha _i} = - {\alpha _{i - 1}}\frac{{{{\textbf{r}}_{i - 1}} \cdot {\rm{\Delta }}{\textbf{r}_i}}}{{{\rm{\Delta }}{{\textbf{r}}_i} \cdot {\rm{\Delta }}{{\textbf{r}}_i}}}\end{align}
  with ${\alpha _i}$ being the Aitken relaxation factor for iteration $i$ , ${\textbf{r}}$ is the residual vector and ${\rm{\Delta }}{{\textbf{r}}_i}$ is the difference between the residuals at iteration $i$ and $i - 1$ . Then the ${\rm{RHS}}$ is updated as described in Equation (9).

For non-linear vortex lattice elements, aerodynamic loads are calculated using the unsteady formulation of the Kutta-Joukowsky theorem, as for the linear case plus a viscous correction on the drag term:

(12) \begin{align}\boldsymbol{d}{\boldsymbol{F}_{i \in s}} = \boldsymbol{dF}_{i \in s}^{sys} + \frac{1}{2}{{ho} _\infty }({\boldsymbol{U}^2}\!\left( {\boldsymbol{r},t} \right)\!d{A_{i \in s}}c_d^s\!\left( {{\rm{cos}}\!\left( \alpha \right)\!\hat{\boldsymbol{t}} + {\rm{sin}}\!\left( \alpha \right)\!\hat{\boldsymbol{n}}} \right)\end{align}

where $i \in s$ stands for a panel $i$ that belongs to the stripe $s$ and $\boldsymbol{dF}_{i \in s}^{sys}$ is the force calculated from the system.

2.1.1 Coupled governing equations

The potential contribution of the velocity field results from the free-stream velocity’s superposition and the induced velocity from the doublet and source singularities, modeling the bound vorticity on the surface of the investigated bodies. Given the free stream velocity, the rotational component of the velocity field, the state of the wake, and the motion of the solid components of the model, the intensity of the surface doublets are the actual unknowns of the potential problem obtained from the surface panel discretisation, lifting line discretisation and the subdividing vortex lattice in the linear $\left( {{v_l}} \right)$ and non-linear case $\left( {{v_{nl}}} \right)$ :

(13) \begin{align}\left\{\begin{array}{l} \left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} \textbf{A}_{s,s} & \textbf{A}_{s,v} & \textbf{A}_{s,v_{nl}} & \textbf{A}_{s,l} \\[5pt] \textbf{C}_{v,s} & \textbf{C}_{v,v} & \textbf{C}_{v,v_{nl}} & \textbf{C}_{v,l} \\[5pt] \textbf{C}_{v_{nl},s} & \textbf{C}_{v_{nl}, v} & \textbf{C}_{v_{nl}, v_{nl}} & \textbf{C}_{v_{nl}, l} \end{array} \right] \left\{ \begin{array}{c} \boldsymbol{\Gamma}_s \\[5pt] \boldsymbol{\Gamma}_{vl} \\[5pt] \boldsymbol{\Gamma}_{nl} \end{array}\right\} = -\left[\begin{array}{c@{\quad}c@{\quad}c} \textbf{B}_{s, s} & \textbf{0}_{s, v_{l}} & \textbf{0}_{s, v_{nl}} \\[5pt] \textbf{D}_{v, s} & -\textbf{I}_{v, v} & \textbf{0}_{v, v_{nl}} \\[5pt] \textbf{D}_{v_{nl}, s} & \textbf{0}_{v_{nl}, v} & -\textbf{I}_{v_{nl}, v_{nl}}\\[5pt] \end{array}\right] \left\{ \begin{array}{c} \sigma_s \\[5pt] \sigma_{vl} \\[5pt] \sigma_{nl} \end{array}\right\} -\left[\begin{array}{c} \textbf{A}_{sw} \\[5pt] \textbf{C}_{v_lw} \\[5pt] \textbf{C}_{v_{nl}w} \end{array}\right] \boldsymbol{\Gamma}_w \\[5pt] \boldsymbol{\Gamma}_l = \boldsymbol{f}_l ( \boldsymbol{\Gamma}, \boldsymbol{\sigma}, \boldsymbol{\Gamma}_w )\end{array}\right.\end{align}

where the vector ${{\boldsymbol{\Gamma }}_w}$ collects the intensity of the wake doublets, ${\boldsymbol{\Gamma }} = \left( {{{\boldsymbol{\Gamma }}_s},{{\boldsymbol{\Gamma }}_{{v_l}}},{{\boldsymbol{\Gamma }}_{{v_{nl}}}},{{\boldsymbol{\Gamma }}_l}} \right)$ collects the intensity of the body doublets on the surface panel, the vortex lattice and the lifting line elements, the vector $\boldsymbol\sigma = \left( {{\boldsymbol\sigma _s},{\boldsymbol\sigma _{{v_l}}},{\boldsymbol\sigma _{{v_{nl}}}}} \right)$ contains the unperturbed relative normal velocity at the collocation points of the surface panels and the vortex lattice elements, ${\textbf{0}_{sv}}$ is a null matrix ${N_s} \times {N_v}$ , ${{\textbf{I}}_{vv}}$ is the identity matrix ${N_v} \times {N_v}$ . The non-linear vortex lattice is considered implicit in the sense that its solution is obtained by iteratively solving the linear system that contains all the elements except the lifting line, whereas the solution of the lifting line problem is decoupled from the other element type.

2.2 Validation of the method over fixed-wing test cases

In order to validate the implementation of the numerical element, this approach was tested on fixed-wing benchmark cases available from literature related to aerodynamic and aeroelastic problems.

Firstly, a simple finite wing benchmark case was considered consisting of a rectangular wing with an aspect ratio equal to 12, with a constant chord length equal to 1 m and a constant NACA 4415 aerofoil section. Figure 1 shows the comparison between the wing lift coefficients computed by DUST and high-fidelity CFD simulations results [22]. In particular, DUST simulation results are related to the use of both the classical LL and VL elements and the non-linear vortex lattice method (NL-VL).

Figure 1. Comparison of the NACA 4415 wing lift coefficient ( ${C_L}$ ) as a function of angle-of-attack: high-fidelity CFD simulation results from Ref. [22] (CFD), (DUST VL) LL and NL-VL simulations results. Mach number = 0.2.

The effectiveness of the non-linear vortex-lattice method is evident in the linear range as well as in the post-stall region, as shown by the remarkable agreement with high-fidelity CFD simulations. In particular, the NL-VL results are quite similar to the ones obtained using the DUST LL approach considering viscous contributions for airloads evaluation. This result confirms the correctness of the implementation based on the integration of the viscous aerodynamic 2D data and shows the capabilities of the method to capture stalled and post-stalled conditions with respect to the purely potential VL model. In addition, thanks to the possibility of the NL-VL elements to evaluate the chordwise pressure distribution, Fig. 2 shows the ${\rm{\Delta }}{C_p}$ distribution evaluated between upper and lower surface of the NACA4415 aerofoil by DUST simulations compared to high-fidelity CFD simulations results from [22] for a nearly stalled and post-stalled condition.

Figure 2. Comparison of ${\rm{\Delta }}{C_p}$ distribution evaluated between upper and lower surface of the NACA4415 aerofoil for high-fidelity CFD simulations from Ref. [22] (CFD) and DUST NL-VL simulations. Mach number = 0.2.

This comparison shows an accuracy similar to CFD of DUST NL-VL method in terms of chordwise pressure difference over the aerofoil for the nearly stalled condition characterised by a low degree of flow separation. On the other hand, due to the limitation of the VL approach for highly separated flow conditions, pressure distribution comparison at very high angle-of-attack shows higher discrepancy with respect to CFD.

Secondly, the Lovell [23] wing test case was considered for a further validation of the DUST methodology. This test case consists in a 30.5 $^{\circ }$ leading edge swept wing (see Fig. 3) used for wind tunnel investigation of the effect of flap deflections. The geometry is simple, with a constant aerofoil geometry along the span and experimental results are available for the isolated wing with different flap deflections. DUST numerical models were built using 50 spanwise elements for both LL and NL-VL while 5 elements in chord were used for the NL-VL model. DUST simulations results are compared to VLM coupled with RANS solutions [16] and experimental data for lift, drag and pitching moment coefficients, see Figs. 4, 5, and 6. Results comparison showed that DUST models reproduce the behaviour of the lift and polar curves obtained by Parenteau et al. [16] with a similar numerical model, thus further validating the implementation of NL-VL model. In particular, the comparison of the ${C_L} - {C_M}$ curve shows a better agreement of the solution obtained with DUST with respect to experimental data.

Figure 3. DUST simulation results of the Lovell wing test case [23]: wake visualisation and sectional load distribution.

Figure 4. Lovell clean configuration: CL. Comparison between experimental results, Parenteau et al. VLM (2D RANS), DUST LL, and DUST NLVL.

Figure 5. Lovell clean configuration: CD. Comparison between experimental results, Parenteau et al. VLM (2D RANS), DUST LL, and DUST NLVL.

Figure 6. Lovell clean configuration: CM. Comparison between experimental results, Parenteau et al. VLM (2D RANS), DUST LL, and DUST NLVL.

Thirdly, an aeroelastic benchmark test was considered, i.e. the Goland’s wing [24]. This test case was widely used in literature to validate aeroelastic codes. The study of this low-aspect-ratio wedged wing with aspect ratio $ \approx 3.33$ is also of interest because it highlights the effect of using a 2D or a 3D aerodynamic model on flutter computation. This aeroelastic model was investigated in a previous work by Savino et al. [12] using DUST aerodynamic models built with classical VL and SP elements coupled to the multibody solver MBDyn code. Indeed, as a matter of fact, the use of LL approach for this kind of coupled simulation failed due to the implicit numerical scheme that characterises this element. The present activity aimed to extend the DUST-MBDyn capabilities by including the use of the non-linear vortex lattice method for the evaluation of wing flutter. In particular, to study the flutter instability, a non-zero angle-of-attack of 0.05° was introduced as a perturbation, as done by Murua et al. [25]. The frequency and damping of the model response were identified from the time history of the wing-tip deflection using the matrix pencil estimation (MPE) method [26]. Table 1 reports a comparison of flutter speed and frequency computed by several authors for this test case. Results from the coupled DUST-MBDyn code are in quite good agreement with those obtained by similar codes using 3D aerodynamic models [25, 27, 28], while a certain difference of flutter speed occurs with respect to calculations obtained with 2D aerofoil models. In particular, results obtained using the NL-VL validated the suitability of this approach to investigate aeroelastic problems. In detail, the discrepancy with results obtained with the same multibody structural model MBDyn, but using its built-in aerodynamic model based on two-dimensional unsteady strip theory, indicates the higher capability of the coupled code based on surface elements as SP and VL for the investigation of aeroelastic problems over low aspect ratio wings.

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

Further validation of the coupled code DUST-MBDyn approach based on the use of different aerodynamic models including the NL-VL approach is provided in Figs. 7 and 8 showing the computed frequency and damping of the first beam torsional mode of the wing as functions of the free-stream speed.

Figure 7. Frequency $f$ of the first torsional modes vs speed ${\boldsymbol{U}_\infty }$ for Goland’s wing. Coupled DUST-MBDyn simulation results (VL, NL-VL and SP mesh) and MBDyn results with 2D strip theory aerodynamic model.

Figure 8. Damping ( $\xi $ ) of the first torsional modes vs speed ${\boldsymbol{U}_\infty }$ for Goland’s wing. Coupled DUST-MBDyn simulation results (VL, NL-VL and SP mesh) and MBDyn results with 2D strip theory aerodynamic model.

Numerical results of the coupled simulations obtained using a panel mesh (SP) show slightly higher aerodynamic damping than those obtained using both VL approaches. An increase in the predicted flutter speed of approximately 3.7% is also observed. Given these minor differences in the results obtained with the two models, a vortex lattice mesh appears to be more convenient than a surface panel one, as the computational cost is reduced significantly with no significant loss in accuracy. In particular, the use of the non-linear VL elements could provide better performance for test cases characterised by more consistent viscous effects. In the considered case, as angles of attack are small and within the linear range of the aerofoil, the viscous correction is not significant. It is instead important to consider a correction for compressibility since flutter occurs at a Mach number around 0.5. In particular, for surface panels and classical vortex lattice models, Prandtl-Glauert correction was used, while for NL-VL model compressibility effect is evaluated considering 2D aerodynamic tables.

3.0 Results and discussion

The non-linear vortex lattice methods were used for the study of two rotorcraft configurations, i.e. a proprotor of the XV-15 tiltrotor operating in different flight configurations and two propellers in tandem with overlapping rotor disks. In order to discuss the advantages and limitations occurring using the NL-VL method, results obtained from simulations employing these elements are hereafter compared with experimental and high-fidelity CFD data for the first test case, while comparison with results obtained using classical elements implemented in DUST as LL and SP are presented for the second test case.

3.1 XV-15 proprotor

The proprotor of the XV-15 tiltrotor with metal blades was simulated using NL-VL elements in different flight conditions, i.e. hover condition, forward flight in helicopter mode, and aircraft mode. For the present test case, the three propeller blades were modeled as rigid blades using non-linear vortex lattice elements only. The tabulated 2D aerodynamic data of the propeller aerofoils were taken from Felker et al. [29]. DUST simulations were performed considering 10 propeller revolutions with a time discretisation of 5 ${^\circ}$ of blade azimuthal angle. In the following, DUST simulations results obtained using the NL-VL approach are compared to the recent high-fidelity numerical simulations results obtained by Jia et al. [30] using a Detached Eddy Simulation (DES) approach and to the experimental data collected in the test campaigns described by Felker et al. [29] and Betzina [31].

3.1.1 Hover flight condition

Figures 9 and 10 show the comparison of the rotor thrust coefficient ( ${C_T}$ ) as a function, respectively, of blade collective angle and of rotor torque coefficient ( ${C_Q}$ ). Moreover, Fig. 11 shows the comparison of the rotor figure of merit ( $FM$ ).

Figure 9. Comparison of the ${C_T}/\sigma $ vs collective angle for the XV-15 proprotor in hover. Experimental data from Ref. [29] (Exp), numerical data from Ref. [30] (CFD) and DUST (NL-VL).

Figure 10. Comparison of the ${C_Q}/\sigma $ vs ${C_T}/\sigma $ for the XV-15 proprotor in hover. Experimental data from Ref. [29] (Exp), numerical data from Ref. [30] (CFD) and DUST (NL-VL).

Figure 11. Comparison of the ${F_M}$ vs ${C_T}/\sigma $ for the XV-15 proprotor in hover. Experimental data from Ref. [29] (Exp), numerical data from Ref. [30] (CFD) and DUST (NL-VL).

The performance curves obtained with DUST describe quite well the behaviour of the experimental data in the whole range of blade collective angles tested. Moreover, the comparison of the curves shows that the DUST approach provides similar capabilities to a DES approach in terms of aerodynamic performance evaluation for such a case, but requires a quite lower amount of computational effort. The figure of merit comparison depicted in Fig. 11 shows a very good agreement with experiments of the overall rotor performance computed by DUST using NL-VL, particularly at low ${C_T}/\sigma $ . Moreover, the lower computational cost required by the DUST approach allows us to perform a decidedly higher amount of simulations with respect to high-fidelity CFD, thus covering with a finer step of collective angle the whole operational range described by the experimental curve. Indeed, the computational time required to complete the simulation of the rotor configuration was about 8 minutes using a workstation with a Dual Intel Xeon Gold 6230R @2.10Ghz with 104 cores processor.

A flow field representation of the present test condition is presented in Fig. 12, showing the helical vortical structure of the proprotor wake in hover computed by DUST highlighted by the iso-surfaces of Q-criterion.

Figure 12. Wake visualisation of the XV-15 proprotor in hover at ${\theta _{075}} = {15^ \circ }$ by means of iso-surfaces of Q-criterion computed by DUST coloured by Mach number.

3.1.2 Forward flight condition

In forward flight conditions, helicopter mode configurations of the XV-15 proprotor were investigated with DUST considering three shaft angle attitudes, i.e. $\alpha = - {5^ \circ }$ , $\alpha = {0^ \circ }$ and $\alpha = {5^ \circ }$ , at advance ratio 0.17. This choice enabled to investigate DUST capabilities in both propulsive and descending forward flight conditions, analogously to what was done in the work by Jia et al. [30] that was here considered together with experimental data for the assessment of NL-VL performance capabilities. Figures 13 and 14 show the comparison of the rotor torque coefficient ( ${C_Q}$ ) as a function, respectively, of rotor thrust coefficient ( ${C_T}$ ) and of lift coefficient ( ${C_L}$ ).

Figure 13. Comparison of the ${C_Q}/\sigma $ vs ${C_T}/\sigma $ for the XV-15 proprotor in forward flight. Experimental data taken from Ref. [31] (Exp), numerical data from Ref. [30] (CFD) and DUST (NL-VL).

Figure 14. Comparison of the ${C_Q}/\sigma $ vs ${C_L}/\sigma $ for the XV-15 proprotor in forward flight. Experimental data from Ref. [31] (Exp), numerical data from Ref. [30] (CFD) and DUST (NL-VL).

A quite good agreement of DUST NL-VL simulations results with experiments is found also for this flight condition. In particular, DUST results slightly underestimate the experimental performance curves in the whole range of collective blade angles considered. Nevertheless, the discrepancies of DUST results from experimental data are relatively lower with respect to CFD simulation results, particularly in the lower range of collective angles tested. On the other hand, at higher thrust levels or lift coefficient ranges, larger discrepancies with respect to experimental data are found for DUST representation of the proprotor aerodynamic performance. This behaviour could be related to the limitation of the DUST approach to reproduce accurately high blade loading conditions characterised by significantly separated flow regions. Moreover, these discrepancies could be also due to high three-dimensional flow effects occurring in this configuration.

Indeed, these flight conditions are also characterised by relevant aerodynamic interactional effects, as shown by the flow visualisation presented in Fig. 15 highlighting the interaction of the tip vortices with the downstream blades.

Figure 15. Wake visualisation of the XV-15 proprotor in advanced flight at 5 ${^\circ }$ collective, and $\alpha = {5^ \circ }$ by means of iso-surfaces of Q-criterion computed by DUST coloured by Mach number.

3.1.3 Aircraft mode flight condition

Aircraft mode configurations of the XV-15 proprotor were simulated using DUST at an advance ratio of 0.337 for different collective blade pitch angles. The capabilities of the DUST NL-VL approach to reproduce the aerodynamic performance of the proprotor in this flight condition are evaluated by comparing the torque coefficient ( ${C_Q}$ ) and propulsive efficiency ( $\eta $ ) as a function of the thrust coefficient ( ${C_T}$ ), as can be seen respectively in Figs. 16 and 17.

Figure 16. Comparison of the ${C_Q}/\sigma $ vs ${C_T}/\sigma $ for the XV-15 proprotor in aircraft mode flight. Experimental data from Ref. [31] (Exp), numerical data from Ref. [30] (CFD) and DUST (NL-VL).

Figure 17. Comparison of the propulsive efficiency eta vs ${C_T}/\sigma $ for the XV-15 proprotor in aircraft mode flight. Experimental data from Ref. [31] (Exp), numerical data from Ref. [30] (CFD) and DUST (NL-VL).

A very good agreement between DUST simulation results and experimental data is found also for this flight conditions. In particular, the discrepancies with respect to experimental curves exposed by DUST simulations results are quite lower with respect to the ones obtained by high-fidelity CFD simulations, thus confirming the suitability of DUST NL-VL approach for an accurate evaluation of propellers’ aerodynamic performance. Indeed, a slight overestimation of the propulsive efficiency evaluated by experiments can be observed from DUST simulation results over almost the whole range of rotor thrust conditions tested, while a higher discrepancy in the order of 20% is found only for the lowest blade loading condition tested.

The flow field representation of an aircraft mode flight condition computed by DUST is shown in Fig. 18 by means of iso-surfaces of Q-criterion, highlighting a quite coherent helical vortical structure of the proprotor wake without interactions due to the free-stream velocity convection characterising this flight condition.

Figure 18. Wake visualisation of the XV-15 proprotor in aircraft mode by means of iso-surfaces of Q-criterion computed by DUST coloured by Mach number.

3.2 Tandem propellers

The second rotorcraft application concerns the aerodynamic interaction between multiple propeller configurations typical of eVTOL vehicles. In particular, two propellers in tandem with different overlapping of the rotor disk were simulated by DUST using LL, SP and the NL-VL elements. DUST simulation results were compared with experimental data available from a wind tunnel campaign performed at the S. De Ponte wind tunnel of Politecnico di Milano. The experimental setup of the tandem propeller models inside the wind tunnel test section is shown in Fig. 19. Details about the propeller models and the experimental campaign can be found in the following works [32, 33].

Figure 19. Layout of the experimental setup of the tandem propellers at S. De Ponte wind tunnel of POLIMI [33].

The axial distance ( ${L_x}$ ) between the two propeller disks was set equal to five rotor radii, while several lateral separation distances ( ${L_y}$ ) between propeller rotation axes were considered. A scheme of the tandem propellers layout is depicted in Fig. 20. DUST simulations reproduced wind tunnel test conditions that consisted of runs performed with tandem co-rotating clockwise propellers both with a rotational speed controlled to 7050 RPM. This RPM target value was considered to reproduce a typical tip Mach number, i.e. ${M_t} = 0.325$ , of full-scale eVTOL aircraft propellers in cruise flight conditions.

Figure 20. Layout of the tandem propellers numerical model built for DUST simulations.

The blade geometry used to build the numerical model of the propeller was digitally created by means of a 3D scanning of the blade model. In particular, CAD software was used to generate the blade geometry from the surfaces provided by the scanning system. The maximum difference between the reconstructed blade CAD geometry and the 3D scanned surfaces was below 0.1 mm. A total number of 12 sections were extracted along the span direction. For each section, the aerofoil geometry was extracted and the distribution of twist, chord, sweep and dihedral was derived along the blade radial coordinate (r), as reported in Table 2. The blade numerical model was built considering aerofoils of the GOE, NACA and MH series reproducing the sectional geometries derived from the scan. The selected aerofoils as well as their spanwise positions along the blade are reported in Table 2. The 2D aerodynamic coefficients of the selected aerofoils were calculated by XFOIL simulations [34] in the range of angles of attack before stall. The method presented in Ref. [35] was used to calculate the post-stall behaviour of the sectional aerodynamic load coefficients in the angle-of-attack range between $ \pm 180^\circ $ .

As previously stated three different numerical models were used for DUST simulations, i.e. using LL elements, surface panels (SP), and NL-VL. For each of the three blades, both LL, NL-VL and SP models have a total of 50 elements in the spanwise direction. In particular, for the NL-VL model, a chordwise discretisation of 5 elements was used, while 30 elements were used for the SP case. The spinner-nacelle surface was modeled with 1212 surface panel elements. In particular, Fig. 21 shows the layout of the propeller mesh built for DUST simulations with blades discretised with vortex lattice elements.

All DUST simulations were performed considering a length of 10 propeller revolutions with a time discretisation of 4 ${^\circ }$ of blade azimuthal angle. The computational time required to complete the simulation of a single propellers configuration was about 10 minutes for the LL mesh and about 20 minutes for the NL-VL using a workstation with a Dual Intel^® Xeon Gold 6230R @2.10GHz processor with $52$ physical cores and $2$ threads for each core. A tandem propellers configuration simulation required about 17 minutes for the LL mesh and about 40 minutes for the NL-VL. Thus, the use of the NL-VL method does not introduce a consistent increase of computational effort, particularly if compared to the time required for the simulation of the same propeller configurations by high-fidelity CFD approach [36].

Before analysing the tandem propeller configuration, DUST simulation results obtained over simulations of the single propeller are presented. Figure 22 shows the comparison of the thrust coefficient ${C_T}$ and power coefficient ${C_P}$ as a function of advance ratio $J$ obtained with the different DUST elements compared to experimental measurements.

Table 2. Aerofoils sections, chord and twist distributions along the propeller blade span

Figure 21. Layout of the propeller model mesh with blades discretised with vortex lattice elements.

Figure 22. Single propeller performance evaluated by DUST compared to experimental data, ${M_t} = 0.325$ .

The performance curves of the single propeller obtained with LL and NL-VL elements are in quite good agreement with experimental results, while a net overestimation of both thrust and power is obtained with the blade SP elements model. This underlines the importance of viscosity effects for this propeller flight condition that are well captured by LL and NL-VL methods. In particular, the NL-VL model performs slightly better than the LL model in terms of thrust coefficient capture at high advance ratios, probably due to a better representation of the pressure distribution on blade surface thanks to chord-wise discretisation, while a very similar behaviour of the models can be found considering the evaluation of the power coefficient in the whole range of advance ratio tested.

Then, the tandem propellers test case was investigated, particularly by focusing the attention on an advance ratio $J$ equal to 0.8 representing a typical eVTOL aircraft cruise flight condition. The rear propeller performance computed by DUST for different overlapping degrees of the propellers disks is presented in Fig. 23 as a function of the lateral distance ${L_y}$ and compared with experimental data [32]. In particular, the average values of the rear propeller thrust and power coefficients were normalised with respect to the values obtained from the simulations of the single propeller configuration. This representation enables to highlight the interactional effect provided by the front propeller slipstream on the rear propeller performance.

Figure 23. Rear propeller performance evaluated by DUST compared to experimental data as a function of the lateral distance ${L_y}$ and ${L_X} = 5$ R, $J = 0.8$ , ${M_t} = 0.325$ .

The experimental data indicate that when the lateral separation distance between propellers is set to ${L_y}$ = 2R, the normalised performance coefficients tend to converge towards unity, indicating that aerodynamic interactions have a negligible effect on rear propeller performance in this configuration. As the degree of overlapping between the propeller disks is increased to ${L_y} = 1R$ , the experimental results show only minor losses in performance with respect to single propeller performance, i.e. below 5%. When the lateral distance between propellers decreases below one radius, the curves’ slope increases significantly due to the negative impact of aerodynamic interaction effects on the rear propeller’s performance. In particular, for the co-axial propeller configuration, the performance of the rear propeller decreases by nearly 30% and 20% for thrust and power coefficients, respectively. The physical explanation of the occurring phenomena causing this loss of rear propeller performance can be provided by the analysis of the local velocities computed on rear propeller blade as described in detail in Ref. [33]. In particular, this analysis showed that the ingestion of front propeller slipstream provides a remarkable increase of the axial velocity component experienced by the rear propeller blade as well as a slight variation of tangential velocity with respect to single propeller configuration. The combination of these effects are responsible of a decrease of the local effective angle-of-attack occurring on rear propeller blade and responsible of the large performance loss observed for this co-axial configuration. This effect on rear propeller performance becomes larger for co-axial configuration, as rear propeller disk is fully invested by front propeller disk, while increasing the vertical distance the effect on rear propeller performance decreases due to the lower degree or overlapping between propeller disks.

The comparison of DUST simulation results shows a quite good agreement of the LL and NL-VL computations with experimental data along the entire range of lateral separation distances tested. On the other hand, the SP model presents greater discrepancies with respect to experimental values, particularly for the configurations characterised by the higher degree of disks overlapping. This confirms the suitability of the models including viscous corrections for the performance evaluation of rotorcraft test cases with massive aerodynamic interactions between wakes. In particular, the NL-VL model performs better than the LL model for the evaluation of the rear propeller thrust loss in co-axial configuration, where the higher degree of interactional effects due to the front propeller slipstream occurs. On the other hand, LL and NL-VL models behave almost similarly concerning the evaluation of rear propeller power loss.

Further assessment of DUST capabilities to investigate the present problem can be provided by the comparison of the flow fields computed using the different numerical models built for the tandem propellers. In particular, Fig. 24 shows the comparison of the average freestream ( $u$ ) velocity components computed for the co-axial propellers configuration, representing the most demanding test case in terms of the interactional effect on performance. The numerical flow fields are here compared also to experimental data obtained by stereo PIV surveys [33].

Figure 24. Comparison of the averaged axial velocity components for tandem propeller configurations ${L_x} = 5R$ and ${L_y} = 0R$ at $J = 0.8$ and ${M_t} = 0.325$ . PIV (left) and DUST (right).

PIV results obtained for the co-axial tandem propeller configuration ( ${L_y} = 0$ ) reveal that the slipstream generated by the front propeller causes the acceleration of the rear propeller outer wake region (see Fig. 24). The acceleration provided by both propellers is accurately captured by DUST simulations with LL and NL-VL models. On the other hand, the SP numerical model provides a quite higher value of flow acceleration related to propeller blade rotation. Thus, flow field comparison further confirms the need for a numerical model including viscous corrections as the NL-VL element for an accurate investigation of complex interactional aerodynamics phenomena typical of novel rotorcraft vehicle configurations. An additional aspect to be considered is the capability of the NL-VL element to provide accurate pressure distributions along both spanwise and chordwise coordinates of the blades. This represents a promising feature that could improve the performance of aeroacoustics computation made with FWH approach with respect to the use of 1-D numerical elements as lifting lines or potential 2D elements as surface panels [36].

4.0 Conclusions

A non-linear vortex lattice element was implemented to improve the capabilities of the mid-fidelity aerodynamic code DUST for rotorcraft applications. The aim of the method is to overcome the limitations of the state-of-the-art numerical elements, such as lifting lines, surface panels and classical vortex lattice elements. Indeed this element relies on a fully implicit numerical scheme providing a higher degree of numerical stability, suitable particularly for the simulation of fluid-structural coupled problems. Moreover, these aerodynamic elements include the capability to capture viscous effects as well as to provide a discrete representation of the mean surface of lifting bodies.

The numerical method was validated and tested against different fixed-wing test cases facing both aerodynamic and aeroelastic problems. The comparison with high-fidelity numerical results obtained for a steady fixed-wing test case showed that DUST simulations performed using this approach allowed great accuracy in capturing the wing aerodynamic performance, also for nearly stalled conditions, while the numerical approach fails for conditions characterised by very large angles of attack and massive flow separations. Moreover, aeroelastic simulations performed with DUST over a wedged fixed-wing confirmed the capabilities of the approach to capture flutter speed without any loss of accuracy with respect to classical surface panels or vortex lattice methods. This result opens the scenario to the use of a stable accurate element including viscous correction for the simulation of complex aeroelastic problems.

The new numerical element implemented in DUST was then used to simulate two different rotorcraft configurations typical of tiltrotors and eVTOLs. DUST simulation results performed over the XV-15 proprotor showed a very good agreement with experimental data for different flight conditions ranging from hover to aircraft mode. In particular, numerical results obtained using the NL-VL elements to model the blades showed the capability to capture rotor aerodynamic performance with a degree of accuracy quite similar to the high-fidelity CFD approach, with a significantly lower computational cost; this further highlighted the potentialities of the implemented approach for the design and investigation of rotorcraft configurations characterised by consistent non-linear effects.

Comparison between DUST simulations results over the tandem propeller configurations typical of eVTOL in aircraft mode conditions highlights the need for numerical elements implementing viscous corrections for a fast and accurate representation of propellers’ aerodynamic performance characterised by strong interactional effects. In particular, the combined capability to provide an accurate discrete representation of pressure on blades’ mean surface opens a new scenario for the use of non-linear vortex lattices for fast and accurate evaluation of the aeroacoustics of complex rotorcraft configurations.