Parametrization Effects of the Non-Linear Unsteady Vortex Method with Vortex Particle Method for Small Rotor Aerodynamics

Abstract

The aim of this article is to investigate the parameter sensitivity of the (Non-Linear) Unsteady Vortex Lattice Method-Vortex Particle Method [(NL-)UVLM-VPM] with Particle Strength Exchange-Large Eddy Simulations (PSE-LES) method on a lower Reynolds number rotor. The previous work detailed the method, but introduced parameters whose influence were not investigated. Most importantly, the Vreman model coefficient was chosen arbitrarily and was not suitable to ensure stability for this lower Reynolds number rotor simulation. In addition, the previous work presented a consistency study where geometry and time discretization were refined simultaneously. The present article starts with a comparative literature review of potential methods used to solve the aerodynamics of an isolated hovering rotor. This review highlights the differences in modeling, discretizations, sensitivity analysis, validation cases, and the results chosen by the different studies. Then, a transparent and thorough parametric study of the method is presented alongside discussions of the observed results and their physical interpretation regarding the flow. The sensitivity analysis is performed for the three free parameters of UVLM, namely Vatistas core size, the geometry and the temporal discretizations, and then for the three additional parameters introduced by UVLM-VPM, which are the Vreman model coefficient, the particle spacing, and the conversion time. The effect of different databases in the non-linear coupling is also shown. The method is shown to be consistent with both geometry and temporal refinements. It is also consistent with the expected behavior of the different parameters change, including the numerical stability that depends on the strength of the LES diffusion controlled by the Vreman model coefficient. The effect of discretization refinement presented here not only shows the integrated coefficients where different errors can cancel each other, but also looks at their convergence and where relevant, the distributed loads and tip singularity position. Finally, the aerodynamics results of the method are compared for different databases and with higher fidelity Unsteady Reynolds Averaged Navier–Stokes (URANS) 3D results on a lower Reynolds number rotor.

1. Introduction

The improvement of next-generation rotorcraft performances is related to the development of higher quality preliminary design tools. To better estimate the rotor aerodynamics, those new design tools must reduce the number of assumptions limiting the predictive capabilities of the physical principles [1]. To be usable in early design phases, the software need to be fast, accurate, and reliable. The faster software currently available are based on the low fidelity Blade Element Momentum Theory (BEMT) method, which fails to capture important features such as wake unsteadiness and aerodynamic interactions. The more accurate high fidelity methods, such as the URANS simulations, yield very good aerodynamic prediction, but are far too expansive to be used in early design. That is why researchers are aiming to develop medium fidelity tools based on tridimensional potential methods. Those methods need only to mesh the geometry, not the surrounding air like the higher fidelity methods, meaning a much simpler mesh generation and a significantly reduced computational cost. Potential methods can produce very good results when their hypotheses are respected, especially for rotary wing aerodynamics because accurate wake development is key in that context. They avoid the artificial dissipation issues requiring higher fidelity method to have a densely refined mesh in the areas of interest in the wake to avoid over-diffusion of vorticity.

Potential methods are appropriate to model helicopter rotor blades at low to medium rotational speed. However, since they are built on the potential aerodynamics assumptions (incompressible and inviscid flow), they should not be used directly for high speed transonic flows or smaller rotors such as those of ever more popular small Unnamed Aircraft Vehicles (UAVs) rotors. Low computational cost coupling algorithm were developed to incorporate locally any desired aerodynamic effect, such as viscosity or compressibility to potential methods [2,3]. This is appealing for the rotorcraft aerodynamics because the rotor wake is subsonic, even for an helicopter rotor at high advance ratio that experiences transonic speed on a small area of the advancing blade. Good hovering aerodynamic prediction were obtained on a small (low Reynolds number) rotor both for the thrust ($C_T$) and the torque ($C_Q$) coefficients [4]. The present method, as presented in ref. [5], is similar to that and has achieved very good aerodynamic prediction on a moderate Reynolds and Mach number rotor.

No matter how good some results have been, more work is still needed on potential methods. One of the issues with potential method is that they often rely on tweaking parameters, such as core size or dissipation and destruction thresholds, whose effects are not always clearly shown by researchers. It then becomes hard to reproduce the results without having the values of the said parameters, or understand how that parameter might affect the simulation. Another issue is the consistency assessment. The objective of this article is to propose improvements to the method presented in ref. [5], named NL-UVLM-VPM, by exposing the effects of its many parameters and performing mesh refinement studies which are lacking in the literature.

The three improvements to the method presented in this article are the computation of the overlapping factor, the possibility to keep wake panel rows before converting the panels in particles, and the ability to convert a panel side in multiple particles to better approximate long vortex lines and better match different time steps particle spacing at the tip while retaining root stability. Improvements 1 and 3 now grant spatial and temporal consistency separately and improvements 2 and 3 help to increase the time step for NL-UVLM-VPM simulations up to 20° compared with the $2.5°$ used in the previous work. The increased time step allows to test the NL-UVLM-VPM stability for far more rotor revolutions than in the previous work.

The seven parameters of the present sensitivity study are listed below with the reason why their sensitivity was not fully analyzed in the previous work [5]:

1.

UVLM Vatistas core size ($\sigma$): It was not needed in the previous work, because the wake panels were instantly converted into particles. In the present work, some rows of wake panels are kept behind the blade to help reduce the near-field discretization error of the velocity induced on the blades because the discrete particles approximation of the straight-line vortex elements are now farther from the rotor plane. Reducing the impact of particles on the blade allows coarser particles discretization than before. However, the prescribed Vatistas core size provided to smooth the wake panels induced velocity singularity needs to be carefully selected to achieve stability without significantly affecting the results;

2.

Geometry discretization refinement (mesh): In the previous work, the geometry discretization refinement was conducted simultaneously with the time discretization refinement. The method appeared to be consistent with refinement. After the publication, it was realized that independently varying the geometry and time discretizations was not consistent. The reasons for this inconsistency and the solution to the problem are detailed at the beginning of Section 3.2;

3.

Time discretization refinement (time step, $\Delta \Psi$):

$$\Delta \Psi =\Delta t·\Omega$$

(1)

where $\Delta \Psi$ is the azimuthal increment per time step in degree, $\Delta t$ is the time step in seconds and $\Omega$ is the rotational speed in degree per second. It was previously performed jointly with the geometry discretization refinement;

4.

Vreman model coefficient ($C_v$): In the previous work, the Vreman model coefficient was set to 0.07 as it is the equivalent value to the theoretical value of the Smagorinsky constant for homogeneous isotropic turbulence. The theoretical value kept all the simulations stable, so it was thought to be appropriate. However, the theoretical value did not keep the VPM stability on the smaller radius and aspect ratio rotor of the current work. Vreman states that to obtain robust simulations for complex practical cases, the value can be higher or lower than the theoretical value [6]. In this work, it is observed that the value needed for this constant varies with the tip particle spacing;

5.

Tip particle spacing ($\Delta x_{tip}$):

$$\Delta x_{tip}=\frac{\Delta \Psi }{N_{tip}}$$

(2)

where $\Delta x_{tip}$ is the tip particle spacing and $N_{tip}$ is the number of streamwise particles per time step at the tip of the blade. In the previous work, a straight-line vortex line was always converted in a single particle, therefore it was not possible to test that parameter independently from the time step. The time step was much smaller, so it was a reasonable simplification. With the larger time step used in this work, it is important to have multiple particles on the longer straight-line elements to avoid discretizing them with a single enormous discrete particle;

6.

Wake-particle conversion time (revolution): In the previous work, the wake panels were instantly converted into particles, so that parameter was not present. In the current work, keeping the wake panels helps to increase the time step by moving farther from the rotor plane the discretization error caused by the discrete particles;

7.

Database for the non-linear coupling: Two different databases were tested along with the linear UVLM in the previous work. In this work, since the parametric study focuses on a single case, four different databases are tested in the Results section of the present article.

The numerical challenges to obtain accurate rotor blades aerodynamic prediction are tightly linked with the flight conditions that are simulated. For instance, if the rotor has a high advance ratio or climbing speed, the wake is convected away from the rotor plane and becomes less important than if it is in hover or descending flight. However, the high advance ratio blade might be subject to transonic effects on the advancing side while experiencing low speed or even reverse flow near the root on the retreating side [7]. Thus, accurately capturing the varying speed aerodynamic loads is more important to the advancing flight condition than the long distance wake development. Inversely, the steady hovering rotor blade is greatly impacted by the slowly developing wake under the rotor and not normally subject to highly unsteady loads. In the mid 2010s, the American Institute of Aeronautics and Astronautics (AIAA) prediction workshop focused its attention on the hover condition because it represents the true value of the helicopter and it is a limiting design point in terms of power requirement [8]. This work also limits its scope to hovering rotors, but there are no hypotheses contained in the present method that should prevent it to be extended to other flight conditions. The current sensitivity study focuses on the hover condition, because accurate wake capturing is crucial in this condition and most of the parameters of the method directly affect the wake development. In addition, we chose an isolated rotor to avoid unnecessary aerodynamic interactions that could theoretically be captured by the present method, but would only add layers of complications in the results interpretation. A small Reynolds number is chosen as it is more challenging.

This article starts, in Section 2, with a comparative literature review of similar methods and their isolated hovering rotor validation cases that shows gaps in our understanding and the lack of a common simplified validation cases for smaller rotor aerodynamics. Section 3 presents a summary of the previous method and the proposed improvements. Section 4 describes a simplified validation case and the stripwise dataset generation. Section 5.1 exposes the parametrization study of the aerodynamic solver (UVLM) and Section 5.2 the parametrization study of the combined UVLM and Vortex Particle Method (VPM) approach. Section 6 presents results for the high-fidelity URANS approach and compares all methods. Finally, the article is concluded in Section 7.

As this paper requires extensive understanding of several scientific concepts including a comprehensive literature review, summaries are presented at the end of the three largest sections to highlight the key takeaways while providing the option to skip the complete descriptions.

2. Literature Comparison

The present literature review attempts to compare some of the most recent literature of potential methods for solving the aerodynamics of an isolated rotor in hover. As it will be seen, direct comparison is hard to achieve since differences are found not only in the methods themselves, but also in their usage and validation. Section 2.1 will compare the different methods and their discretization. Since most of these studies are validated on a different validation case, the validation cases are then presented in Section 2.2. After that, the comparison of the results produced by the different methods presented in Section 2.1 with their validation case detailed in Section 2.2 are shown in Section 2.3. To better appreciate the differences between the different work, tallies have been constructed and are presented in those three sections. Care has been taken to present the data as accurately as possible and to compare the methods on a similar footing, but some information was not directly available and had to be computed, if possible, from the data at hand. The reader is invited to consult the original references. Finally, a conclusion of the literature comparison is provided in Section 2.4.

2.1. Similar Methods Summary

Table 1. Summary of the similar methods and their discretization for solving isolated hovering rotors aerodynamics from the literature. Information shown with a question mark (?) is not clearly stated in the referenced work and is inferred from context, as detailed in the text for each occurrence.

Authors	Method	Compressibility	Viscosity	Mesh	$\mathsf{\Delta }\mathsf{\Psi }$ (°)	Sensitivity	Convergence (Rev)
Colmenares et al. [9]	UVLM	No	No	6 × 15 cos	10?	No	12
Perez et al. [10]				7 × 33	10	$C_T$	9
Alvarez and Ning [11]	VPM	Prandtl–Glauert	Lookup	1 × 50	5	No	20
Yucekayali [12]		No		1 × 30	7.5	No	6
Zhao and He [13]	LL-VPM	No need with look-up	Lookup	1 × 50	3.75	Downwash at fixed $C_T$	No
Singh and Friedmann [14]	UVLM- VPM	Karman–Tsien	Sectional drag	2 × 8	12	Not shown	6
Ferlisi [15]		No	No	4 × 15? cos	10	No	10
Tan and Wang [16]	UPM-VPM	No	No	60 ×  20	5	No on hover	10
Tugnoli et al. [17]	NL-LL- VPM	No need with non-linear	Alpha or Gamma	1 × 16	9	Not shown	20
Samad et al. [18]	NL-UVLM	Prandtl–Glauert	Alpha	10 × 25	15	No	24
Jo et al. [4]	NL- UVLM- VPM	No need with non-linear	Gamma	15 × 30	5	$C_T$ and $C_Q$	30
Lee et al. [19]				20 × 45? cos-cos	5?	No	20?
Previous work [5]			Alpha	8 × 20 cos	2.5	$C_T$ and FM	24
Current work			Alpha	16 × 40 cos	10	$C_T$, FM, $C_t$, $C_q$ and tip vortex	28

presents a summary of similar methods and their discretization for solving isolated hovering rotors aerodynamics. The information found in this table is the authors’ names and the work reference, the method used in that work, whether the method is corrected for compressibility and viscous effects, the blade mesh discretization (number of chordwise elements x number of spanwise elements), the time discretization (the blade azimuthal increment per time step $\Delta \Psi$, as defined in Equation (1)), if that work attempted any sensitivity study such as a refinement study or the sensitivity of some parameters, and the number of revolution simulated to obtain convergence. The following paragraphs will explain these columns in more details.

The method column shown in Table 1 refers to the methods used to model the lifting blades and their wake. The two can be the same, but do not need to be in general, explaining why certain studies have a single method and others have two. However, most potential rotor blade aerodynamic solvers are using Vortex Particle Method (VPM) wake representation, which allows the particles to be convected freely in separate directions without causing unrealistic filament or panel stretching. The VPM can also be used to model the blade bound circulation at discrete locations [11,12]. The lifting line (LL) blade representation replaces the particles by straight line vortex elements ensuring a definition of the circulation over the whole span of the blade instead of the discrete particles location. The UVLM is the surface extension of the lifting line, allowing circulation to vary chordwise as well as spanwise. The Unsteady Panel Method (UPM) is composed of source and doublet panels that allow the thickness of the blade to be modeled by having panels on the upper and lower sides of the blades instead of the thin camber surface used by the UVLM. The effect of thickness, camber, as well as non-linear effects not captured by potential flow hypotheses such as compressibility and viscosity can be incorporated at small computational cost in the 3D potential methods via 2D sectional aerodynamic coefficients. To do so, simple lookup method or non-linear (NL) coupling algorithm can be used. The non-linear coupling can be performed using the circulation ($\Gamma$) [20] or the angle of attack ($\alpha$) [2] coupling variable. NL algorithms differ from simple lookup with the subiterations that effectively captures the difference in downwash caused by the change of lift of the sections within a time step. Since this article focuses on isolated rotors in hover reaching a steady state, upon convergence, the lookup methods should converge to the same equilibrium as the non-linear methods, since the time iterations play a similar role as the non-linear subiterations. However, the lookup methods have a lag in capturing the change in downwash, so they might be less appropriate in unsteady loading. There are not many reasons for not doing the alpha coupling algorithm, aside from simplicity, since the computational cost is negligible compared with wake development. Compared to the Gamma coupling, the alpha coupling is unambiguous even in post-stall and requires significantly less under-relaxation. The non linear implementations of Jo et al. [4] and Lee et al. [19] are using a Gamma coupling and an arbitrary point along the section chord to compute the induced velocities to obtain the effective angle of attack. The problem with that is that the singularities have an enormous effect on their near field, so if at a given time step a singularity is close to one of the control points, the effective angle of attack will be drastically affected.

By construction, the potential methods without viscous database do not capture the compressibility nor viscous effects. There exists corrections based on the Mach number that can be used directly in the method as done by Singh and Friedmann [14] using Karman–Tsien correction to account for compressibility. The viscous effects cannot be modeled by such a simple correction. That work introduces a viscous drag component in the method via power series approximation assuming a non stalled high dynamic pressure flow and knowledge of different profile aerodynamic coefficients. Unfortunately, that does not take into account viscous effects on the lift of the sections. As stated previously, it is possible to simultaneously capture the compressibility and viscous effects on the lift and drag of the blade sections via lookup tables, with or without the non-linear coupling algorithm. The tabulated 2D data might be obtained from higher fidelity computation or from experimental data at the desired condition (Mach and Reynolds) as done in the alpha coupled NL-UVLM-VPM of the previous work [5] with URANS and experimental correlation for the tabulated data. Alvarez and Ning [11] still use a compressibility correction despite the use of tabulated data, because their database is computed with the incompressible 2D Xfoil software. Jo et al. [4] use polars computed with Standford University Unstructured (SU2) flow solver with the Spalart–Allmaras turbulence model and Bas–Cakmakcioglu transition model. They could have used incompressible database since their validation case is subsonic, but there is no harm in using a compressible software for subsonic aerodynamic. The work of Samad et al. [18] uses a Prandtl–Glauert compressibility correction to the UVLM as well as compressible database generated with RANS software SU2 and Spalart–Allmaras transition model. One has to be careful in doing so, because using a compressibility correction in the UVLM changes the $C_{l\alpha }=2\pi$ lift slope of a UVLM section. It is simpler to have the compressibility applied only to the database. It is unfortunate that some studies (Yucekayali [12], Zhao and He [13], and Lee et al. [19]) do not detail the origin of their tabulated data, because as will be seen in the present work, the database has a strong impact on the final results. Tugnoli et al. [17] use RANS database, and it seems like it is with the Spalart–Allmaras transition models computed with ROSITA, however it is not clear in the article if that clarification applies only to the 3D high fidelity validation.

Regarding the mesh and time discretization, some of the information in Table 1 is presented with question marks (?) because the information is not clearly stated in the referenced work. For instance, Colmenares et al. [9] give the mesh discretization for the isolated rotor in hover, but only give the time discretization of $\Delta \Psi =10°$ for the second validation case (the two-rotor aircraft in hover). The best guess is then to assume the same time discretization for the isolated rotor. On the contrary, Ferlisi [15] writes the time discretization, but not the blade discretization. However, looking at previous figures in the thesis (Figure 1.1 and Figure 2.6), we can count 4 × 15 panels with what looks like a cosine tip refinement. The mesh and time discretization used by Lee et al. [19] is not stated, but a previous work of the authors states a 20 × 40 mesh with cosine refinements at both ends of the blade and $\Delta \Psi =5°$. They most likely used a similar discretization. The studies modeling the lifting blades with VPM or LL are shown with 1 chordwise element as they do not capture chordwise variation. The work of Tan and Wang [16] is using significantly more chordwise elements compared with the other methods, partly because the UPM requires elements on the upper and lower sides of the profile. We can then see that most of the studies are using between 4 and 20 chordwise elements and between 15 and 50 spanwise elements and between $\Delta \Psi =5-10°$ per step. The work of Singh and Friedmann [14] is the one using the coarsest blade discretization with only 2 chordwise elements by 8 spanwise elements as well as one of the coarsest time discretization with $\Delta \Psi =12°$ per step, with only the work of Samad et al. [18] using a coarser time discretization with $\Delta \Psi =15°$ that is more appropriate with the continuous wake trailed vortices of the UVLM compared with the discrete particle wake of the VPM.

The work of Singh and Friedmann [14] do not show any refinement or sensitivity studies to justify their discretization. They state that reducing the time step do not change the results without giving further details and without addressing a change in mesh discretization. As can be seen in the sensitivity column, most of the studies fail to present a satisfactory sensitivity analysis. Indeed, 8 of the 13 studies did not show any sensitivity analysis at all. The work of Perez et al. [10] only assess the effect of mesh and time refinement on the global $C_T$ coefficient, without presenting the effect of the refinement on $C_Q$ or on the spanload. It is unfortunate since their final results do show $C_Q$. The work of Zhao and He [13] presents a sensitivity analysis of the wake cutoff, the time-step, the overlapping parameter, and the resolution parameter on the rotor-induced downwash. That analysis does not include the blade mesh discretization. Results show some sensitivity to those parameters, except for the overlapping parameter, but the authors state that the difference is not significant without quantifying it. The choice of downwash makes it is hard to appreciate how the aerodynamic loads are affected by the parameter sensitivity, especially since the simulations were run at constant $C_T$, not constant collective. That means that if a parameter affects the thrust generated by the blade, the collective angle is tuned to match the total force. Then, the blade loading is mostly unaffected, unless the loading shape is significantly different while integrating to the same amount. With the same blade loading, the particles shed into the wake have a very similar strength regardless of the parameter effect. Consequently, the downwash in close proximity of the rotor is expected to be very similar. Thus, the comparison of the downwash sensitivity at constant $C_T$ does not show the full extent of the simulation sensitivity to the parameters. Tan and Wang [16] performed a sensitivity study of the wake cutoff distance, the time step size, the spanwise number of elements and smoothing parameter, but unfortunately chose to perform such an analysis on the forward flight validation case, where the rotor is moving away from the wake, reducing the long time effect of the wake on the blade aerodynamics. In Tugnoli et al. [17], it is said that the parameters represent a choice based on the dependence study performed in a previous work. However, the cited reference does not show a refinement or parametric study. The work of Jo et al. [4] presents a detailed sensitivity analysis, showing the insensitivity of the blade and time discretization on the global $C_T$ and $C_Q$ coefficients (showing that 2 of the 3 global coefficients is sufficient to recover the third if needed). They also present the mean and RMS pressure loading on the blade surface for the different blade mesh and time discretizations. The distributed loads are not considered and are most likely still affected by grid refinement when looking at the pressure loading variation between the refinement levels, despite the overall coefficients convergence and converged pressure loading for the temporal refinement. Indeed, the zone of maximum mean pressure near the leading edge tip appears to be about 50% longer in the spanwise direction on the finest grid compared to the medium grid for about the same width in the chordwise direction. The previous work [5] also compares the effect of mesh and time discretization on global $C_T$ and FM refining both mesh and time at the same time instead of perturbing mesh and time in separate analyses. The effect of refinement on spanload is not shown either. The current work shows the effect of all the parameters of the NL-UVLM-VPM with PSE-LES on the global and distributed loads and, where relevant, on the tip vortex position.

Finally, the last column of Table 1 shows how many rotations were simulated to obtain steady state results. The information is available for all of the studies, except for Zhao and He [13] and Lee et al. [19]. For the latter, the value given in a previous work is shown with a question mark. Except for the studies of Yucekayali [12] and Singh and Friedmann [14] that consider their results converged after 6 revolutions, all of the work have simulated between 9 and 30 revolutions. The convergence in time can also be affected by the number of slow started revolutions, where Colmenares et al. [9], Perez et al. [10], Samad et al. [18], the previous [5], and current studies respectively used 1, 1.11, 2, 3, and 2 slow started revolutions. Lee et al. [19] does not mention slow started revolutions, but used 1 slow started revolution in a previous work. The other studies do not mention slow started revolutions.

2.2. Similar Methods Validation Cases

Table 2. Summary of the validations used by the similar methods for solving isolated hovering rotors aerodynamics from the literature.

Validation	Exp. Year	$\mathit{AR}$	${\mathit{N}}_{\mathit{b}}$	RPM	Mach Tip	Reynolds
Harrington [22]	1951	20.4	2	477	0.44	$R{e}_{tip}\approx R{e}_{75}=1.2$ M
		6.67	2	250	0.23	$R{e}_{tip}=2.3$ M
				374	0.35	$R{e}_{tip}=3.5$ M
Boatwright [23]	1972	16.2	2	245	0.40	$R{e}_{tip}=2.2$ M
				340	0.56	$R{e}_{tip}=3.1$ M
Caradonna–Tung [21]	1981	6.00	2	1250	0.44	$R{e}_{tip}=1.8$ M
				1500	0.52	$R{e}_{tip}=2.2$ M
				1750	0.61	$R{e}_{tip}=2.5$ M
Hariharan et al. [8]	1985 (2016)	14.6	4	1484	0.60	$R{e}_{tip}=1.1$ M
Shinoda [24]	1996	14.6	4	293	0.60	$R{e}_{tip}=5.6$ M
Droandi et al. [25]	2015	5.50	4	1120	0.32	$R{e}_{tip}=987$ k
Zawodny et al. [26]	2016	6.67	2	5400	0.20	$R{e}_{tip}=70$ k
Zhou et al. [27]	2017	5.73	2	4860	0.18	$R{e}_{70}=62$ k
Perez et al. [10]	2019	6.59	4	4500	0.25	$R{e}_{75}=110$ k
Present work	NA	8.00	2	1000	0.15	$R{e}_{tip}=170$ k

summarizes all the validation cases used by the numerical studies presented in Table 1. There are fewer validation cases than numerical studies because many numerical studies used the popular Caradonna–Tung [21] experiments as their benchmark case. All of the other validation cases are used by one numerical work in Table 1. Before detailing which validation is used by which numerical work in

Table 3. Summary of the similar methods results with respect to their validation test case when solving the aerodynamics of an isolated hovering rotor from the literature. * All the numerical studies validated with Caradonna–Tung [21] only used RPM = 1250, except for Tan and Wang [16] that used RPM = 1250, 1500, and 1750 at constant collective. ^† Errors computed with the $C_Q$ vs. $C_T$ trend line.

Authors	Validation	Validation Points	|$\mathsf{\Delta }{\mathit{C}}_{\mathit{T}}$| (%)	|$\mathsf{\Delta }{\mathit{C}}_{\mathit{Q}}$| (%)	Trimmed ${\mathit{C}}_{\mathit{T}}$	Span Load	Tip Vortex	Convergence Figure
Colmenares et al. [9]	Caradonna–Tung [21]	3	6	No	No	$C_t$	No	$C_T$
Perez et al. [10]	Perez et al. [10]	1	CFD: 8 Exp: 12	CFD: 18 Exp: 5	No	$C_p$ contour	No	$C_T$ and $C_Q$
Alvarez and Ning [11]	Zhou et al. [27]	1	0.6	No on hover	No	No	No	$C_T$
Yucekayali [12]	Shinoda [24] and Caradonna–Tung [21]	5 with [24] 1 with [21]	No	7 † with [24]	Yes	$C_l$ with [21]	Yes	$C_T$
Zhao and He [13]	Boatwright [23] and Caradonna–Tung [21]	1 with each ref	No	No	Yes	$C_l$ with [21]	No	No
Singh and Friedmann [14]	Harrington [22]	22	No	[4–12] †	Yes	No	No	No
Ferlisi [15]	Caradonna–Tung [21]	3	4	No	No	No	Yes	$C_T$
Tan and Wang [16]	Caradonna–Tung [21] *	3	No	No	No	$C_p$ sections	No	No
Tugnoli et al. [17]	Droandi et al. [25]	11	No	Exp: 9 †	Yes	No	No	$C_T$
Samad et al. [18]	Caradonna–Tung [21]	3	6	No	No	$C_l$	No	$C_T$
Jo et al. [4]	Zawodny et al. [26]	1	CFD: 3.6 Exp: 0.4	No	No	No	No	No
Lee et al. [19]	Caradonna–Tung [21]	3	Small	No	No	$C_t$	Yes	No
Previous work [5]	Hariharan et al. [8]	6	CFD: 3 Exp: 8	CFD: 2 † Exp: 5 †	Yes	$C_t$ and $C_q$ $C_p$ sections	No	$C_T$ and FM
Present work	Present work	1	[0.3–5.3]	[0.3–8.4]	No	$C_t$, $C_q$ and $C_m$	No	$C_T$ and FM

, the different validation sources are briefly explained.

All of the validation cases presented in Table 2 are experimental results, except for the present work and the Hariharan et al. [8], which is a compilation of the results from different higher fidelity 3D URANS software gathered during the 2nd Hover Prediction Workshop compared with the experimental results of Balch [28,29]. The information found in Table 2 is the name of the authors of the validation case and the reference, the year the experiment results were published, the estimated Aspect Ratio ($AR$) of the blades, the rotor number of blades ($N_b$), the rotational speed of the blades in Rotations Per Minutes (RPM), the Mach number at the tip of the blades, and the Reynolds number at a given location (generally at the tip of the blade or at 75% of the span). At first glance, the four last validation test cases in that table are different than the previous benchmark cases with a significantly lower Reynolds number and Mach number well below compressibility effects. That is because of the increased interest in small UAVs rotors, compared with helicopter-like rotor models or full scale rotors. The benchmark case of Droandi et al. [25] is also a little bit different in nature as the rotor is optimized for tiltrotor aircraft. Reynolds and Mach were computed as truthfully as possible when not provided in the references, assuming standard air when temperature and pressure are unknown. The $AR$ is easily calculated for rectangular blades, but is more complicated for complex geometry blades such as the commercial small UAVs rotors and is shown here as estimated with the available information about the geometries. The $AR$ was calculated as:

$$AR=\frac{R_{max}-R_{min}}{c_{75}}$$

(3)

It is shown to give the reader an idea of the relative slenderness of the blades, but the reader is warned that in many cases, it is an estimate, not a precise value.

The Harrington [22] is the oldest of the validation cases presented in Table 2 by a significant margin. The experiments measured the steady state thrust ($C_T$) and torque ($C_Q$) of rotors in isolated and coaxial configurations. The experiment does not report the collective angle of the blades, so that validation case does not allow to validate a software ability to correctly predict the thrust at a given collective. Only the rotor efficiency can be compared. The experiment does not report spanload or wake measurements. The rotor 1 blades were made of wood and had a tapered chord and thickness, while the rotor 2 blades were made of metal and only had a tapered thickness (rectangular platform). Both rotors blades used NACA four-digit symmetrical airfoil sections and were untwisted. The isolated rotor 1 results are reported at a blade tip speed of 500 ft/s (Mach = 0.44), resulting in a Rey7nolds number of 1.2 M at 75% of the span as well as at the blade tip. The larger inertia rotor 2 results are given at tip speeds of 262 ft/s and 392 ft/s (Mach = 0.23 and 0.35), resulting in tip Reynolds numbers of 2.3 M and 3.5 M, respectively. This means that the rotor 1 blade sections are strictly subsonic up to 68% of the blade and not too drastically affected by compressibility effects, even at the tip of the blades, while the vast majority of rotor 2 blades are subsonic. The Reynolds numbers are moderate for rotor 1 and close to that of full scale helicopter blades for rotor 2. Therefore, the hypothesis of incompressible and inviscid flows of the potential methods are not much challenged on that validation case.

The Boatwright [23] experiment was mainly focused on obtaining the 3D wake velocity components under a full scale OH-13E helicopter rotor. Three different combinations of blade pitch and rotor speed resulting in 2 different $C_T$ were tested. Since this experiment is mostly focused on wake measurements, there is little information about the actual loading of the blades. $C_T$ vs. collective is given at the three test conditions and a rotor efficiency curve ($C_T$ vs. $C_Q$) is shown with more data points. The blade geometry is not properly described in the article and the OH-13E is an old helicopter whose blade geometry details are harder to find. However, assuming no root cutout and linear chord variation, the reported blade area is retrieved and $AR$ can be estimated. The sole work in Table 1 using this experiment as a validation case (Zhao and He [13]) only compares wake velocities at the lowest Mach condition and do not show their method capacity to predict $C_T$ vs. collective or $C_T$ vs. $C_Q$.

The Caradonna–Tung [21] experiment results are still often used to validate hover prediction capabilities of potential software. Indeed, nearly half of the numerical studies presented in Table 1 present some sort of comparison with those results. The experiments gathered steady state thrust, pressure coefficient at 60 points along 5 radial sections as well as hot wire tip vortex position at different blade collective pitch and rotational velocities. However, the torque of the rotor was not measured, meaning that this validation case cannot be used for rotor efficiency prediction. The two-bladed rotor had untapered and untwisted NACA0012 profile blades. This experiment is the only one in Table 2 using rectangular untwisted and untapered rotor blades. The geometry is therefore more academic than representative of modern helicopters or UAVs rotor blades, which better spreads the thrust over the span in order to smooth the thrust drop-off near the blade tip. Consequently, the tip vortex is very strong compared with the inner vortex sheet, resulting in large interactions between the rotor and the wake and ultimately large induced drag. It is therefore a challenging configuration where accurate prediction of the wake development is key to an accurate aerodynamic prediction. Despite the many rotational velocities reported, most of the potential software validation is done using only the lowest rotational velocity ($M_{tip}=0.44$) for which the coefficient of thrust is given, resulting in only 3 collective tests results to compare. In doing so, the hypothesis of incompressible and inviscid flows of the potential methods are still very reasonable. Only Tan and Wang [16] among those using the Caradonna–Tung validation found in Table 1 use data at RPM = 1500 and 1750, but at constant collective angle. None are validated with the data at RPM = 2268 ($M_{tip}=0.79$). The simplicity of the geometry combined with the strength of the potential methods to capture accurately wake evolution without dissipation are most likely the main reasons why this validation case became and appears to remain popular.

For the AIAA 2nd Hover prediction workshop held in 2015, several existing benchmark cases were considered. The scaled 1/4.71 S-76 rotor test of Balch [28,29] was chosen because of the public availability of the relatively modern 4-bladed rotor aerodynamic design and the effect of different tip shapes. The experiments of Balch reports collective, $C_T$ and FM for different tip Mach numbers and different tip shapes in and out of ground effect. However, the spanload and the tip vortex position were not measured as part of the experiment. The compilation of results of Hariharan et al. [8] helps in this regard, because it compares not only the global coefficients, but also the thrust and torque spanwise distributions, pressure coefficient at 3 radial sections along the blade and the tip vortex position predicted with the state-of-the-art software from across industry, government agencies, and academia. Of the 7 participants, most methods were based on the RANS equations, where 2 participants had hybrid RANS-Lagrangian solvers and 1 solver was based on the Detached Eddy Simulation (DES) equations. In our opinion, this is the most thorough validation case presented in Table 2. It is not surprising that many potential solvers decides to avoid this validation case as the Mach is higher and the Reynolds is lower than the previous validation cases, both effects deviating from potential assumptions. The previous work [5] clearly showed the improvement of using non-linear coupling on this test case with potential methods.

The experimental work of Shinoda [24] is similar to that of Balch [28,29] because it uses the same geometry, but this time at full scale meaning a much larger Reynolds number. It has the same tip Mach number, well in the compressible aerodynamic range. This experiment is mostly focused on forward flight and acoustic and is also doing hover tests. It also reports collective, $C_T$, and FM, but this time only out of ground effect at one tip Mach number and tip geometry. The main difference with the Balch experiment remains not having the numerical validation provided by the Hover prediction workshop to fill the gaps concerning spanload and tip position.

The validation benchmark presented by Droandi et al. [25] is interested in the aerodynamic interactions between wing and rotors for tiltrotor aircraft. The isolated rotor had collective, $C_T$, and $C_Q$ global coefficients measured experimentally as well as the wake velocities under the rotor obtained with Particle Image Velocimetry (PIV), thereby providing tip vortex position. The article does not provide a $C_T$ vs. collective figure or table. Instead, the collective of the limiting case and the PIV configuration are reported in the text. The spanload was not measured as part of the experiment. The multi-objective tilt rotor optimization process produced a complex blade geometry with 7 profile changes as well as non linear chord and twist distributions over the span. Backward sweep was added at the tip after the optimization to reduce compressibility effects in propeller mode and forward sweep was introduced inward to keep the aerodynamic center at the feathering axis. This experiment has a Reynolds number slightly lower than the previous experiments, but a Mach number low enough to ensure an incompressible airflow.

The work of Zawodny et al. [26] is focused on measuring the hover thrust and acoustics of two isolated rotor representative of small-scale rotary-wing unmanned aircraft systems (UAS). The aerodynamic measurements are limited to global thrust. The experiments are conducted at different rotational speeds on the two 2-bladed rotors. As it is typically the case for UAVs rotors, the two rotors have fixed-pitch blades. The first rotor geometry is a replica of the DJI 9443 and the second rotor is the APC 11X4.7SF. Despite the experiments being conducted at different rotational speeds, only $R{e}_{tip}=70$ k is shown in Table 2 since it is the only validation point considered by the work of Jo et al. [4] using this experiment as validation. The much lower Reynolds number as compared with previous validation cases means that the flow around the blades fall into the laminar-transitional regime, so it should not be treated with potential methods inviscid assumption.

The experimental case presented by Zhou et al. [27] investigates the effects of rotor-to-rotor interactions on the aerodynamic and aeroacoustic performances of a small UAVs rotors. Thrust measurements, PIV, and Stereoscopic PIV were performed. The rotor torque and the blade spanload were not measured. Only a single measurement point is reported for the isolated rotor. The rotor blade geometry is composed of a single profile, with a chord varying with the optimal chord length equation and a pitch down twist from 20° to 10°. Unfortunately, the twist distribution is not specified. The rotor Reynolds number is slightly higher than the one of the work of Zawodny et al. [26] if it were transferred from 70% of the span to the tip station. Reynolds is still low enough to ensure laminar-transitional regime over the blades.

The work of Perez et al. [10] is a joint experimental-numerical analysis. They performed their own flight test using a commercial quadcopter while also performing the potential UVLM and the higher fidelity CFD simulations. The chosen quadcopter was the ARAKNOS V2 made by ADVECTOR. The atmospheric pressure and the average temperature were measured with a meteorological station. The RPM-PWM correlation for the motors was measured at the ground station with strobe light tachometer Monarch. In addition, a voltage/current sensor was installed on the quadcopter. Keeping the quadcopter in hover and knowing its weight, $C_T$ was computed. The torque of the motor was not measured, but estimated using the current, angular velocity, and the motor efficiency. No measurement of the blade spanload or the wake were performed. For the quadcopter to remain in hover, knowing that the blades collective cannot be changed, only one rotational speed could be tested. Obviously, since this is a quadcopter, aerodynamic interactions between rotors and also with the fuselage are present in this experiment. Of the three small UAVs rotors validation cases, this is the only one that reports torque despite not being a direct measurement. Potential methods therefore have few numerical-experimental benchmark to assess their ability to accurately predict the torque on these lower Reynolds, smaller aspect ratio rotors. Once again, Reynolds is low enough to need some viscous treatment in potential methods to obtain good results.

For the sensitivity study of the present work, a lower Reynolds rotor with a simpler geometry was desired. Therefore, a new test case is defined in Section 4. Since there are no experimental results on this simpler test case, the present method is compared with higher fidelity in Section 6.2.

2.3. Similar Methods Validation Comparison

Comparing the same numerical studies as in Table 1 and Table 3 presents the studies validation and a summary of the results. The columns show the name and reference of the numerical studies, the name and reference of the validation case, the number of validation points simulated with their numerical method, the absolute value of the relative error for the global $C_T$ (|$\Delta C_T$| (%)) and $C_Q$ (|$\Delta C_Q$| (%)), whether the results are shown at the same $C_T$ as the validation, what type of spanload is shown, whether the tip vortex position is shown, and if a convergence of global coefficients is shown.

The error for $C_T$ and $C_Q$ are not necessarily given as a quantifiable number in the numerical work references and is computed as follows. For $C_T$, since the simulations are performed at the same blade collective angle as the validation, the relative error can be computed at every simulated point using the validation as the reference value and the mean error is reported. For $C_Q$, direct comparison with validation is not as simple since the $C_T$ is generally different at a given collective. The method employed in the present work to compute a quantitative error for $C_Q$ for the references of Yucekayali [12], Singh and Friedmann [14], Tugnoli et al. [17], and the previous work [5], all reported with a ^† superscript in Table 3, is as follows. The data is extracted from the $C_Q$ vs. $C_T$ plots. Using the $C_T$ as the independent variable, the trend line are constructed with 2nd degree polynomial, resulting in $R^2$ in the order of 0.9999. Then, the numerical work $C_Q$ is interpolated using the trend line at each of the validation $C_T$ value that falls between the maximal and minimal $C_T$ values simulated with the numerical work (avoiding extrapolation of the trend line). The error can then be computed at each valid validation point and finally averaged.

The UVLM of Colmenares et al. [9] is validated with the Caradonna–Tung [21] experiments. As for most numerical studies that chose that validation case, $C_T$ is compared for three different simulations, but $C_Q$ cannot be compared as it is not provided by the experiment. $C_t$ is compared with the experimental measurements without matching the global $C_T$. Therefore, the collective angle of 8 degree spanload is the closest to the experiment, because it has the closest $C_T$ of the three simulations. The $C_T$ convergence is fair, but a few more rotations would help ensuring convergence.

The UVLM of Perez et al. [10] is validated by comparing the results of a single simulation with their own experiment and with higher fidelity CFD. One should take any numerical results validated on a single simulation with caution. The error with experiments or higher fidelity methods might vary considerably with the input parameters. Validating on a single point might have a lucky or unlucky set of parameters. In this case, the error with both CFD and experiment is considerable, especially for $C_Q$, because their UVLM overestimates the thrust and under-predicts torque at that higher thrust, compounding the effect on torque vs. thrust efficiency prediction. Since their implementation does not handle viscosity, it is normal that torque prediction is inaccurate. The blade pressure coefficient is compared between CFD and UVLM, showing clear differences near the tip of the blades. Unsteady loading is visible with the darker blue patches near the trailing edge of the right UVLM blade as it is not symmetric with the left blade. A qualitative wake visualization between their CFD Q-Criterion and the UVLM vortex panel positions below the rotor is shown, but no quantitative comparison of vortex position. The $C_T$ and $C_Q$ convergence figures show some minor oscillations, but the mean values do not appear to vary too much despite still having a slight downward tendency at the end.

The VPM of Alvarez and Ning [11] is validated with the experiments of Zhou et al. [27]. That experiment has a single measurement of the isolated rotor in hover to provide the validation. Note that these works have more data points for the advance flight and multirotor cases, but those are disregarded for the present comparison. $C_T$ compares very well between the numerical and experimental studies, but the torque or the spanload cannot be compared since they are not provided by the experiment. The thrust coefficient convergence is satisfactory.

The VPM of Yucekayali [12] is first validated in Figure 4.7 of the thesis with three different references. That figure presents the $C_Q$ vs. $C_T$ efficiency of the rotor in hover. However, of the three references, only the experiments of Shinoda [24] has hover without advance ratio results, so the present analysis ignores the other two references. Unfortunately, not all the data from that experiment are shown in Figure 4.7 and the process of elimination is not explained. For the present error computation, all of the hover data from Shinoda [24] has been used, except for run 52 data, which is an obvious outlier compared with the other runs. This results in an error of about 7% for $C_Q$, with a clear tendency for the VPM to under-predict the increase of $C_Q$ as $C_T$ increases. It is worth noting that the results presented in the thesis appear to be closer at higher $C_T$ to the other references that were excluded in the present analysis. The vortex position is compared with a higher fidelity result on this same test case. The thrust convergence plot appears to be very well converged after only 5 revolutions. The thesis also compares the VPM spanload with the one measured experimentally by Caradonna–Tung [21] and obtained by other numerical methods at trimmed thrust for a single simulation (RPM = 1250, collective not given). The VPM results appear to be better than other lower fidelity methods and not as good as higher fidelity results, but the error is not quantified.

The LL-VPM of Zhao and He [13] presents a single simulation result comparison for the downwash velocity at different planes below the rotor with the experiments of Boatwright [23] and a single lift coefficient distribution along the blade with the experimental measurements of Caradonna–Tung [21], both at trimmed $C_T$. There is no convergence plot to assess numerical behavior.

The UVLM-VPM of Singh and Friedmann [14] has the most simulated validation points with 22 since it presents the rotor efficiency $C_Q$ vs. $C_T$ for the two isolated rotors tested by Harrington [22]. However, since it was the only data reported from that experiment, validation cannot be achieved for $C_T$ vs. collective, span loading, or tip vortex positioning. The error for $C_Q$ is shown as a range because it is quite different for the two rotors, that is, 11.5% for rotor 1 and 3.9% for rotor 2. The error for rotor 1 seems high for the qualitative agreement between the UVLM-VPM and the experimental results. That is because the experimental results are most populated at lower $C_T$ where the difference between UVLM-VPM and experimental results is the largest. The error would drop to 4.6% if results at $C_T<5\times {10}^{-4}$ were ignored, corresponding to $C_T$ lower than 14% of the maximum $C_T$ tested. Rotor 2 was not tested at such lower $C_T$, therefore did not suffer from that lower $C_T$ difference. Only 6 revolutions were computed.

The validation of the UVLM-VPM of Ferlisi [15] is similar to that of Colmenares et al. [9], showing a slightly lower $C_T$ error and no $C_Q$ because of the unavailability of this information from the experiment. Ferlisi [15] decided to validate the tip vortex position with the experiment, but not the spanload. The $C_T$ convergence is unambiguous.

The UPM-VPM of Tan and Wang [16] is validated only with coefficient of pressure ($C_p$) sections compared with higher fidelity and experimental results. They decided to show results at the same blade collective angle (8°) and 3 different RPM (1250, 1500, and 1750). Experimental data are provided (r/R = 0.5, 0.68, 0.80, 0.89, and 0.96) on each case. They show some sections for each case, without stating how those sections are chosen. For the RPM = 1250, r/R = 0.68 and 0.96 are shown. Three sections are shown at RPM = 1500: r/R = 0.68, 0.80, and 0.89. The RPM = 1750 shows the same sections as the RPM = 1250 simulation. The validation does not include any convergence plot.

The NL-LL-VPM of Tugnoli et al. [17] is validated with the experiments of Droandi et al. [25]. The validation is limited to the $C_P$ vs. $C_T$ efficiency of the rotor. FM vs. $C_T$ is also reported, but this is redundant information given $C_T$ and $C_P$. The experimental data is chosen for the present error estimate since the experiment follows the CFD results really well and they are plotted over them, making them easier to recover. The average error is 8.8%, where, similarly to the results of UVLM-VPM of Singh and Friedmann [14], the highest error is observed at the lowest collective. However, the error is much more evenly distributed and is not as dramatically affected by removing a few points. This time, the error remains higher than the mean error up to 42% of the maximum $C_T$. The 20 revolutions was more than enough to ensure the convergence of $C_T$.

The NL-UVLM of Samad et al. [18] validation is similar to that of Colmenares et al. [9] and Ferlisi [15]. The $C_T$ error with the experimental is similar to the others. The $C_l$ distribution compares well with the experimental data, however it does not appear to be presented at trimmed trust. For all of the 3 simulations, $C_T$ does not appear to be converged after 12 revolutions.

The NL-UVLM-VPM of Jo et al. [4] has the lowest error to the experimental data for $C_T$ prediction, but this is on a single validation point. $C_p$ contour is provided only for the sake of the refinement study. There is no convergence plot.

The NL-UVLM-VPM validation of Lee et al. [19] is similar to the others using the experiments of Caradonna–Tung [21]. The $C_T$ error was difficult to extract in the present analysis because the figure in the article has wide axes and the points are partially on top of each other. The distributed $C_t$ produced by the NL-UVLM-VPM of Lee et al. [19] is also compared with the experimental results. While the integrated $C_T$ is very good, the distributed $C_t$ does not seem to capture the dip in $C_t$ at r/R = 0.8 in the experimental results that is sometimes captured by other methods like Colmenares et al. [9], Zhao and He [13], and Samad et al. [18]. The comparison of the tip vortex position appears to be good, however the vertical axis is again quite wide. The work does not show a convergence plot.

The NL-UVLM-VPM of the previous work [5] is the only work presented in Table 3 to present both $C_t$ and $C_q$ distributions and also $C_p$ sections, as it uses the only validation that provides all of that information via higher fidelity results. It is also one of the two to provide both integrated $C_T$ and $C_Q$ as well as their convergence, the other being the work of Perez et al. [10], which has a much higher error. Most of the error for $C_T$ of the NL-UVLM-VPM of the previous work [5] was observed at lower $C_T$, whereas, similarly to the VPM of Yucekayali [12], it does not predict enough added torque at higher $C_T$.

The present work is more focused on the sensitivity study of the method since the validation with other studies was performed in the previous work. Therefore, the comparison with higher fidelity results is performed on a single point. The $C_T$ and $C_Q$ differences can be found in Section 6.2. They are reported as a range because they depend on the inputted database. Section 6.2 compares the thrust, torque, and twisting moment span loading with the higher fidelity results. The tip vortex position is compared with the Kocurek empirical formulation [30] in the sensitivity study, but is not validated with the higher fidelity results. $C_T$ and FM convergence figures are shown for each parametric change in this study.

2.4. Conclusion of the Literature Comparison

The published work presenting potential methods for solving the aerodynamics of an isolated rotor in hover is hard to compare one-to-one because of the differences between all of those work. The differences are the methods themselves (how they model the blade and the wake), the additional physics that the methods attempt to capture (compressibility and viscosity), the space and time discretization (with or without a sensitivity study to justify it), the chosen validation case (which itself differs in the measurements performed, the blades geometry, Mach and Reynolds and the number of measurement available), the chosen results type that are shown, the number of results shown, and when relevant if the results were generated at the same $C_T$ as the validation.

To remove some of these differences in the comparison, there should be a clear and consistent validation case to compare hover prediction of potential methods as it was done for the higher fidelity methods with the hover prediction workshops. There should be predefined results definition to ensure direct comparisons between methods. Such results should at least include $C_T$ and $C_Q$ over a certain range of collective angles and the spanload over the blade, ideally providing thrust and torque distributions separately. Tip vortex position over the first few rotations could also be an interesting metric. The validation case reported by Hariharan et al. [8] is the only validation case presented in Table 2 that provides all these results.

3. Method

This sections presents a summary of the previous method and the proposed improvements.

3.1. Previous Method Summary

The aerodynamic modeling of the lifting blades is accomplished with the tridimensional Unsteady Vortex Lattice Method (UVLM), a potential aerodynamic method [31]. The blades are represented by thin vortex panels placed on the camber line. The boundary conditions imposed at each time step are the no-through slip-flow on each panel and the Kutta condition at the trailing edge. At each time step, a row of vortex panels is shed at the trailing edge into the wake when the lifting geometry is displaced. This is a free-wake method, because at each time step, all the wake elements are displaced following the local flow velocity (induced plus freestream flows). This operation is by far the most computationally expensive of the UVLM because of the quadratic nature of the process on all the N elements in the wake ($O\left(N^2\right)$), where N increases at each time step.

Hovering blade rotor aerodynamic simulations are challenging for the free-wake potential methods, because the wake rolls up and clusters around the tip vortex forming below the rotor plane. The Boundary Element potential methods are also known to be singular methods, because the velocity induced by a potential element (vortex, doublet, or source) is singular when the distance of evaluation approaches zero. The singularity needs to be addressed to ensure a stable rotor simulation. As in the previous work, the Vatistas smoothing kernel [32] is used instead of the singular Biot–Savart equation for the computation of the velocity induced by a vortex where a core size $\sigma$ smooths the singularity when the distance between the straight-line vortex element and the evaluation point tends to zero (see Equation (4)).

Another issue that can arise with UVLM simulations of rotor blades, especially in the vicinity of the ground or an obstacle, is for panels control points to be convected in different directions [15]. This can lead to unrealistically stretched vortex filaments that can extend through the rotor plane or through the wake, affecting the simulation results and stability. This issue can be resolved by replacing vortex panels by their Vortex Particle Method (VPM) approximation [5,11,15,16,33,34,35,36]. The vortex particles need an additional equation to be solved for their straining, but are free to move in any direction with respect to their neighbors. The additional equation and the often increased number of particles compared with the panel wake simulations causes a significant computational cost increase. The VPM is well suited for the Fast Multipole Method (FMM) [37], an algorithm designed to reduced the quadratic computational cost ($O\left(N^2\right)$) to a linear computational cost ($O\left(N\right)$) with negligible error. Another difficulty that arises with the free-wake vortex particle simulations is the numerical instability. This happens because the VPM is inviscid, so vortex stretching does not have a counter acting mechanism to dissipate turbulent energy. Therefore, viscosity can be introduced by different methods [38] in the VPM to dissipate energy and keep the simulation stable. In this study, the viscous diffusion is modeled using the particle strength exchange (PSE) originally developed by Degond and Mas-Gallic [39] with an added eddy viscosity from LES to account for the subgrid scale [5,40,41]. The subgrid scale model (SGS) eddy viscosity developed by Vreman is used because of its highly desired properties compared with the standard Smagorinsky SGS model [6]. VPM with PSE-LES does not need arbitrary particle destruction thresholds such as the distance from the rotor plane or the time since the particle was created to keep the simulation stable.

Because of the potential flow assumptions, UVLM and UVLM-VPM simulations capture no viscous nor compressible effects, and are therefore suitable only to low Mach and very high Reynolds simulations. Fortunately, low computational cost non linear algorithms (NL-UVLM) can incorporate stripwise higher fidelity simulations results (2D or 2.5D [42,43,44]) or experimental measures to introduce locally the effects of compressibility and viscosity. While there exists the $\Gamma$ and the $\alpha$ methods, the latter is generally preferred as it is unambiguous even in post-stall situations and requires significantly less relaxation to achieve convergence in the non-linear region [2,3]. In the previous work [5], the NL-UVLM-VPM proved to be useful at moderate Reynolds (192 k $\le Re\le$ 1.09 M) and Mach (0.1 $\le M\le$ 0.65). The goal of the present study is to investigate the parameter sensitivity of the NL-UVLM-VPM at lower Reynolds (27 k $\le Re\le$ 170 k) and low Mach (M $\le 0.15$).

The complete description of the previous method used in this work can be found in ref. [5]. The fundamental mathematical formulation of the method used in this work is unchanged from the one presented in the previous work, so only pertinent equations are recalled when analyzed.

3.2. Method Improvements

Our initial implementation of the UVLM-VPM had trouble with discretization refinement. Figure 1 shows this issue, wherein keeping everything else equal, changing the time step changed the global coefficients convergence significantly. In the previous work [5], the refinement study was performed simultaneously for the blade mesh and the time step thus avoiding that issue. The refinement problem in our previous implementation was caused by an error of understanding concerning the overlapping factor.

Figure 2 shows the wake behind a blade for two different time discretizations ($\Delta \Psi =20°$ on the left and $\Delta \Psi =10°$ on the right). The blade is meshed with a rather coarse 4 × 10 panels with cosine tip refinement. There is 1 row of wake panels trailing behind the blade and then the panels start being converted into particle. As the conversion begin, only 3 sides of each panel are converted into particles, creating a “ghost” panels row with only the side next to the wake panels active (the panel sides that have particles are not contributing to the flow field). Behind the “ghost” panels row are the particles created at the previous time steps. The particles created from panels sides that were parallel to the stream direction are the trailed particles (in yellow) and the particles created between wake panels rows are the shed particles (in purple). The particle strength is proportional to the difference in the circulation of the panels on each side of the edge. Therefore, for the axisymmetric hover case, upon convergence, the shed particles are of lesser importance because their strength approaches zero. Finally, a schematic of the distance between particles near the root of the blade is shown shortly after their creation. For clarity, the schematic is only done for the trailing particles because the shed particles play a lesser role in hover, but the same logic is also true for them. The initial distance between particles is important because it sets the core size value of the particles with the overlapping condition and the core size has a major role in the stability and the accuracy of the method. As a first implementation, the distance shown in red was chosen. That is, for every trailing particle, the distance was the maximum between the streamwise and the spanwise nearest trailing particles. This seemed like a reasonable choice as it also guarantees the inclusion of the nearest shed particles and thus ensures a greater stability of the method. Unfortunately, this choice caused issues when performing the mesh independent analysis. As presented in the figure, when the time step is halved, the core size of the represented particle remains unchanged because it was set by a mesh constraint, not the time constraint. That creates a bond between the mesh and time discretization aspect ratio where the consistency of the method can only be obtained when both are refined simultaneously. That introduces an arbitrary aspect ratio constraint with the results, where changing only the time step or the blade mesh changes the simulated field. Figure 1 shows that issue where for the same geometry mesh, the time step significantly changes the coefficient of thrust convergence. Returning to Figure 2, the distance shown in green is simply the streamwise distance of the nearest trailing particles (which would be spanwise for the shed particle). Therefore, the bond between the spanwise and the chordwise direction is lifted and the particle core size is accurately reflecting the changes in that given direction, as shown on the right side of the figure. The unidirectional core size computation does not guaranty the inclusion of the neighbor shed particles, nor the spanwise nearest trailed particles, but that has not caused stability issues for the present test case using PSE-LES.

Another issue that arises when the time discretization is refined is the different wake particle spacing. For instance, Figure 3 shows a constant 8 × 20 blade mesh and 1 wake panels row before transforming sides into 1 particle with 3 different time steps ($\Delta \Psi =20°$ on left, $\Delta \Psi =10°$ in the middle and $\Delta \Psi =5°$). From left to right, the wake panels have been transformed into particles for 1, 2, and 4 iterations, respectively, therefore having a constant 20° arc length modeled with particles. As one can see, reducing the time step will cause a greater number of particles to be generated in the same area. Therefore, the effect of time step cannot be isolated directly.

The simplest solution to try to solve this issue is to increase the number of particles generated in the streamwise direction as exemplified in Figure 4, remembering that the shed particles (purple) are negligible for the hover test case. Unfortunately, this does not work so well as the particles cluster near the root can cause stability issues at larger time steps.

This unnecessary clustering can be reduced for the larger time steps as shown in Figure 5, but not for the smallest time steps because each panel side mush shed at least one particle to avoid losing vorticity. The idea is to match the tip region particle number and spacing of the finest time discretization and to gradually reduce the number of particles to keep a more even distance as the panels’ streamwise sides reduce approaching the root of the blade. The more even distribution helps re-stabilize the method at larger time steps and reduces the computational cost compared to the direct matching of the number of streamwise particles while still ensuring a finer distribution needed at the tip. Therefore, the number of particles used to replace the panels VPM parameter is instead taken as the number of tip particles (nTip) used. In Figure 5, nTip would be 4, 2, and 1, respectively. It can also be thought as a tip particle spacing ($\Delta x_{tip}$ as defined in Equation (2)) which in this case would be constant for the three simulations at 5°.

Regarding the number of wake panels rows, instead of prescribing a number of panel rows, it makes more sense to define a number of panel revolutions to keep the simulations consistent as the time step is changed. Knowing that the particles are punctual approximations of a straight line vortex element (a panel edge), it seems reasonable to keep a certain amount of panels in the blade’s wake to accurately simulate the near field vorticity in the wake, especially as it interacts with the following blade and to eventually convert the panels in particles before they experience too intense stretching. This cutoff is initially set to 2 revolutions and will be tested like the other parameters.

4. Test Case

This section presents the definition of the simplified smaller rotor test case and the generation of the 2D RANS database used by the 3D NL-UVLM-VPM.

4.1. Definition

The isolated hovering rotor test case used in this article is the EMpEROR detailed in

Table 4. EMpEROR hovering rotor test case definition.

Property	Value
Number of Blades	2
Root cut-out ratio ($r_{min}/R$)	15.79%
Radius	475 [mm]
Twist	No
Taper	No
Chord	50 [mm]
Profile	Rectangle
Thickness	3 [mm]
Rotational velocity ($\Omega$)	1000 [RPM]
Collective angle ($\theta$)	5 [°]
Reynolds number (root–tip)	26,884–170,263
Mach number (root–tip)	0.02–0.15

. EMpEROR is an experimental rotor setup that is currently in development at ISAE-SUPAERO and should provide additional validation in a future work. The geometry was chosen to be as simple as possible to facilitate numerical and experimental reproducibility. This was meant to be similar to the popular model rotor of Caradonna–Tung [21], but for smaller Reynolds close to that of UAVs rotors. The two-bladed rotor has a radius 475 mm, with a root cutout of 75 mm. The blade geometry is rectangular (no taper, sweep, nor twist) with a chord of 50 mm and a span of 400 mm, leading to an $AR$ of 8. The blades are clamped at the middle of the chord, that is the feathering axis is at the half-chord. The profile section is also rectangular with a thickness of 3 mm. The RPM is set to 1000 and the preset collective angle is 5°.

4.2. Database Generation

The databases generated for the non-linear algorithm are summarized in

Table 5. Summary of the databases generated for the non-linear algorithm.

Software	Flow Solver	Mesh	Turbulence Model
StarCCM+ (Commerical)	2D Steady Reynolds Averaged Navier–Stokes (RANS) Compressible coupled	182 k elements (quad and hex)	$k-\omega$ SST [45] (KwSST) (Fully turbulent)
			KwSST-$\gamma R{e}_{\theta }$ [46,47] (GRT) (Transitional)
CHAMPS (In-house)		92 k elements (quad)	Spalart–Allmaras [48] (SA) (Fully turbulent)
			SA Low-Reynolds [49] (SA Low Re) (Fully turbulent)
NSCODE (In-house)			Spalart–Allmaras [48] (SA) (Fully turbulent)

. Additional information can be found on the database generation in Appendix A.

5. Parametrization

This section details the parametrization of the method. To do so, each parameter identified in the introduction is varied and their effect on the results is analyzed. The parametrization begins with the UVLM that has a wake composed solely of panels. Once the UVLM parameters are determined, the vortex particles are added to the simulations to have wake panels and particles.

5.1. Parametrization of the UVLM

In this section, the parametrization of the UVLM is done for the Vatistas core size and refinement of the mesh and time. The parameters related to the vortex particles need not be analyzed for the UVLM because the wake is composed panels only.

5.1.1. UVLM Vatistas Core Size

The first step to finding the appropriate parameters for the NL-UVLM-VPM is to start with only the UVLM parameters. The three main parameters are the value of the core size used in the Vatistas smoothing kernel (the value of $\sigma$ in Equation (4)) and the discretizations, both spatial along the geometry and the time step. The Vatistas smoothing kernel [32] is:

$${\overrightarrow{u}}_P=\frac{\Gamma }{4\pi }\frac{{\overrightarrow{r}}_1\times {\overrightarrow{r}}_2}{{\left({\left|{\overrightarrow{r}}_1\times {\overrightarrow{r}}_2\right|}^{2n}+{\left|\sigma {\overrightarrow{r}}_0\right|}^{2n}\right)}^{\frac{1}{n}}}{\overrightarrow{r}}_0·\left(\frac{{\overrightarrow{r}}_1}{\left|{\overrightarrow{r}}_1\right|}-\frac{{\overrightarrow{r}}_2}{\left|{\overrightarrow{r}}_2\right|}\right)$$

(4)

where $\Gamma$ is the unknown panel circulation, ${\overrightarrow{r}}_0$, ${\overrightarrow{r}}_1$, and ${\overrightarrow{r}}_2$ are vectors set by the position of the panel and the evaluation point and $\sigma$ is the vortex core size.

First, a mesh of 8 chordwise and 20 spanwise (Sine 8 × 20) panels with a tip cosine refinement [31] and a time step allowing the blade to advance 10° ($\Delta \Psi =10°$) is chosen from previous studies [5,50] to test the core size effect on the EMpEROR geometry. The RPM ramping is initially set to 10 revolutions to avoid large vortexes and will be reduced when the computational cost will become more important. This particular case is a much smaller rotor with a less uniform blade loading at a lower collective angle than what has previously been tested, all factors that are expected to increase the sensitivity to the UVLM parameters compared to previous studies. The value of the core size was set to around 20% of the chord in the previous studies. If the core size is chosen too small, the simulations results can experience major oscillations, making the results harder to interpret. On the other hand, if the core size is taken too large, the simulations will be unnecessarily dampened, leading to too low $C_T$ and FM. The expected behavior is obtained in Figure 6 where the convergence of these two coefficients is shown with the number of revolutions of the rotor. The core size is given as a percentage of the chord of the blades. One odd observation in those figures is the distinct jump on the $C_T$ converged value when the core size is increased from 20% to 40% and from 60% to 80%, but not when the core size is increased from 40% to 60% or from 80% to 100%. FM is showing a similar behavior. It is suspected that this strange behavior is the result of numerical errors caused by a too coarse geometrical mesh paneling.

Table 6. UVLM mean and standard deviation of global coefficients over the last 10 rotor rotations for different Vatistas core size values. The difference with the chosen reference $\sigma =60\%$ is also shown.

$\mathit{\sigma }$ (%)	${\mathit{C}}_{\mathit{T}\mathit{mean}}\left({10}^{-3}\right)$	${\mathit{C}}_{\mathit{T}\mathit{std}}$ (%)	$\Delta {\mathit{C}}_{\mathit{T}\mathit{\sigma }=60\%}$ (%)	${\mathit{FM}}_{\mathit{mean}}\left({10}^{-1}\right)$	${\mathit{FM}}_{\mathit{std}}$ (%)	$\Delta {\mathit{FM}}_{\mathit{\sigma }=60\%}$ (%)
20	1.97	0.43	2.86	8.35	0.72	3.29
40	1.91	0.42	−0.31	8.08	0.62	−0.10
60	1.92	0.24	-	8.09	0.25	-
80	1.80	0.34	−6.04	7.82	0.14	−3.37
100	1.81	0.09	−5.78	7.67	0.14	−5.13

presents the mean and standard deviation for $C_T$ and FM averaged over the last 10 revolutions of each simulation. All the simulations reached 54 revolutions except $\sigma =100\%$, which only achieved 52 revolutions. The difference with $\sigma =60\%$ simulation is also shown. The simulation with $\sigma =60\%$ is chosen as the reference to continue the parametric study as it has a core size large enough to have a relatively smooth convergence of the global coefficients (standard deviation is about 0.25% for both the thrust coefficient and the figure of merit) and remains closer to the smaller core size mean values than the simulations with a larger core size value.

The tip singularities position can also be analyzed. The tip singularities position is obtained by tracking all the singularities that are released at the tip of each blade at every iteration. Their age, typically represented in degrees, corresponds to the distance the blades have traveled since that singularity was shed in the wake:

$$age=n·\Delta \Psi$$

(5)

where n is the number of time steps since the singularity was shed from the blade tip.

One has to be careful when looking at these results, because it does not thoroughly track the physical tip vortex that can be seen when looking at the velocity field under the rotor. The singularities released at the tip of the blades could very well be moving around the tip vortex as they move away from the rotor. Nonetheless, this information about the wake is easily obtained compared to the actual velocity field and rigorous tip vortex tracking. Figure 7 shows the mean non dimensional tip singularities position with and without the standard deviation.

The figures are generated using the tip singularities data of the last 10 revolutions of the rotor. On the graph is also shown the Kocurek empirical formulation [30] for the tip vortex position, computed for both the minimum and maximum $C_T$ ($1.800\times {10}^{-3}$ and $1.970\times {10}^{-3}$, respectively). The curves on the top of the figures show the radial (r/R) contraction of the tip singularities position where the tip singularities move inward as they descend in the wake. The bottom curves show the axial displacement (z/R) where the tip singularities move down into the wake with time. The figures show the expected behavior, where the tip singularities’ contraction and descent rates, given by the slope of these curves, increase as the thrust is increased (or when the core size is reduced). The axial position shows an expected acceleration after the first blade passage at 180° for this 2-bladed rotor. From these figures, it can be seen that the mean position oscillates more for the two smallest core size values and is very stable for the three highest core size values. It appears that the core size 60% matches the empirical formulation of Kocurek best for wake contraction in the first blade rotation and all core size values perform similarly for older vortexes. For the axial displacement, none of the core size values capture very well the initial tip vortex position before the first blade passage, but after that point, the three lowest core size simulations appear to have a reasonable slope and the 2 highest core size clearly under-predict axial displacement of the tip vortex. On the right figure, the standard deviation of the position shows an increasing tendency as the core size is reduced, especially as the wake gets older. These indicate stronger asymmetry in the near wake, explaining the higher variation observed in the global coefficients. Once again, the core size 60% appears to be a reasonable compromise between variance and smoothness of its results.

The spanwise distribution of the force and torque can then be compared for the different core size values. Figure 8 shows on the left the local contribution to the coefficient of thrust ($C_t$) and on the right the coefficient of thrust ($C_q$), both non dimensionnalized with the section’s local velocity to better see differences in the results for the whole span of the blade:

$$C_t=\frac{T}{\frac{1}{2}{\rho }_{\infty }{\left(\Omega r\right)}^{2}S}$$

(6)

$$C_q=\frac{Q}{\frac{1}{2}{\rho }_{\infty }{\left(\Omega r\right)}^{2}SR}$$

(7)

where T is the thrust of the section, ${\rho }_{\infty }$ the freestream density, $\Omega$ the rotational velocity, r the radial position of the section, S the surface of the section, Q the torque of the section, and R the rotor radius. $C_t$ is the same as $C_l$ of the section, but with the axis given by the rotational axis instead of normal to the local 2D stream. $C_q$ is comparable to $C_d$ with the force in the rotational plane, but also multiplied by the distance with the rotational axis to obtain a local torque contribution. The inboard part of $C_t$ is showing the expected behavior where a reduced core size value causes the vortex sheet to move faster downstream, thus reducing its induced velocity on the blade and therefore increasing the thrust in that region. The peak loading region does not seem to be well resolved. The loading shape is not smooth and it seems strange that the mid core size value results in a wide much higher peak. It is suspected that this is caused by a close interaction with the tip vortex of the preceding blade. Upon closer inspection, it can be seen that the points experiencing the maximum deviation in the peak loading area are both situated at 92.2% of the span, for the core size 60% and 80% curves. For both of these simulation, the tip vortex stays very close to the rotation plane for the first 180°, that is when the blade encounters the preceding blade tip vortex. The main difference between the core size 60% and 80% simulations is in the radial position, where for the former, after 180°, the tip vortex is around 91.7% of the blade span and for the latter, it lies around 93.0%. Therefore, in both cases, the tip vortex has a strong influence because of proximity on the section at 92.2%. For the core size 60%, the tip vortex lies inboard of the section, creating an upwash that increases the angle of attack and the thrust of the section. The opposite happens at core size 80% where the outboard tip vortex creates a downwash. For the two lowest core size cases, the tip vortex is further inboard and downward compared to the core size 60% case, thus reducing the upwash effect. The main takeaway of that figure is that it is fortuitous that the core size 60% simulation has a higher peak loading, most likely caused by a too coarse spanwise discretization for the others to capture the effect. The torque distribution follows closely the thrust distribution in its essence, but also suffers from the too coarse mesh. The torque coefficient is linearly proportional to the distance from the rotational axis and directly proportional to the drag coefficient oriented in the disk plane. Because this is an incompressible non viscous simulation, only the induced drag is computed. The torque increases as the thrust increases, except in the peak thrust region because the increase of thrust that increases the drag is partially balanced by the tip vortex upwash that reorients the aerodynamic force forward, reducing the drag proportion of the aerodynamic force in that region. Note that this is a linear simulation computing only the induced torque. When the viscous coupling is added, the increase of angle of attack that increases the coefficient of drag has a stronger contribution and will make the maximum thrust and torque coincide.

5.1.2. UVLM Mesh Refinement

The mesh convergence test is done next because it is suspected that the mesh currently used is too coarse to accurately capture the spanwise distribution of the force. The mesh convergence is carried with the 60% core size value because of its relatively low variance and smoothing in the preliminary analysis and also because it was found to be one of the most sensitive cases with the spanload results, where the tip vortex is not moving downward as fast as the lower core sizes and is less smoothed than the larger core sizes tested.

The mesh convergence can quickly become numerically expensive when dealing with the UVLM $O\left(n^2\right)$ computational cost. It is therefore decided to conduct the geometry mesh refinement on a time step of $\Delta \Psi =20°$ instead of $\Delta \Psi =10°$. This change causes a difference of less than 1% for $C_T$ and less than 2% for FM compared with the simulation of the previous section. Figure 9 shows the coefficient of thrust and the figure of merit for different mesh discretizations, where the first number is the number of chordwise panels and the second the number of spanwise panels per blade. Note that the simulation accuracy and computational complexity is highly dependent on the number of spanwise panels and almost unaffected by the number of chordwise panels. Those figures clearly show that the mesh independent results for the global coefficients requires a mesh of at least 24 × 60, thus confirming the suspicion of too coarse mesh for the initial core size analysis.

Table 7. UVLM mean and standard deviation of global coefficients over rotations 60 to 80 for different blade mesh. The difference with the finest mesh is also shown.

Mesh	${\mathit{C}}_{\mathit{T}\mathit{mean}}\left({10}^{-3}\right)$	${\mathit{C}}_{\mathit{T}\mathit{std}}$ (%)	$\Delta {\mathit{C}}_{\mathit{T}\mathit{fine}}$ (%)	${\mathit{FM}}_{\mathit{mean}}\left({10}^{-1}\right)$	${\mathit{FM}}_{\mathit{std}}$ (%)	$\Delta {\mathit{FM}}_{\mathit{fine}}$ (%)
4 × 10	1.90	0.55	3.87	8.24	0.65	−0.42
8 × 20	1.90	0.26	3.53	8.23	0.26	−0.46
12 × 30	1.77	0.14	−3.47	8.42	0.19	1.81
16 × 40	1.80	0.15	−1.63	8.32	0.29	0.56
24 × 60	1.83	0.10	0.15	8.28	0.31	0.13
32 × 80	1.84	0.21	0.47	8.26	0.14	−0.18
48 × 120	1.83	0.11	-	8.27	0.16	-

shows the average, standard deviation, and difference with the finest discretization for the coefficient of thrust and the figure of merit, averaged between rotations 60 and 80.

The mean coefficient of thrust varies rapidly with discretization for the coarser meshes, going from a 3.53% overprediction with a mesh of 8 × 20 to an underprediction of −3.47% with a mesh of 12 × 30, that is, with only 1.5 as many panels. The $C_T$ difference is much closer for the finer meshes where doubling the number of panels from 24 × 60 to 48 × 120 only causes a difference of $C_T$ of 0.15%. A slightly larger difference is found between 32 × 80 and 48 × 120, but that is still under half a percent. The figure of merit reveals a similar story but is apparently less sensitive on the mesh. For both coefficients, the standard deviation tends to be slightly larger for the coarsest meshes and smaller for the finest mesh with few exceptions.

Figure 10 shows the mean span distribution of the thrust and torque with their standard deviation again for the revolutions 60 to 80. It is very interesting to see that the over-prediction of the peak thrust for the 8 × 20 mesh, as it was found for the 10° time discretization, is not amplified with the mesh refinement. As the mesh is refined, the distribution of the coefficient of thrust converges towards the finest mesh and the deviation is generally reduced, except for the 32 × 80 that experiences more deviation than the 24 × 60 and 48 × 120 near the peak loading. This might be caused by some interaction with the previous tip vortex as highlighted in the previous section. Still, the three finest meshes have a distribution of thrust that is almost indistinguishable. The distribution of torque shows the 2 coarsest meshes over-predict the torque in the drop area after the peak and the 12 × 30 and 16 × 40 meshes under-predict it. The 24 × 60 is hidden under the 32 × 80 and 48 × 120 meshes for most of the span especially after the peak torque.

The tip vortex is not shown here as it is not very affected by the change of discretization, aside that the deviation is reduced as the mesh is refined, as for the spanload. The panel mesh convergence was also confirmed using the NL-UVLM with the 16 × 40, 24 × 60 and 32 × 80 meshes and reducing the ramp up revolutions from 10 to 2 to reduce the computational cost. The 24 × 60 mesh is again fine enough with those changes.

5.1.3. NL-UVLM Time Refinement

The non-linear algorithm is now added to the UVLM algorithm for the time refinement analysis. It is hypothesized that the non linear algorithm has a marginal effect on the parameter selection of the UVLM, because it only modifies the thrust distribution on the blades. In this case, it increases $C_T$ by 10% and modifies the thrust distribution by at most 15% over the span of the blade. The discretization and wake parameters appropriate for either UVLM or NL-UVLM should also be appropriate for the other. Note that the more significant added drag by the non linear algorithm is only taken into account in the output, thus not affecting the core UVLM wake structure.

The NL-UVLM time refinement is conducted on the previously selected 24 × 60 mesh. Only 3 time discretizations are tested, because coarser than 20° per step ($\Delta \Psi =20°$) seems too prohibitive and the 5° per step ($\Delta \Psi =5°$) is already very computationally expensive with the $O\left(n^2\right)$ complexity (already 2483 iterations done, with iterations at the end of the simulation close to 1090 s on 40 cores). Figure 11 shows the convergence of the coefficient of thrust and the figure of merit for the 3 simulations. It can be seen in those figures that the finest case has been stopped at a questionably converged state, because of the cost of running the simulation longer. Therefore, the comparison with the finest case is to be taken lightly. Nonetheless, both the $\Delta \Psi =20°$ and $\Delta \Psi =10°$ are well converged.

Table 8. NL-UVLM mean and standard deviation of global coefficients over rotations 29 to 34 for different time discretizations. The difference with the finest time is also shown.

$\Delta \Psi$ (°)	${\mathit{C}}_{\mathit{T}\mathit{mean}}\left({10}^{-3}\right)$	${\mathit{C}}_{\mathit{T}\mathit{std}}$ (%)	$\Delta {\mathit{C}}_{\mathit{T}\mathit{fine}}$ (%)	${\mathit{FM}}_{\mathit{mean}}\left({10}^{-1}\right)$	${\mathit{FM}}_{\mathit{std}}$ (%)	$\Delta {\mathit{FM}}_{\mathit{fine}}$ (%)
20	2.02	0.10	−3.80	1.02	0.14	−5.01
10	2.06	0.28	−1.73	1.05	0.39	−2.31
5	2.10	0.42	-	1.07	0.59	-

compares the mean and standard deviation for the coefficient of thrust and the figure of merit, for the revolutions 29 to 34.

The higher deviation for the finest case is expected as the simulation would have ideally been converged for a longer time. The deviation of the $\Delta \Psi =10°$ could have been reduced by averaging the values slightly further, where the coefficients reach a better convergence after revolution 32. The trend seems to be consistent in this table, where halving the time step seems to increase the mean coefficient of thrust and the figure of merit in the order of 2% and 2.5%, respectively.

Analysis of the spanwise distribution of force for the time refinement case must be accomplished carefully, because the spanload can take longer than the overall coefficients to stabilize and the smallest $\Delta \Psi$ simulation did not achieve a full convergence of the spanload. To illustrate this, Figure 12 shows the spanwise distribution of thrust for different time steps and averaging windows. The $\Delta \Psi =20°$ has results up to revolution 58, $\Delta \Psi =10°$ up to 48 and $\Delta \Psi =5°$ up to 34. Therefore, all time steps can be compared for the 29–34 revolutions window but only $\Delta \Psi =20°$ and $\Delta \Psi =10°$ can be compared at revolutions 38–48 and so on. Focusing on the peak area, the $\Delta \Psi =20°$ rev 29–34 has the lowest peak. Then, 3 simulations have a similar peak: $\Delta \Psi =20°$, rev 38–48 and $\Delta \Psi =5°,10°$, rev 29–34. Finally, 2 simulations predict similarly the highest peak: $\Delta \Psi =20°$, rev 48–58 and $\Delta \Psi =10°$, rev 38–48. This means that the peak loading is increased as the revolutions increase and as the time step is better resolved. For both averaging windows that $\Delta \Psi =20°$ and $\Delta \Psi =10°$ share (rev 29–34 and rev 38–48), they are superimposed just inboard of the peak, in the 83–90% of blade span region. Finally looking at the inboard region, at a given $\Delta \Psi$, the thrust is increased as the flow converges, but the increase remains smaller than the offset caused by the change in the time step. The magnitude of the variation shown here is rather small (about 3% between the peak min and max) and that is why it does not affect the overall coefficient so much, but is still to bear in mind when comparing span distribution of force.

The distribution of torque is not shown here because it did not show enough differences to be readable. It will be shown in the next section. The tip vortex position is also not shown because it showed little differences with respect to the chosen time step. For computational reasons, the time step of $\Delta \Psi =10°$ is chosen for the NL-UVLM and the time step will be analyzed again with the addition of the particles.

5.1.4. NL-UVLM Vatisas Verification

The effect of core size is revisited with the finer discretization and the non linear coupling. More specifically, should the Vatistas core size value be reduced now that an appropriate blade and time discretization have been selected?

Figure 13 shows the global coefficients with different Vatistas core size values. Simulations with core size 60% and 40% appear to converge better than the simulations with the same core size values previously shown in Figure 6 with the longer ramp up phase, the coarser discretization and without the non linear coupling. Indeed, the standard deviation for those core size values reported in

Table 9. NL-UVLM mean and standard deviation of global coefficients over the last 10 rotor rotations for different Vatistas core size values. The difference with the chosen reference $\sigma =60\%$ is also shown.

$\mathit{\sigma }$ (%)	${\mathit{C}}_{\mathit{T}\mathit{Mean}}\left({10}^{-3}\right)$	${\mathit{C}}_{\mathit{T}\mathit{Std}}$ (%)	$\Delta {\mathit{C}}_{\mathit{T}\mathit{\sigma }=60\%}$ (%)	${\mathit{FM}}_{\mathit{Mean}}\left({10}^{-1}\right)$	${\mathit{FM}}_{\mathit{std}}$ (%)	$\Delta {\mathit{FM}}_{\mathit{\sigma }=60\%}$ (%)
0	2.26	0.34	8.94	1.18	0.47	11.83
20	2.16	0.58	4.40	1.11	0.76	5.87
40	2.14	0.35	3.44	1.10	0.47	4.52
60	2.07	0.12	-	1.05	0.16	-

, once again averaged over the last 10 revolutions, is now lower than that previously shown in Table 6.

However, the differences in the results between the simulations with Vatistas core size value 60% and with the lower core size values has increased compared to previously. For example, when comparing with core size value of 20%, the differences for $C_T$ and FM are now 4.40% and 5.87%, respectively, when previously they were 2.86% and 3.29%. The simulation with Vatistas core size value 40%, which previously was close to the 60%, is now closer to the 20%.

As it can be seen in Figure 14, the spanwise thrust distribution now shows a more consistent behavior with core size variation. The smaller the core size, the faster the wake is convected downward. This reduces the downwash in the inboard part of the blade, leading to a higher aerodynamic force in that region. It also explains the changes in the peak loading, where a farther tip vortex will create less upwash in the tip region, leading to a flattened peak. The 20% core size is the simulation showing the most deviation in the peak region. Both 40% and 60% core size values show a well behaved peak loading, narrowing as the core size is increased. The deviation in the inboard region consistently increases with the reduction of the core size. It is quite impressive how much smaller the standard deviation is for the 60% core size in the first half of the blade, especially near the root. The torque distribution now follows more closely the thrust distribution as the non linear algorithm adds the profile 2D drag coefficient containing the pressure and skin friction components that are more impactful than the induced drag. For a very similar peak $C_t$, the peak $C_q$ is almost one order of magnitude larger than that computed only with the induced drag as shown in Figure 8. Therefore, the region of maximum thrust is also the region of maximum torque, as the tip vortex upwash increases the angle of attack, which in turn increases the drag coefficient more significantly than the forward rotation of the thrust force.

The tip vortex is not shown as it has a very similar behavior as the core size test conducted with the 8 × 20 mesh.

The main takeaway from this section is to show that the spanwise loading, just like the global coefficients prediction, has a predictable response to core size, given that the mesh is sufficiently discretized. The goal of this study is not as much to determine the optimal core size value that should be taken in all circumstances, but to highlight the effects of the choice of the core size value that could become a tuning parameter depending on the physics wanting to be captured. The study is continued with the core size 60% for its increased stability, reduced deviation, and also because knowing from previous experience, the NL-UVLM and NL-UVLM-VPM tend to over-predict the thrust at a given collective angle.

5.1.5. Summary of the UVLM Parametrization

The first parameter analyzed is the Vatistas core size given a baseline discretization of 8 × 20 panels on the geometry and a time step of 10°. As expected, the $C_T$ and FM convergence are smoothened and lowered with an increasing core size. The Vatistas core size of 60% appears as a reasonable compromise, but the distribution of thrust seems to suggest that the mesh discretization is too coarse to accurately capture the peak loading.

Therefore, the second parameter analyzed is the geometry mesh discretization. As suspected, the previous 8 × 20 mesh discretization was too coarse since the results clearly show that a mesh of 24 × 60 is necessary to become mesh independent.

Finally, the third parameter tested is the time step. Unfortunately, the simulation using the time step of 5° is stopped after 34 revolutions which was before reaching full convergence because of the computational cost of that simulation. Nonetheless, the trend indicates that the time step of 10° appears sufficiently fine to be reasonably converged with respect to time step refinement.

Finally, the effect of Vatistas core size is revisited with the refined mesh discretization and the same core size value appears to be appropriate.

5.2. Parametrization of the UVLM-VPM

The previous UVLM parametric analysis was done on a wake composed only of wake panels. As it will be seen in this section, there are some advantages to transform the vortex panels into particles. This conversion introduces new parameters that also need to be accounted for. Namely, there is the number of wake panels rows that are kept before being transformed into particles, the number of particles used to replace the panels, and the Vreman model coefficient if PSE-LES is used to keep the simulations stable. Note that the core size of the particles is fully defined by the spacing of the particles using the well known overlapping factor. If that factor is set below unity, the core size of the particles does not encompass their neighbor and the simulation quickly becomes unstable. As is done in most of the VPM work, the overlapping factor is set to 1.3.

5.2.1. UVLM-VPM Vreman Model Coefficient Stability at Constant $\Delta x_{tip}$

In this work, the vortex particle stretching is stabilized by adding a turbulent LES viscosity from the Vreman model to the fluid viscosity in the viscous diffusion computed with the Particle Strength Exchange. The Vreman eddy viscosity [6] is defined as:

$${\nu }_T=C_v\sqrt{\frac{{\Pi }_{\beta }}{{\overline{\alpha }}_{ij}{\overline{\alpha }}_{ij}}}$$

(8)

where ${\overline{\alpha }}_{ij}=\frac{\partial {\overline{u}}_j}{\partial x_i}$ is the tensor that represents the gradient of the filtered velocity $\overline{u}$ and ${\Pi }_{\beta }={\beta }_{11}{\beta }_{22}-{\beta }_{12}^2+{\beta }_{11}{\beta }_{33}-{\beta }_{13}^2+{\beta }_{22}{\beta }_{33}-{\beta }_{23}^2$ is the invariant of the tensor ${\beta }_{ij}={\sum }_{m=1}^3{\overline{\Delta }}_m^2{\overline{\alpha }}_{mi}{\overline{\alpha }}_{mj}$. The filter width $\Delta$ is chosen to be the particle core size $\sigma$. In this formulation, the Vreman model coefficient $C_v$ is a loose parameter that needs to be set. In the previous study [5], it was set to 0.07 without testing the sensitivity of that parameter. That value was chosen because it is the equivalent value to the theoretical value of the Smagorinsky constant for homogeneous isotropic turbulence [6]. It was somehow fortuitous that the first tested Vreman model coefficient value was able to stabilize the simulations. Unfortunately, this value is unstable for the current validation case.

Even though the previous section suggested that the mesh and time convergence for the UVLM required discretization of 24 × 60 and $\Delta \Psi =10°$, respectively, the sensitivity of the UVLM-VPM to the Vreman model coefficient is first investigated on a coarser discretization (8 × 20 and $\Delta \Psi =20°$) to first gain insight on the stability and behavior with respect to that parameter at lowered computational cost. The following analysis has been conducted for $\Delta x_{tip}=20°$ (nTip = 1) and $\Delta x_{tip}=5°$ (nTip = 4) to extract the general trend. The $\Delta x_{tip}=20°$ is not shown here because it never achieved a clean convergence regardless of the Vreman model coefficient value, so it was hard to interpret, but suggested a value in the order of 0.007 to reduce the oscillations. It is believed that approximating a 20° long straight line vortex element by a single particle is simply too coarse a discretization to obtain convergence. The $\Delta x_{tip}=5°$ is thus more interesting and is detailed here.

First, an order of magnitude test is conducted to have an idea of the values that might be well behaved. Figure 15 shows the convergence of $C_T$ and FM for different Vreman model coefficient values. As it can be expected, $C_v$ = 0.0, which is equivalent to no LES viscosity in the PSE, quickly becomes unstable (around 15 revolutions). Increasing the Vreman model coefficient consistently increases the stability region, keeping the simulation stable up to 40 revolutions at $C_v$ = 0.005 and 50 revolutions at $C_v$ = 0.011. However, starting at $C_v$ = 0.05, the converged global coefficients gradually deviate from the lower $C_v$ values near revolution 30 and very large $C_v$ causes an over-prediction and oscillations of the global coefficients. It is interesting to notice that the $C_v$ = 0.005 and $C_v$ = 0.11 yield very similar results in their stability region. Another interesting observation is that when the Vreman model coefficient is chosen correctly, the use of UVLM-VPM removes the long panel convergence causing oscillations in the global coefficients for UVLM (panel only) in the 10–30 revolutions region, therefore needing less revolutions to obtain converged results. An additional advantage of the VPM formulation is the reduction of the computational cost from $O\left(n^2\right)$ to $O\left(n\right)$ with the use of FMM. The VPM requires more equations to be solved, making the initial computational cost higher than UVLM, but becomes much faster as n becomes large, that is after a high number of revolutions.

A finer analysis is then conducted for Vreman model coefficient values in the range $0.011<=C_v<=0.050$ to see if the stability can be further increased without affecting the results as much. Figure 16 shows the same coefficients with a finer Vreman model coefficient delta between the curves. The important observation is to see that the trend is the same as the previous figures, where increasing the Vreman model coefficient consistently increases stability and if set too large eventually starts affecting the results. The

Table 10. UVLM-VPM tally of the simulation instability in revolutions for different Vreman model coefficient values ($C_v$).

${\mathit{C}}_{\mathit{v}}\left({10}^{-2}\right)$	0.0	0.3	0.5	0.7	0.9	1.1	1.5	2.0	>= 2.5
Instability (Rev)	13	36	40	42	49	51	52	61	>70

summarizes the stability margin of the simulations for the different $C_v$ values tested. The exact instability point of $C_v>=2.5\times {10}^{-2}$ has not been found despite running up to 70 revolutions.

Table 11. UVLM-VPM standard deviation of global coefficients over rotations 60–70 for different Vreman model coefficient values. The difference with the UVLM (Panel only wake) reference is also shown.

${\mathit{C}}_{\mathit{v}}\left({10}^{-2}\right)$	${\mathit{C}}_{\mathit{T}\mathit{std}}$ (%)	$\Delta {\mathit{C}}_{\mathit{T}\mathit{UVLM}}$ (%)	${\mathit{FM}}_{\mathit{std}}$ (%)	$\Delta {\mathit{FM}}_{\mathit{UVLM}}$ (%)
2.5	0.11	0.10	0.14	2.94
2.8	0.12	0.08	0.14	2.91
3.0	0.10	0.32	0.11	3.14
3.5	0.10	0.65	0.11	3.47
5.0	0.08	1.04	0.10	3.82

shows the standard deviation and the difference in mean $C_T$ and FM with respect to the panel reference case, where the average is performed over revolutions 60–70 for the stable cases. The use of UVLM-VPM reduces the deviation by about half compared to the panel case (previously about 0.25% for both $C_T$ and FM), and the $C_T$ compares very well for the lowest stable Vreman model coefficients with the panel reference. However, the induced drag is under-estimated leading to over-estimation of FM of about 3% for the lowest stable Vreman model coefficients. With the use of the non linear algorithm, that difference will become much smaller in the next section when the viscous effects are taken into account. It is impressive that the stability region of the UVLM-VPM is increased by PSE-LES from 14 revolutions to more than 70 with a relatively large time step (which typically reduces the stability) without relying on any arbitrary threshold like destroying the particles after an arbitrary cutoff while giving results in very good agreement with UVLM. The disadvantage of having an additional parameter to parametrize with PSE-LES is compensated for by the increased stability and result agreement.

Figure 17 shows the spanwise distribution of the coefficient of thrust and torque for the different Vreman model coefficient values, averaged over rotations 40–50 where the Vreman model coefficients 0.011 and larger are stable. The increase in the Vreman model coefficient tends to slightly reduce the peak loading and increase the inboard force. The standard deviation seems to be about the same for most cases. The differences are deemed small considering the range of Vreman model coefficients (almost 5 times larger from the smallest to the largest Vreman model coefficient).

The differences in the peak loading might be explained by the convergence of the flow field. As the Vreman model coefficient is increased, the more viscosity dissipates energy in the wake, especially in the starting vortex, increasing the time it takes to get a fully converged flow. Figure 18 shows the spanwise loading averaged over revolutions 40–50 and 60–70 for the smallest stable Vreman model coefficient tested over theses ranges and for the $C_v$ = 0.050 identified as the upper bond from the order of magnitude analysis. The UVLM spanload from Section 5.1.1 labeled “Panel” is also shown as a reference. The spanload of the simulations $C_v$ = 0.011, Rev = 40–50 (red) and $C_v$ = 0.025, Rev = 60–70 (green) are very close to the panel reference in the inboard part of the blade, except maybe at the root where the UVLM-VPM simulations capture an increased load. The $C_v$ = 0.011, Rev = 40–50 (red), despite having the highest peak $C_t$ in the previous Figure 17, is under-predicting the peak loading slightly compared with the panel reference. This might be because the spanload is not fully converged at revolutions 40–50, as it can be suspected when looking at a larger Vreman model coefficient simulation stable for more revolutions. For instance, the $C_v$ = 0.050 (teal) has a smaller peak loading than the $C_v$ = 0.011 (red) for the same averaging window (rev = 40–50), because it dissipates more energy via PSE-LES, thus slowing convergence. When the averaging is performed for revolutions 60–70, the peak loading of the $C_v$ = 0.050 (blue) is increased from the lowest to the highest on that plot. It is interesting to notice that both simulations averaged over revolutions 60–70 are superimposed in the peak region, converging to a slightly increased peak loading compared with the panel reference. The higher Vreman model coefficient appears to affect more the inboard part of the blade where the $C_v$ = 0.050 is over-predicting the loads in that region regardless of the averaging revolutions compared to the panel reference. All UVLM-VPM simulations induced torque distribution agree well with the panel reference for most of the blade except in the peak thrust region where it is under-estimated.

Figure 19 shows the position of the tip singularity as a function of age (as defined in Equation (5)) for the different Vreman model coefficients and the panel reference. For the UVLM-VPM cases, the panels are converted into particles after 2 rotations (720°), explaining the sudden kink in the radial position. It can be seen that the UVLM-VPM cases have a very similar behavior regardless of the Vreman model coefficient and predict a tip singularity position in very good agreement with the panel reference for the first 1.5 revolutions (540°). The differences in the near field of the particles compared to that of panels has an increasing effect approaching the replacement location. The VPM wake region is acting similarly to the panels with a higher core size, having a slower descent rate and a smaller standard deviation than the panel reference.

The $C_v$ = 0.028 is chosen for this $\Delta x_{tip}=5°$ case because it is stable and leads to the smallest differences with the panel reference case for the global coefficients. It seems that the Vreman model coefficient has to be increased as the particle spacing is decreased ($C_v$ = 0.028 for $\Delta x_{tip}=5°$ and about 0.0007 for $\Delta x_{tip}=20°$). Noting that there is a factor 4 for both the Vreman model coefficient and the spacing, it is hypothesized that there is a linear relationship between the two such that:

$$C_v=C_{v,min}\frac{\Delta x_{tip,max}}{\Delta x_{tip}}$$

(9)

where $C_{v,min}=0.007$ and $\Delta x_{tip,max}=20.0$ for this specific configuration. It cannot be stressed enough that the value of $C_{v,min}$ is case specific and expected to be very different in other contexts. The simulation $\Delta \Psi =20°$, $\Delta x_{tip}=10°$ is tested at the predicted Vreman model coefficient of 0.014, leading to a simulation still stable at the end of 86 revolutions. Halving that constant (0.007) leads to very similar results for the first 70 revolutions before it starts oscillating and diverges. The same was done at $\Delta \Psi =10°$ for spacing $\Delta x_{tip}=10°,5°$, where the predicted Vreman model coefficient is compared with this value doubled and halved, to find that the halved value was too small to keep the oscillations damped after a long time and the doubled value was starting to affect the results. This means that the choice of Vreman model coefficient depends on the particle distance, not the time step. The takeaway here is that knowing the correct Vreman model coefficient for a given particle spacing, the Vreman model coefficient is easily predicted for different spacing, and also that the results are not too sensitive on the exact Vreman model coefficient choice as long as it is chosen in the neighboring of the correct value.

The same verification is finally conducted with the NL-UVLM-VPM 24 × 60 blade mesh at $\Delta \Psi =20°$, $\Delta x_{tip}=20°,10°$ and $\Delta \Psi =10°$, $\Delta x_{tip}=10°$ confirming the appropriateness of the predicted Vreman model coefficient values with a finer mesh.

5.2.2. NL-UVLM-VPM Tip Vortex Particle Refinement ($\Delta x_{tip}$) at Constant Time Step ($\Delta \Psi$)

Using the NL-UVLM-VPM with the 24 × 60 blade mesh and the Vreman model coefficient predicted from the previous section, the effect of the number of particle is now investigated for $\Delta \Psi =20°$ and $\Delta \Psi =10°$. Because of the finer blade mesh that increases the computational time, the ramp up revolutions are limited to 2 and the total number of revolutions is reduced to around 30 to 40 depending on the computational complexity. The averaging is reduced to 2 revolutions instead of 10, to avoid taking data that are not yet converged.

Figure 20 shows the convergence of the global coefficients with $\Delta \Psi =20°$ for different $\Delta x_{tip}$. All simulations predict very similar maximum $C_T$ and FM at the end of the ramp up revolutions, but the farther the particles are (larger $\Delta x_{tip}$), the more pronounced the undershoot and the subsequent overshoot of the coefficients are before stabilizing on a converged value. Except for the $\Delta x_{tip}=20°$ that struggles to converge to a stable value, the other 3 simulations all converge to very similar final values. As $\Delta x_{tip}$ is reduced, the more particles are placed in the tip vicinity and the less rotations the convergence appear to require.

Figure 21 shows the same plot but this time with a $\Delta \Psi =10°$, which explains why the curve $\Delta x_{tip}=20°$ is not present. It is interesting to note the similarities between the same $\Delta x_{tip}$ at the different $\Delta \Psi$ in both the convergence and the final value. A given $\Delta x_{tip}$ seems to be slightly more stable and smooth in the second overshoot for $\Delta \Psi =10°$ compared to the $\Delta \Psi =20°$.

Table 12. NL-UVLM-VPM standard deviation of global coefficients over rotations 35–37 for different tip vortex particle spacing. The difference with the finest $\Delta x_{tip}$ and the UVLM (panel only wake) reference is also shown.

$\Delta {\mathit{x}}_{\mathit{tip}}$	${\mathit{C}}_{\mathit{T}\mathit{std}}$ (%)	$\Delta {\mathit{C}}_{\mathit{T}\mathit{fine}}$ (%)	$\Delta {\mathit{C}}_{\mathit{T}\mathit{UVLM}}$ (%)	${\mathit{FM}}_{\mathit{std}}$ (%)	$\Delta {\mathit{FM}}_{\mathit{fine}}$ (%)	$\Delta {\mathit{FM}}_{\mathit{UVLM}}$ (%)
20	0.14	−2.90	−3.14	0.19	−3.53	−3.74
10	0.01	0.38	0.13	0.01	0.52	0.30
5	0.02	0.25	−0.01	0.03	0.31	0.10
2.5	0.03	-	−0.25	0.04	-	−0.22

shows the standard deviation and the difference with both the finest $\Delta x_{tip}$ discretization and the UVLM reference from Section 5.1.3 for the global coefficients on the $\Delta \Psi =20°$ simulation. The standard deviation is much smaller than previously partly because the blade mesh is better refined and also because the averaging window is only 2 revolutions and the coefficients are varying slowly in the convergence (no distinct high frequency). Except for the $\Delta x_{tip}=20°$, which is simply too coarse, the other refinements show very close agreement on the converged results with respect to the finest case. Except for $\Delta x_{tip}=20°$, the agreement with the UVLM reference is still very good for $C_T$ and is now much better than that of the previous section for FM with the addition of the non linear algorithm. Very similar results are observed with $\Delta \Psi =10°$.

The convergence of the spanload with respect to the number of revolutions is then investigated. Figure 22 shows the $C_t$ distribution for the $\Delta \Psi =20°$, $\Delta x_{tip}=10°$ in the last 10 revolutions (revolutions 34–44), averaged over 5 windows of 2 revolutions each. As it can be seen in the figure, apart from the peak loading that varies about 1% over those 10 revolutions, the rest of the span loading appears to be well converged with time.

The convergence of the spanload is then assessed with the particle distance. The spanwise distribution of thrust and torque in Figure 23 are averaged between rotations 35 and 37 for NL-UVLM-VPM $\Delta \Psi =20°$ and the panel reference. As the distance between the particles is decreased, the peak loading is consistently increased showing convergence towards the panel reference in the tip region. The simulation with $\Delta x_{tip}=2.5°$ is even predicting a slightly higher peak than the panel simulations, but the reader is reminded that the spanload for the panels simulations was still increasing slightly with time. As the number of particles is increased, the simulations appear to capture another lower peak near the root, as the panel $\Delta \Psi =5°$ did. Where the global coefficients suggested that $\Delta x_{tip}=10°$ is sufficiently fine to capture the overall coefficient, these figures show that $\Delta x_{tip}=10°$ is too coarse to capture accurately the details of the spanwise loading.

Figure 24 shows the mean position of the tip singularity averaged between rotations 35 and 37 for NL-UVLM-VPM $\Delta \Psi =20°$ and the panel reference. As the tip distance is reduced, the axial and radial position of the tip singularity tends towards the panel reference. The two finest discretizations are showing oscillations in the tip singularity position, but these oscillations are reduced at smaller time step.

The takeaway of this section is that the NL-UVLM-VPM converges toward the NL-UVLM reference as the number of simulated revolutions increases and as the number of particles placed at the tip is increased. Also, converged global coefficients is not a guaranty of converged distributed spanload.

5.2.3. NL-UVLM-VPM Mesh Convergence

As it was done for the UVLM mesh convergence, the mesh convergence for the UVLM-VPM is conducted on the $\Delta \Psi =20°$ time discretization. Also, to reduce the computational complexity for the finest mesh, the study is conducted with $\Delta x_{tip}=10°$, even thought the span loading peak might be slightly under-predicted.

Figure 25 and

Table 13. NL-UVLM-VPM standard deviation of global coefficients over rotations 37–39 for different blade mesh. The difference with the finest mesh is also shown.

Mesh	${\mathit{C}}_{\mathit{T}\mathit{std}}$ (%)	$\Delta {\mathit{C}}_{\mathit{T}\mathit{fine}}$ (%)	${\mathit{FM}}_{\mathit{std}}$ (%)	$\Delta {\mathit{FM}}_{\mathit{fine}}$ (%)
4 × 10	0.03	−0.06	0.04	0.03
8 × 20	0.02	−0.32	0.02	−0.48
16 × 40	0.01	0.12	0.01	0.15
24 × 60	0.02	−0.05	0.02	−0.07
32 × 80	0.02	-	0.02	-

compare the global coefficients for different blade mesh discretizations. It can be seen that the convergence is a little bit different for the coarsest meshes, but that the converged results are very close for all cases. The deviation and the difference with the finest case is much lower than for the panel case (see Section 5.1.2). Figure 26 also confirm that the spanwise variation of thrust and torque are negligible for the three finest meshes and still close to the coarsest discretization. This is in clear contrast with the UVLM (panel) case where a mesh of 24 × 60 was needed to ensure mesh convergence. The UVLM-VPM is therefore much less sensitive to mesh coarsening than UVLM. The mesh 24 × 60 is still kept in this study to ensure direct comparison with the panel simulations.

5.2.4. NL-UVLM-VPM Time Convergence

The NL-UVLM-VPM time step convergence is done at constant tip distance ($\Delta x_{tip}$). The mean coefficients are compared in Figure 27 for $\Delta x_{tip}=5°$ and in

Table 14. NL-UVLM-VPM standard deviation of global coefficients over rotations 40–42 for different time discretizations with $\Delta x_{tip}=10°$. The difference with the finest time discretization is also shown.

$\Delta \Psi$	${\mathit{C}}_{\mathit{T}\mathit{std}}$ (%)	$\Delta {\mathit{C}}_{\mathit{T}\mathit{fine}}$ (%)	${\mathit{FM}}_{\mathit{std}}$ (%)	$\Delta {\mathit{FM}}_{\mathit{fine}}$ (%)
20	0.01	−0.89	0.01	−1.14
10	0.05	-	0.06	-

and

Table 15. NL-UVLM-VPM standard deviation of global coefficients over rotations 26–28 for different time discretizations with $\Delta x_{tip}=5°$. The difference with the finest time discretization is also shown.

$\Delta \Psi$	${\mathit{C}}_{\mathit{T}\mathit{std}}$ (%)	$\Delta {\mathit{C}}_{\mathit{T}\mathit{fine}}$ (%)	${\mathit{FM}}_{\mathit{std}}$ (%)	$\Delta {\mathit{FM}}_{\mathit{fine}}$ (%)
20	0.06	−2.18	0.07	−2.83
10	0.02	−1.07	0.02	−1.40
5	0.03	-	0.04	-

for $\Delta x_{tip}=10°,5°$, respectively. Whether $\Delta x_{tip}=10°$ or 5°, the refinement study shows a consistent increase of about 1% for $C_T$ and 1.4% for FM when the time step is halved.

Figure 28 shows the spanwise distribution of thrust and torque differences at $\Delta x_{tip}=5°$. Finer time steps lead to an increased inboard thrust, as if it is converging in fewer revolutions than the larger time steps. Only $\Delta \Psi =20°$ captures a peak loading that is slightly too low. The overall agreement between the 3 time step is very good, especially for torque prediction.

The tip singularity position is not shown here because it captures very little difference between the time steps, aside that the oscillations for the 20° vanish at smaller time step.

It is interesting to obtain such time convergence at $\Delta \Psi =10°$ without having to change the blade mesh, where in the previous work [5], the time step was limited to $2.5°$ because of the lack of wake panels behind the blade and the time convergence was only achieved keeping a constant aspect ratio of the wake panels (as explained in the beginning of this section). The UVLM-VPM is even showing a cleaner time convergence than the UVLM without VPM, which was also hindered by its enormous computational cost alleviated in the VPM with the use of FMM.

5.2.5. NL-UVLM-VPM Wake-Particle Conversion Revolution

Figure 29 shows the effect of the wake-particle conversion revolution on the global coefficients. As it can be expected, increasing the conversion revolution makes the UVLM-VPM converge towards the UVLM (conversion = ∞) results.

The stability of the simulations was found to depend on the convert time. The

Table 16. NL-UVLM-VPM tally of the simulation instability in revolutions for different numbers of wake panels rotations before being transformed in particles.

Conversion (Rev)	0.0	0.5	1.0	2.0	4.0	∞
Instability (Rev)	34	12	∼25	>42	>32	>48

summarizes the stability of the simulations as a function of the convert revolution. The earliest numerical instability was found at half a rotation of wake panels (conversion = 0.5 Rev), that is conversion from wake panels to particle at the first blade passage. The simulation with conversion after 1 revolution, or at the second blade passage, was stopped because it was slowly diverging unlike the other unstable cases that rapidly become unstable in typically less than 1 revolution. The simulation converting instantaneously the wake panels in particles (conversion = 0.0 Rev) was stable for slightly longer up to 34 revolutions. The other simulations converged before showing instability.

The simulation with panel to particle conversion after 4 revolutions was stopped after 32 rotations because of the computational cost. The computational cost increases as the number of wake panels increases because the $O\left(n^2\right)$ of the UVLM and also because the vortex panel and vortex particles interactions are not part of our current FMM implementation. The cost can be expressed as an order $O(n·m)$ which is reasonable if either the number of panels or the number of particle is small, but can quickly escalate and in the worst case reach an $O\left(n^2\right)$ cost if their number are even.

The global coefficients are averaged between revolutions 30 and 32 to compare as many stable simulations as possible, excluding conversion after 0.5 and 1.0 that never reached converged state.

Table 17. NL-UVLM-VPM standard deviation of global coefficients over rotations 30–32 for different numbers of wake panels rotations before being transformed in particles. The difference with the panel reference is also shown.

Conversion (Rev)	${\mathit{C}}_{\mathit{T}\mathit{std}}$ (%)	$\Delta {\mathit{C}}_{\mathit{T}\mathit{UVLM}}$ (%)	${\mathit{FM}}_{\mathit{std}}$ (%)	$\Delta {\mathit{FM}}_{\mathit{UVLM}}$ (%)
0	0.09	2.38	0.12	2.95
2	0.03	1.13	0.04	1.61
4	0.06	0.42	0.09	0.63
Panel (∞)	0.24	-	0.32	-

compares the overall coefficients deviation and difference with the panel reference. The expected trend is found, where the coefficients are converging toward the UVLM reference as the panel-particle conversion is performed later and the standard deviation is lower for the UVLM-VPM simulations than for the UVLM reference.

Figure 30 shows the mean spanwise distribution of thrust and torque. Only the conversion immediately behind the blades (conversion = 0.0) is significantly different than the panel reference. Both the conversions after 2.0 and 4.0 revolutions capture a similar blade loading as the panel reference, while still under-predicting slightly the peak loading, as expected from the $\Delta x_{tip}$ chosen in the previous section.

Figure 31 shows the mean position of the tip singularity for the first 3 revolutions (1080°). The conversion after 4.0 revolutions captures very well the tip singularity position of the panel reference for the first 3 rotations, while the conversion after 2.0 revolutions only captures well the first 1.5 rotations, as explained previously. Only the instant conversion in particles is significantly different in the first 1.5 rotations, explaining the observable difference on the spanwise loading.

Aside for the tip singularity position in the farther wake, the results between 2 and 4 wake panels revolutions are very similar to each other and close to the panel reference. Because of the increased computational cost with keeping more panel revolutions, it is decided to keep as little wake panels as needed to keep simulation stability and accurate results at higher time steps. It seems fortuitous that the optimal number of wake panels rotations to keep is the 2 revolutions initially imposed. This might be because of a bias in the selection of the other parameters, especially the Vreman model coefficient, caused by this initial hypothesis. Nonetheless, the similarity in the results using 2 and 4 revolutions and the panel reference is comforting as the simulation is eventually converging to a unique solution when enough wake panels rows are kept.

5.2.6. Summary of the UVLM-VPM Parametrization

The first parameter tested is the Vreman model coefficient of the PSE-LES, keeping 2 revolutions of wake panels and then creating the particles with a tip spacing of $\Delta x_{tip}=5°$. If the Vreman model coefficient is selected too small, the computed turbulent viscosity is too small to efficiently diffuse the large vorticity of the strongest particles and the simulations are numerically unstable. Increasing the Vreman model coefficient consistently increases the stability of the simulations, allowing them to reach more rotor revolutions before becoming numerically unstable. For the larger Vreman model coefficient values, the simulations are still stable after 70 revolutions without requiring any arbitrary particle destruction. The converged results appear to be mostly unaffected by the increase in Vreman model coefficient for a considerable interval of values. Increasing further the Vreman model coefficient affects the converged results more and more up to a point where they become chaotic. The Vreman model coefficient value is found to be different for each tip particle spacing regardless of the time step. For this specific test case, the appropriates values found in this work appear to be linear with tip spacing: $C_v=0.007$ for $\Delta x_{tip}=20°$, $C_v=0.014$ for $\Delta x_{tip}=10°$ and $C_v=0.028$ for $\Delta x_{tip}=5°$.

The second parameter tested is the particle tip spacing $\Delta x_{tip}$. It is found that $\Delta x_{tip}=20°$ is too coarse, $\Delta x_{tip}=10°$ is fine enough for global coefficient convergence and $\Delta x_{tip}=5°$ is required for accurate peak spanload. Increasing the number of particles decreases the differences with UVLM.

The third parameter investigated is the mesh refinement for the UVLM-VPM. It is clear that the UVLM-VPM is much less sensitive to mesh discretization than UVLM, where 16 × 40 appears to be sufficiently refined.

The fourth parameter that is varied is the time step. Unlike for the UVLM, the three tested time steps can be fully converged in a reasonable amount of computational time for the UVLM-VPM. The time step refinement consistently increases $C_T$ by about 1% and FM by about 1.4%.

The fifth parameter tested is the conversion time. It appears that reducing conversion time below 2 revolutions affects the simulation stability and increasing it reduces the difference with UVLM while also increasing the computational time. It is possible that choosing the Vreman model coefficient first introduced a bias towards conversion after 2 revolutions. Nonetheless, the spanload obtained with instantaneous conversion (0 revolutions) is clearly unable to capture the peak loading, whereas the conversion after 2 and 4 revolutions are very close to each other and not too far from the UVLM.

Overall, the addition of particles in the UVLM-VPM adds three parameters compared with the UVLM, but it also has four main advantages. Namely that the long wake panel convergence is eliminated (and thus requires less rotor revolution to obtain convergence), the computational cost scales linearly instead of quadratically with the number of element (so it is faster with a large number of element), the converged results show less deviation, and the mesh convergence is obtained at a coarser mesh than UVLM. The use of VPM is even more advantageous for other situations such as a hovering rotor in ground effect where UVLM alone becomes unstable.

6. Results

This section briefly describes the generation of the higher fidelity URANS 3D results along with the refinement study and then compares the results obtained with the present NL-UVLM-VPM method with the higher-fidelity URANS 3D.

6.1. URANS 3D

Appendix B presents the refinement study performed with the higher fidelity 3D URANS simulations that are used as the reference for the comparison with the NL-UVLM-VPM in Section 6.2. The higher fidelity 3D URANS results are generated with StarCCM+ version 13.04 (Siemens Digital Industries Software) using the $k-\omega$ SST turbulence model [45] with $\gamma$ transitional model. Results are generated for 5 different meshes at constant time step and for 5 different time steps on the medium mesh. In both cases, the global coefficients are compared and the spanload is compared for the mesh refinement. Preliminary experimental results suggest good correlations with the 3D URANS results for the global coefficients, but the detailed comparison with the experimental results is not yet available and will be part of a future work.

6.2. NL-UVLM-VPM and URANS 3D Comparison

Figure 32 compares the convergence of the global coefficients of the NL-UVLM-VPM with different databases presented in Section 4.2 with the convergence of the 3D URANS for the first 30 revolutions. The linear FM is not shown because it is so large that displaying it causes the rest of the results to look collapsed on each other.

Table 18. NL-UVLM-VPM standard deviation of global coefficients over rotations 28–30 for different viscous databases. The difference with the URANS averaged over rotations 40–50 reference is also shown.

Viscous Database	${\mathit{C}}_{\mathit{T}\mathit{std}}$ (%)	$\Delta {\mathit{C}}_{\mathit{T}\mathit{URANS}}$ (%)	${\mathit{FM}}_{\mathit{std}}$ (%)	$\Delta {\mathit{FM}}_{\mathit{URANS}}$ (%)
Linear	0.12	−9.29	0.11	622
GRT	0.05	1.99	0.06	15.3
KwSST	0.03	11.0	0.04	31.3
SA	0.04	5.05	0.06	8.85
SA Low Re	0.03	6.54	0.04	12.7

shows the standard deviation and the difference between the converged coefficients and the URANS reference. The results obtained with the 3D URANS reference are $2.060\times {10}^{-3}$ for $C_T$ and 0.093 for FM. The GRT database is the closest for the thrust prediction with a 2.0% over-prediction, followed by the SA and SA Low Re with their 5.0% and 6.5% over-prediction. The linear simulation is actually closer to the URANS thrust with an under-prediction of 9.3% than the KwSST with an over-prediction of 11.0%. Unfortunately, all of the NL-UVLM-VPM over-predict the thrust and under-predict the torque, leading to large differences in the figure of merit. Only the SA simulation is below 10% difference and SA Low Re and GRT databases are closer to 15% difference.

The sensitivity of the NL-UVLM-VPM to the database can be observed. In the worst case, the $C_T$ prediction has as much as a 9% difference between the GRT and the KwSST databases results. The figure of merit is even more affected by the choice of the database where comparing the KwSST and SA FM results in 23% difference for FM. However, the non linear coupling algorithm actually damps some of the differences in the inputted database lift coefficient. For instance, the lift slope of the KwSST is 25% higher than the slope of the GRT database near the tip and 56% higher near the root, leading to a difference in $C_T$ of only 9% between the NL-UVLM-VPM simulations using those database. Similarly, the SA Low Re slope is 13% higher than GRT near the tip and 22% higher near the root, leading to a $C_T$ difference of less than 5%. The link between $C_d$ of the database and the figure of merit is more complex to establish, but seems to be consistent with the database. At fairly similar $C_T$ and with similarly shaped $C_d$ curves, the SA and GRT torque coefficient and FM predictions are different by about 8% and 7%, respectively, while their $C_{d0}$ is about 13% different at the tip and about 3% at the root.

To understand the source of the large differences between the NL-UVLM-VPM and the URANS global coefficients results, the spanload is compared in Figure 33, Figure 34 and Figure 35. The differences in the spanload are greater near the root because the blade supports are modeled in the URANS 3D but not in the NL-UVLM-VPM. Indeed, the thrust is negative and the torque is much higher up to the end of the support at about 26% of the blade for the URANS reference compared with the NL-UVLM-VPM. The thrust coefficient distribution of the URANS reference reaches a peak very similar to the SA and SA Low Re NL-UVLM-VPM databases, albeit slightly more inboard. Because of the negative thrust near the root and the slow recovery of thrust in that vicinity, the integrated $C_T$ of the URANS reference is smaller than the GRT despite a higher and wider peak loading. The importance of the blade support on the overall coefficients is even more evident for the distributed torque coefficient. The URANS torque distribution is equal or below that of the SA and SA Low Re for the vast majority of the blade span except for the excess torque caused by the support near the root. Still, the URANS has a higher integrated $C_Q$ than all of the NL-UVLM-VPM simulations. One can easily imagine that the integrated torque would be between that of the SA databases and the StarCCM databases if the support region was ignored. However, the NL-UVLM-VPM spanloads are qualitatively different because they do not capture the dips in the torque coefficients happening at about 83% and 94%.

The twist about the half chord shown in Figure 35 does not show the SA Low Re curve, because it had an output problem, but it would be expected to be very close to the SA simulation since it has a very similar lift and moment polar as input. The linear result expectedly shows more twisting moment about the half chord, because just like the NL-UVLM-VPM simulations, it generates the lift at the quarter chord, but does not have a negative moment (nose down) interpolated from the database to reduce the resulting moment. Except for the KwSST, the NL-UVLM-VPM results predict a similar peak twisting moment, but wider and more outward than the URANS reference.

One might wonder what would happen to the overall coefficient comparison if the blade support region is excluded from the integration. It is not as easy as simply integrating from the end of the support to the end of the blade, because of the elliptic nature of the low speed aerodynamics. The mere presence of the support affects the airflow over the whole blade. Nonetheless, it makes sense to step a certain distance away from the support and assume the support influence becomes small. It is decided to do the integration 2 chords away from the support (that is from r/R = 47%) to the end of the blade. This reduced integration domain causes the URANS $C_T$ to be 4.7% lower and the figure of merit to be 5.0% higher.

Table 19. NL-UVLM-VPM global coefficients integrated from spanload, using r/R of 47% up to the tip. The difference with the URANS global coefficients also integrated over the same domain is shown.

${\int }_{\mathit{r}/\mathit{R}=47\%}^{\mathit{r}/\mathit{R}=100\%}\frac{\mathit{d}\mathit{Coef}}{\mathit{dr}}\mathit{dr}$
Viscous Database	$\Delta {\mathit{C}}_{\mathit{T}\mathit{URANS}}$(%)	$\Delta {\mathit{C}}_{\mathit{Q}\mathit{URANS}}$(%)	$\Delta {\mathit{FM}}_{\mathit{URANS}}$(%)
Linear	−14.3	−87.6	540
GRT	−3.26	−8.18	3.63
KwSST	5.32	−8.39	18.0
SA	−0.31	1.68	−2.11
SA Low Re	1.12	0.26	1.43

shows the comparison of the global coefficients computed from that truncated spanload data. It then becomes evident that the spanload of the NL-UVLM-VPM better matches the URANS spanload away from the support. The agreement of the coefficient of thrust is improved for the SA, SA Low Re and the KwSST databases by a significant amount. The GRT database now under-predicts the thrust, which is consistent with the observations in the peak loading area and has a similar error magnitude as before (from 2% to 3%). Only the linear solution really suffers from the reduced integration domain because the under-prediction it has in the peak loading area was partially canceled by the higher thrust in the support region. The difference with URANS is even more significantly reduced for the figure of merit for all of the NL-UVLM-VPM databases. The smallest difference with URANS results is obtained with SA Low Re, that has 1.1% and 1.4% for $C_T$ and FM, respectively. SA and GRT databases also give satisfactory agreement with URANS, but KwSST has significantly more error with 5.3% over-prediction of $C_T$ and 18% over-prediction of FM. Linear UVLM-VPM is of course the worst performing simulation on this lower Reynolds rotor, where viscous effects are important to consider.

The main takeaways of this sections are that the aerodynamic prediction of the NL-UVLM-VPM can agree well with that of a 3D URANS simulation on a lower Reynolds rotor blade in regions where the geometric modeling is the same but the results greatly depend on the quality of the database.

7. Conclusions

In conclusion, a comparative literature review presented the different methods and validation of some of the most recent published potential methods for solving the aerodynamics of an isolated rotor in hover. The literature comparison highlighted the difficulty of comparing the methods one-to-one because of the differences in the results shown and the choice of the validation cases. It also discusses the current unavailability of a complete validation cases for potential methods on smaller rotors such as those of small UAVs. Then, the present NL-UVLM-VPM potential method was briefly summarized, including three improvements from the previously published method, namely the vortex particles initial core size, their initial spacing, and their creation time (keeping more wake panels rows). Next, a simple low Reynolds and Mach 2-bladed hovering rotor test case was defined and used for a detailed parametric and refinement study to evaluate the sensitivity of all the inputs of the current method. The sensitivity study first investigated the wake of panels only (UVLM) and then introduced the wake vortex particles (UVLM-VPM). The UVLM parametrization has found that the method was consistent with the expected behavior of both the temporal and spatial refinements and with respect to the core size selection. It was determined that a blade mesh of at least 24 chordwise by 60 spanwise panels, a time step of 5° and a core size of 60% of the blade chord was needed. For the UVLM-VPM parametrization, it was first determined that the Vremann coefficient model had a large interval of values that could stabilize the simulation while leading to similar results. That parameter was found to vary with particle spacing, but not with time step. Then, it was found that a tip particle spacing of 10° is sufficient to converge on the global coefficient, but 5° was needed to converge on the distributed loads. The mesh refinement revealed that the NL-UVLM-VPM simulations were less sensitive to mesh coarsening than UVLM, showing convergence with a mesh of 16 chordwise by 40 spanwise panels. The time step refinement was easier to perform in a reasonable computational time with the use of vortex particles, because the increased computation cost caused by smaller time steps was accelerated by the use of Fast Multipole Method (FMM). Finally, the conversion time was analyzed, showing some advantages to keep between 2 and 4 revolutions of panels before converting them into particles. The use of UVLM-VPM introduces three more parameters compared to UVLM, but also have four advantages: (1) reduced number of revolutions to simulate to obtain convergence, (2) decreased computational cost with high number of elements, (3) less variance in the results and (4) smaller sensitivity to mesh coarsening. The article finally compares the results obtained with the higher fidelity 3D URANS and those obtained with the present method using different databases for the non-linear strip wise coupling. The results show good agreement between the higher fidelity method and the present NL-UVLM-VPM method aside from the blade support region that is not modeled in the present method. From all the detailed studies above, the key results are that the NL-UVLM-VPM is consistent with discretization refinement both in time and space, is stable for long time running (80 revolutions and more) with the appropriate Vreman model coefficient for the PSE-LES, and provides accurate aerodynamic loading at a fraction of the computational cost when compared with 3D URANS, even at lower Reynolds, where viscous effects are important. These results constitute important steps towards providing new generation tools for early aerodynamic design of rotorcraft. Future work should explore more flight conditions (e.g., forward or descending flight) and more realistic geometrical configurations (e.g., fuselage, rotor hub). The method could be used for multidisciplinary design optimization or for aeroelasticity applications.