An Experimental-Numerical Investigation of the Wake Structure of a Hovering Rotor by PIV Combined with a Γ₂ Vortex Detection Criterion

Abstract

The rotor wake aerodynamic characterization is a fundamental aspect for the development and optimization of future rotary-wing aircraft. The paper is aimed at experimentally and numerically characterizing the blade tip vortices of a small-scale four-bladed isolated rotor in hover conditions. The investigation of the vortex decay process during the downstream convection of the wake is addressed. Two-component PIV measurements were carried out below the rotor disk down to a distance of one rotor radius. The numerical simulations were aimed at assessing the modelling capabilities and the accuracy of a free-wake Boundary Element Methodology (BEM). The experimental and numerical results were investigated by the Γ₂ criterion to detect the vortex location. The rotor wake mean velocity field and the instantaneous vortex characteristics were investigated. The experimental/numerical comparisons show a reasonable agreement in the estimation of the mean velocity inside the rotor wake, whereas the BEM predictions underestimate the diffusion effects. The numerical simulations provide a clear picture of the filament vortex trajectory interested in complex interactions starting at about a distance of z/R = −0.5. The time evolution of the tip vortices was investigated in terms of net circulation and swirl velocity. The PIV tip vortex characteristics show a linear mild decay up to the region interested by vortex pairing and coalescence, where a sudden decrease, characterised by a large data scattering, occurs. The numerical modelling predicts a hyperbolic decay of the swirl velocity down to z/R = −0.4 followed by an almost constant decay. Instead, the calculated net circulation shows a gradual decrease throughout the whole wake development. The comparisons show discrepancies in the region immediately downstream the rotor disk but significant similarities beyond z/R = −0.5.

1. Introduction

During the generation of the required thrust, the helicopter rotor blades produce a complex wake system which is characterized by spanwise shed vortices and trailing vortices having a strength varying along the blade span. In particular, a helical system of strong blade tip vortices develops because of the rotor revolution This vortex system can interact with the main rotor, the tail rotor, and the airframe. Hence, vortex formation and development are important factors influencing the aerodynamics of the entire helicopter. In undisturbed hover, the blade pitch setting does not vary with time and constant lift and induced velocity distributions are therefore generated during the rotor revolution. This leads to a relatively simple vortex system characterized by vortices keeping a constant strength over the azimuth. The characteristics and development of the tip vortices are aspects of fundamental importance for the understanding of the rotor wake evolution.

Over the past twenty years, Particle Image Velocimetry (PIV) has made the detailed investigation of the properties of the blade tip vortices possible by providing planar velocity data. Many PIV measurements on hovering helicopter rotors have been performed. Martin et al. [1] applied phase-resolved stereoscopic PIV measuring on a sub-scale two-bladed helicopter rotor in hovering up to one revolution age. The results were compared with three-component Laser Doppler Velocimetry (LDV) measurements using the same seeding medium. Information on the required spatial resolution to resolve the tip vortices was provided. In the same year, Heineck et al. [2] used a stereo PIV set-up on a two-bladed rotor test stand to capture blade tip vortices up to an age of 270 deg, focusing on the data processing and the effect of the vortex wandering on the vortex core size. McAlister [3] investigated the rotor wake of a two-bladed model helicopter up to a wake age of 390 deg. Conclusions on the temporal development of the maximum swirl velocity as well as the core growth were drawn using a stereo PIV setup. Richard and van der Wall [4] published an analysis of the two-component and stereo PIV on a four-bladed rotor in hover condition up to the half revolution, which encompassed a three-dimensional reconstruction of a vortex filament based on the PIV data and indication about the necessary spatial resolution to resolve properly the vortex characteristics. Many studies have seen the use of PIV for the investigation of the rotor aerodynamics or the full helicopter configuration [5] and a complete overview is given by the paper of Raffel et al. [6].

Several methodologies have been developed and applied during the years to the numerical investigation of the helicopter aerodynamics and the related wake structure, each with their level of sophistication and limitations. A historical review of these methodologies can be found in Leishman [7] and Johnson [8]. Two examples can be mentioned to highlight the pros and cons in the application of these methodologies: considering the sophisticated CFD-based tools, they solve the governing fluid dynamic equations in a region of space surrounding the configuration to be analysed with a level of accuracy that is dependent on the resolution of the volume grid applied for the numerical discretization, and this requires to be extremely fine around the blades and/or the airframe surfaces and in the wake downstream, where the viscous stresses are the highest, to avoid an excessive non-physical numerical dissipation. For this reason, these methods are usually computationally onerous. In addition, the accuracy of the solution also depends on the turbulence model applied. Conversely, considering the lower-fidelity free-wake panel methods, they require limited computational resources because are based on the simplifying assumptions of inviscid and incompressible flow and need just a surface discretization of the wake and body geometries. They are free from numerical dissipation but need the use of suitable regularization models to avoid unphysical singularities produced by application of the Biot–Savart law for the evaluation of the velocities induced by the wake vortex filaments and appropriate models to model the wake dissipation process caused by the natural flow viscosity.

Regardless of the applied methodology, it appears that the accurate numerical reproduction of blade tip vortices is challenging and requires validation [9]. Duraisamy et al. [10] published a direct comparison of a computational fluid dynamics (CFD) simulation and respective experiments focusing on the physics of vortex formation from a single-bladed hovering rotor. They found PIV to be a valid method for a qualitative and quantitative comparison. Unfortunately, the amount of experimental data concerning the vortex formation, development, and decay for the validation of the numerical approaches in rotating systems mostly do not investigate beyond one rotor revolution. In the past, Caradonna [11] investigated the tip trajectory up to a wake age of 1080 deg. by flow visualization and these data were used for validation [12,13]. More recently, the large interest in the brownout phenomenon drove many investigations to characterise the rotor wake. Lee [14] presented remarkable flow visualization results of a hovering rotor out of ground effect but without providing quantitative data of the blade tip vortices evolution.

The current work stems from research activity in the framework of the GARTEUR Action Group 22 Forces on Obstacles in Rotor Wake (Visingardi et al. [15]) to evaluate the mutual effects numerically and experimentally between a small-scale helicopter rotor in hover flight and a cylindrical sling load located at different positions below the rotor disk. The experimental investigations focused on the development of the rotor wake far from the rotor disk. The numerical simulations were performed with the main purpose to assess the modelling capabilities and the accuracy of a free-wake Boundary Element Methodology. The results were presented in a paper of the 44th European Rotorcraft Forum discussing the effect of the wake on the sling load [16].

Following this activity, new investigations focused on the blade tip vortex characteristics and the dissipation phenomena when moving away from the disk of the isolated rotor (i.e., without cylinder) up to a wake age of three rotor revolutions. Particular attention was dedicated to the vortex detection criterion adopted on the PIV data. The most widely used local methods for vortex detection are founded on the velocity gradient tensor ∇u and its three invariants. Examples are the ∆ criterion introduced by Dallmann [17], Vollmers et al. [18], and Chong et al. [19]; the Q criterion by Hunt et al. [20]; λ₂ criterion by Jeong and Hussain [21]. These local vortex-detection criteria are not always suitable for noisy PIV data affected by spurious vectors, thus resulting in high-velocity gradients. The Γ₂ criterion proposed by Graftieaux et al. [22] offered a possible solution and later it was successfully applied to complex wind tunnel measurements by Mulleners and Raffel [23].

The current work illustrates the results of these new numerical and experimental investigations. The paper is organized into sections. Section 2 describes the main characteristics of the four-bladed rotor rig, the PIV system including the data evaluation procedure and a description of the adopted vortex identification criterion. The numerical methodology is illustrated in Section 3, whereas the numerical/experimental comparison of the results obtained is fully documented in Section 4. The conclusions are finally reported in Section 5.

2. Experimental Setup and Test Conditions

A dedicated rotor test rig was built based on an existing commercial radio-controlled helicopter model (Blade 450 3D RTF), (Figure 1a). The original two-bladed rotor was replaced by a four-bladed one with collective and cyclic control. The rotor blades were untwisted, with rectangular planform and parabolic tip. The blades presented a radius of $\mathrm{R}=0.36\mathrm{m}$ and a chord length of $\mathrm{c}=0.0327\mathrm{m}$. The root cut-out was located at 16% of the radius. A NACA0013 airfoil was used throughout the blade span. The blade planform and the geometry of the airfoil are shown in Figure 1b. The rotor solidity was equal to $\mathsf{\sigma }=\left({\mathrm{N}}_{\mathrm{b}}\mathrm{c}\right)/\mathsf{\pi }\mathrm{R}=0.116$ and the aspect ratio of the blades was $\mathrm{AR}=\mathrm{R}/\mathrm{c}\approx 11$. The rotor spun clockwise when seen from above, the collective pitch angle ${\mathsf{\theta }}_0$ varied from $1°$ to $12.2°$ and the maximum speed was $\mathsf{\Omega }=1780\mathrm{RPM}$.

The force and moments produced by the rotor rig were measured by a six components balance (ATI MINI40). The detailed characteristics in terms of full scale and accuracy are summarized in

Table 1. Balance characteristics.

	Fx (N)	Fy (N)	Fz (N)	Mx (Nm)	My (Nm)	Mz (Nm)
Full scale	±20	±20	±60	±1	±1	±1
Accuracy (%FS)	0.25	0.25	0.60	0.0125	0.0125	0.0125

.

The rotor test rig was located at a distance of 5 R from the floor and 3 R from the ceiling to avoid any influence of the surrounding walls. One hall-effect sensor was located on the shaft gear for measuring the rotating speed and for providing a trigger TTL signal at a prefixed azimuth angle to allow phase-locked measurements.

The current investigation addressed the wake downwash generated by the four-bladed rotor in hover conditions at a constant collective angle of ${\mathsf{\theta }}_0=11.8°$. The angular velocity was set to $\mathsf{\Omega }=1740\mathrm{RPM}$, which leads to a blade tip velocity of V_tip = 66 m/s and a thrust value of T = 12 N. The resulting blade loading was ${\mathrm{C}}_{\mathrm{T}}/\mathsf{\sigma }=\mathrm{T}/\left({\mathsf{\rho }\mathrm{A}\mathsf{\Omega }}^2{\mathrm{R}}^2\mathsf{\sigma }\right)=0.052$ with the density $\mathsf{\rho }=1.114\mathrm{kg}/{\mathrm{m}}^3$, and the rotor area A = 0.41 m². The Mach and Reynolds numbers are given at the radius tip (

Table 2. Parameter of the measured test case.

Ω/2π (Hz)	Vtip (m/s)	Mtip	Retip
29	65.6	0.19	1.47*105

).

2.1. PIV System and Evaluation

A fixed frame of reference was defined and having the origin located in the rotor centre with the x-axis horizontally oriented along the rotor blade, the y-axis orthogonal to the x-axis and lying in the rotor disk plane, and the z-axis vertically and upward directed. The rotor wake characteristics were investigated by a standard two-component PIV measurement system composed of a double head Nd-Yag laser with a maximum energy of 300 mJ per pulse at 532 nm and a single double frame CCD camera (2048 by 2048 px) with a dynamic range of 14 bits. The camera was mounted on two components translating system to cover the full region of interest in the xz-plane. The light sheet was vertically oriented and aligned with the rotor blade at the azimuth angle of Ψ = 180°. A standard PIV layout was adopted with the camera line of sight orthogonal to the laser sheet (Figure 2).

The wake downwash was investigated up to 2.25 rotor radii from the rotor disk. The camera, equipped with a 50 mm lens, recorded a measurement region of 320 mm by 320 mm, with a spatial resolution of 6 pixel/mm in the image plane. Separate measurement regions needed to cover the full wake (Figure 3a). The data post processing discussed in [16] provided a spatial resolution of $\Delta x=2.5\mathrm{mm}$ not sufficient for the characterization of the tip vortices.

To track the blade tip vortices in proximity to the rotor disk, the camera was equipped with a lens featuring a fixed focal length of f = 200 mm, and the f-number was set to f_# = 2.8. The measurement region was located immediately below the rotor disk plane on a vertical plane radially ranging between x/R = 0.68 and x/R = 1.08 to identify the trajectories of the trailing tip vortices. Five PIV measurement regions with the size of about 120 mm by 120 mm, partially overlapped, were used to measure the rotor wake, and in particular the blade tip vortices characteristics up to one radius downstream the rotor disc, Figure 3b. This yielded a spatial resolution of about 17 px/mm in the image plane. The delay time between the two laser pulses was set to 25 μs according to the highest expected velocities in the flow. As tracer particles, sprayed diethylhexylsebacate (DEHS) oil was used. A seeding generator with 20 Laskin nozzles provided oil droplets with an average size less than 1 μm. The full test room was seeded to have a homogenous concentration of particles. More than 150 image pairs were recorded at a frame rate of 3 Hz for each region of interest (ROI) over about 1450 rotor revolutions. The PIV images were pre-processed by applying a background grey-level subtraction. PIV-View 3.60 was used to process the images. The analysis consisted in a Multi-grid scheme with a B-Spline of 3rd order image deformation ending at 32 × 32 px² and 75% overlap. Correlation maps were calculated by FFT multiple correlations (Hart [25]) of 2 windows and a 3-point Gaussian peak fit was used to obtain the displacement.

The results presented a velocity spatial resolution of $\Delta x=0.47\mathrm{mm}$. The random noise of the PIV cross-correlation procedure can be estimated as 0.1 px as a rule–of–thumb (Raffel et al. [26]). Using the current values for the optical resolution (17 px/mm) and the laser double–pulse delay (25 μs), this leads to a velocity error of $\Delta V$ of $\approx 0.23\mathrm{m}/\mathrm{s}$ for the PIV measurements. In the proximity to the rotor disk, the core radius r_c of the tip vortices was measured. The core radius is defined as the distance from the vortex centre to the radial position where the maximum swirl velocity is reached (Figure 4). Values of r_c between 3 and 3.3 mm were measured giving a ratio $\Delta x/r_c$ of about 0.15–0.14 according to the value of $\Delta \mathrm{x}/{\mathrm{r}}_{\mathrm{c}}\le 0.2$ indicated by Martin et al. [1] to guarantee a correct vortex characterization.

2.2. Vortex-Identification Criterion

To overcome the difficulty in identifying the centre of the vortices in the presence of spurious vectors, a method based on the velocity field topology, without using velocity derivatives was chosen. This was the ${\mathsf{\Gamma }}_2$-criterion proposed by Graftieaux et al. [22]. The function ${\mathsf{\Gamma }}_2$ is defined in discrete form as:

$${\mathsf{\Gamma }}_2\left(\overrightarrow{{\mathrm{x}}_{\mathrm{i}}}\right)=\frac{1}{\mathrm{M}}{\displaystyle \sum }_{{\mathrm{x}}_{\mathrm{j}}\in {\mathrm{S}}_{\mathrm{i}}}\frac{\left\{\left(\overrightarrow{{\mathrm{x}}_{\mathrm{J}}}-\overrightarrow{{\mathrm{x}}_{\mathrm{i}}}\right)\times \left(\overrightarrow{{\mathrm{u}}_{\mathrm{j}}}-\overrightarrow{{\mathrm{u}}_{\mathrm{mean}}}\right)\right\}·\overrightarrow{\mathrm{n}}}{\left|\overrightarrow{{\mathrm{x}}_{\mathrm{J}}}-\overrightarrow{{\mathrm{x}}_{\mathrm{i}}}\right|\left|\left(\overrightarrow{{\mathrm{u}}_{\mathrm{j}}}-\overrightarrow{{\mathrm{u}}_{\mathrm{mean}}}\right)\right|}$$

(1)

with ${\mathrm{S}}_{\mathrm{i}}$ representing a two-dimensional circle around $x_{\mathrm{i}}$ with radius D, $\mathrm{M}$ the number of grid points ${\mathrm{x}}_{\mathrm{j}}$ inside ${\mathrm{S}}_{\mathrm{i}}$ with $\mathrm{j}\ne \mathrm{i},$ u_mean the average velocity vector within the region S, $\overrightarrow{\mathrm{n}}$ the unit vector normal to the PIV plane and ${\mathrm{u}}_{\mathrm{j}}$ the velocity at ${\mathrm{x}}_{\mathrm{j}}$. The l radius $\mathrm{D}$ of the domain ${\mathrm{S}}_{\mathrm{i}}$ is expressed in terms of grid spacing. ${\mathsf{\Gamma }}_2$ is a 3D non-dimensional scalar function, with $-1\le {\mathsf{\Gamma }}_2\le 1$. The zones delimited by $\left|{\mathsf{\Gamma }}_2\right|>\frac{2}{\mathsf{\pi }}$ identify the vortices present in the measurement region. The vortex centre is identified as the maximum value of the absolute of ${\mathsf{\Gamma }}_2$ in the delimited zone. In cases of PIV data affected by spurious vectors or data void, it is suggested to use the weighted centroid of the ${\mathsf{\Gamma }}_2$ for the identification of the centre of the vortex. In the current work, the weighted centroid is chosen to identify the centres. The choice of the domain radius $\mathrm{D}$ influences the accuracy of the centre detection and the dimension of the identified vortices. De Gregorio and Visingardi [27] indicated that the selection of a domain radius $\mathrm{D}$, equal to the larger core radius existing within the investigated flow field, provides the best measurement accuracy for single and multiple vortices affected by spurious vectors. In the case of elliptical vortices, a domain radius $\mathrm{D}$ equal to the semi-major axis is recommended for reducing centre detection errors, whereas the magnitude of the domain radius is recommended to be set between two to six times the core radius size in the case of significant data void in the vortex core.

In the processing of the PIV images, the lack of particles caused a number of spurious velocity vectors. A normalised median test, presented by Westerweel and Scarano [28], coupled to a Dynamic Mean Test and a Global Histogram Filter (described by Raffel et al. [26]) were used to detect and then remove the spurious vectors. After deleting the outliers, the missing data were replaced using bi-linear interpolation. Taking into account the particle void, the domain radius D was fixed equal to 10 grid spacing, which was almost twice the value of the vortex core radius ($r_c\approx 6$ grid point as shown in Figure 5). Once the centre was detected, the vortex characteristics were calculated on concentric circles moving away from the centre. The raw or filtered data have shown scattering data in the core, while the interpolated data provided a satisfactory agreement with the Vatistas [29] theoretical curve (Figure 5). The agreement between the linear interpolated data and the Vatistas vortex is explained by the nature of a real vortex core where the tangential velocity behaviour is linear and can be expressed as $V_{\theta }=\Omega ·r$, thus justifying the application of a bi-linear interpolation to account for the missing vectors.

3. Numerical Methodology

The numerical simulations were carried out by using the code RAMSYS [30], which is an unsteady, inviscid and incompressible free-wake vortex lattice Boundary Element Method (BEM) solver for multi-rotor, multi-body configurations developed at CIRA. It is based on Morino’s boundary integral formulation for the solution of the Laplace equation for the velocity potential $\mathsf{\phi }$ [31].

3.1. Governing Equation

The fluid is assumed to be inviscid, incompressible and irrotational. Under these hypotheses there exists a function $\mathsf{\phi }\left(x,\mathrm{t}\right)$, called velocity potential, such that the velocity vector $v=\nabla \mathsf{\phi }$ and the continuity equation reduces to Laplace equation:

$${\nabla }^2\mathsf{\phi }\left(\mathrm{x},\mathrm{t}\right)=0$$

(2)

The application of the Green’s function to Equation (2) yields a boundary integral representation of the velocity potential:

$$\mathsf{\phi }\left(x,\mathrm{t}\right)={{\displaystyle \int }}_{{\mathrm{S}}_{\mathrm{B}}}^{}\left(\mathrm{G}\frac{\partial \mathsf{\phi }}{\partial \mathrm{n}}-\mathsf{\phi }\frac{\partial \mathrm{G}}{\partial \mathrm{n}}\right)\mathrm{dS}\left(y\right)-{{\displaystyle \int }}_{{\mathrm{S}}_{\mathrm{W}}}^{}\mathsf{\Delta }\mathsf{\phi }\frac{\partial \mathrm{G}}{\partial \mathrm{n}}\mathrm{dS}\left(y\right)$$

(3)

where ${\mathrm{S}}_{\mathrm{B}}$ and ${\mathrm{S}}_{\mathrm{W}}$ are body and wake surfaces, respectively, and $\mathrm{G}=-1/4\mathsf{\pi }\left|\right|y-x\left|\right|$ is the free-space fundamental solution of the 3D Laplace equation.

3.2. Boundary Conditions

Far from the body, the velocity potential is null at infinity:

$$\underset{x\to \infty }{\lim }\mathsf{\phi }\left(x,\mathrm{t}\right)=0$$

(4)

The impermeability boundary condition on ${\mathrm{S}}_{\mathrm{B}}$ yields:

$$\frac{\partial \mathsf{\phi }}{\partial \mathrm{n}}={\mathrm{V}}_{\mathrm{B}}·n$$

(5)

where ${\mathrm{V}}_{\mathrm{B}}$ denotes the velocity of body points and $n$ denotes the outward unit normal vector on ${\mathrm{S}}_{\mathrm{B}}$.

The boundary condition on the wake is expressed as:

$$\frac{{\mathrm{D}}_{\mathrm{W}}\Delta \mathsf{\phi }}{\mathrm{Dt}}=0$$

(6)

according to which the potential jump $\Delta \mathsf{\phi }$ across the wake surface is constant:

$$\Delta \mathsf{\phi }\left(x_W,\mathrm{t}\right)=\Delta {\mathsf{\phi }}^{\mathrm{TE}}\left(\mathrm{t}-\tau \right)$$

(7)

with $\tau$ denoting the time taken by the wake material point $x_W$ to move from the trailing edge to its current position.

3.3. Novel Boundary Integral Formulation

The novel formulation proposed by Gennaretti et al. [32] is applied in RAMSYS to avoid the instabilities arising in the numerical formulation when wake panels are too close to or impinge the body. In this formulation, the velocity potential $\mathsf{\phi }$ is split into a scattered potential ${\mathsf{\phi }}_{\mathrm{s}}$, generated by sources and doublets over ${\mathrm{S}}_{\mathrm{B}}$ and doublets over the part of the wake surface in contact with the blade trailing edge (Near wake, ${\mathrm{S}}_{\mathrm{W}}^{\mathrm{N}}$), and an incident potential ${\mathsf{\phi }}_{\mathrm{I}}$, generated by doublets over the part of the wake not in contact with the trailing edge (Far wake, ${\mathrm{S}}_{\mathrm{W}}^{\mathrm{F}}$), such that ${\mathrm{S}}_{\mathrm{W}}={\mathrm{S}}_{\mathrm{W}}^{\mathrm{N}}{\displaystyle \cup }{\mathrm{S}}_{\mathrm{W}}^{\mathrm{F}}$ and $\mathsf{\phi }={\mathsf{\phi }}_{\mathrm{s}}+{\mathsf{\phi }}_{\mathrm{I}}$.

The application of this novel formulation provides a new expression of Equation (3), which is then replaced by:

$$\begin{array}{cc}{\mathsf{\phi }}_{\mathrm{s}}\left(x,\mathrm{t}\right)=& {{\displaystyle \int }}_{{\mathrm{S}}_{\mathrm{B}}}^{}\left[\mathrm{G}\left(\chi -{\chi }_{\mathrm{I}}\right)-{\mathsf{\phi }}_{\mathrm{s}}\frac{\partial \mathrm{G}}{\partial \mathrm{n}}\right]\mathrm{dS}\left(y\right)\\ & -{{\displaystyle \int }}_{{\mathrm{S}}_{\mathrm{W}}^{\mathrm{N}}}^{}\Delta {\mathsf{\phi }}_{\mathrm{s}}{}^{\mathrm{TE}}\left(\mathrm{t}-\tau \right)\frac{\partial \mathrm{G}}{\partial \mathrm{n}}\mathrm{dS}\left(y\right)\end{array}$$

(8)

and by:

$${\mathsf{\phi }}_{\mathrm{I}}\left(x,\mathrm{t}\right)=-{{\displaystyle \int }}_{{\mathrm{S}}_{\mathrm{W}}^{\mathrm{F}}}^{}{\mathsf{\Delta }\mathsf{\phi }}_{\mathrm{I}}\frac{\partial \mathrm{G}}{\partial \mathrm{n}}\mathrm{dS}\left(y\right)$$

(9)

where the boundary conditions are given by:

$$\frac{\partial {\mathsf{\phi }}_{\mathrm{s}}}{\partial \mathrm{n}}=V_{\mathrm{B}}·n-v_{\mathrm{I}}·n=\chi -{\chi }_{\mathrm{I}}$$

(10)

and:

$$\Delta {\mathsf{\phi }}_{\mathrm{I}}\left(x_W,\mathrm{t}\right)=\Delta {\mathsf{\phi }}_{\mathrm{s}}{}^{\mathrm{TE}}\left(\mathrm{t}-\tau \right)$$

(11)

with the velocity induced by the far wake obtained from the gradient of Equation (9) combined with Equation (11),

$$v_{\mathrm{I}}\left(x,\mathrm{t}\right)=-{\nabla }_{\mathrm{x}}{{\displaystyle \int }}_{{\mathrm{S}}_{\mathrm{W}}^{\mathrm{F}}}^{}\Delta {\mathsf{\phi }}_{\mathrm{s}}{}^{\mathrm{TE}}\left(\mathrm{t}-\tau \right)\frac{\partial \mathrm{G}}{\partial \mathrm{n}}\mathrm{dS}\left(y\right)$$

(12)

3.4. Numerical Solution

The numerical solution of Equation (8) is obtained by defining boundary elements, i.e., by discretizing ${\mathrm{S}}_{\mathrm{B}}$ and ${\mathrm{S}}_{\mathrm{W}}$ into quadrilateral panels, assuming ${\mathsf{\phi }}_{\mathrm{s}}$, $\partial {\mathsf{\phi }}_{\mathrm{s}}/\partial \mathrm{n}$, $\Delta {\mathsf{\phi }}_{\mathrm{s}}$ to be piecewise constant, and imposing that the equation is satisfied at the centre of each body element (collocation method). Specifically, dividing the blade surface into M panels, ${\mathrm{S}}_{{\mathrm{B}}_{\mathrm{i}}}$, and the wake surface into N panels, ${\mathrm{S}}_{{\mathrm{W}}_{\mathrm{j}}}$, the discretized version of Equation (8) gives the linear algebraic system:

$${\mathsf{\phi }}_{\mathrm{s}}{}_{\mathrm{k}}={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{M}}{\mathrm{B}}_{\mathrm{ki}}\left[\chi \left(\mathrm{t}\right)-{\chi }_{\mathrm{I}}\left(\mathrm{t}\right)\right]+{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{M}}{\mathrm{C}}_{\mathrm{ki}}{\mathsf{\phi }}_{\mathrm{s}}{}_{\mathrm{i}}\left(\mathrm{t}\right)+{\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{N}}{\mathrm{F}}_{\mathrm{kj}}{\mathsf{\Delta }\mathsf{\phi }}_{\mathrm{s}}{}_{\mathrm{j}}\left(\mathrm{t}\right)$$

(13)

where ${\mathsf{\phi }}_{\mathrm{s}}{}_{\mathrm{i}}\left(\mathrm{t}\right)={\mathsf{\phi }}_{\mathrm{s}}\left(x_{\mathrm{i}},\mathrm{t}\right)$, ${\chi }_{\mathrm{i}}\left(\mathrm{t}\right)=\chi \left(x_{\mathrm{i}},\mathrm{t}\right)$, ${\chi }_{\mathrm{I}}{}_{\mathrm{i}}\left(\mathrm{t}\right)={\chi }_{\mathrm{I}}\left(x_{\mathrm{i}},\mathrm{t}\right),{\mathsf{\Delta }\mathsf{\phi }}_{\mathrm{s}}{}_{\mathrm{j}}\left(\mathrm{t}\right)=\Delta {\mathsf{\phi }}_{\mathrm{s}}{{}_{\mathrm{j}}}^{\mathrm{TE}}\left(\mathrm{t}\right)$ whereas source/sink and doublet coefficients are given by:

$$\begin{gathered} {\mathrm{B}}_{\mathrm{ki}}=-\frac{1}{4\mathsf{\pi }}{{\displaystyle \int }}_{{\mathrm{S}}_{{\mathrm{B}}_{\mathrm{i}}}}^{}\left(\frac{1}{\left|y-x_{\mathrm{k}}\right|}\right)\mathrm{dS} \\ {\mathrm{C}}_{\mathrm{ki}}=-\frac{1}{4\mathsf{\pi }}{{\displaystyle \int }}_{{\mathrm{S}}_{{\mathrm{B}}_{\mathrm{i}}}}^{}\frac{\partial }{\partial \mathrm{n}}\left(\frac{1}{\left|y-x_{\mathrm{k}}\right|}\right)\mathrm{dS} \\ {\mathrm{F}}_{\mathrm{kj}}=-\frac{1}{4\mathsf{\pi }}{{\displaystyle \int }}_{{\mathrm{S}}_{{\mathrm{W}}_{\mathrm{j}}}}^{}\frac{\partial }{\partial \mathrm{n}}\left(\frac{1}{\left|y-x_{\mathrm{k}}\right|}\right)\mathrm{dS} \end{gathered}$$

(14)

The solution of the algebraic system is obtained by the application of the GMRES iterative method.

Once the potential field is known, the surface pressure distributions are evaluated by applying the unsteady version of the Bernoulli equation:

$$\frac{\partial {\mathsf{\phi }}_{\mathrm{s}}}{\partial \mathrm{t}}+\frac{\partial {\mathsf{\phi }}_{\mathrm{I}}}{\partial \mathrm{t}}-V_{\mathrm{B}}·\left(\nabla {\mathsf{\phi }}_{\mathrm{s}}+v_{\mathrm{I}}\right)+\frac{\left|\right|\nabla {\mathsf{\phi }}_{\mathrm{s}}+v_{\mathrm{I}}{}^2\left|\right|}{2}+\frac{\mathrm{p}}{\mathsf{\rho }}=\frac{{\mathrm{p}}_{\infty }}{\mathsf{\rho }}$$

(15)

where the incident potential ${\mathsf{\phi }}_{\mathrm{I}}$ is obtained from integration of the incident velocity field by using Equation (12) with the inclusion of the vortex-core model. Finally, the evaluation of the forces and moments is obtained by the integration of Equation (15).

3.5. Vortex Core Model

To account for the viscous diffusion of the wake vortex elements, the Vatistas vortex core model was used, according to which the swirl velocity is expressed as:

$${\mathrm{V}}_{\mathsf{\theta }}=\frac{{\mathrm{r}\mathsf{\Gamma }}_{\mathrm{v}}}{2\mathsf{\pi }{\left({\mathrm{r}}^{2\mathrm{n}}+{\mathrm{r}}_{\mathrm{c}}{}^{2\mathrm{n}}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\mathrm{n}$}\right.}}$$

(16)

where the coefficient $n$ has been set equal to “1”, as suggested by Scully [33].

The applied diffusion model is the one described by Squire [34]. In this model, the growth with the time of the core radius ${\mathrm{r}}_{\mathrm{c}}$ is given by:

$${\mathrm{r}}_{\mathrm{c}}=\sqrt{{\mathrm{r}}_{\mathrm{c}0}^2+4\alpha \delta \nu \mathrm{t}}$$

(17)

where the term ${\mathrm{r}}_{\mathrm{c}0}$ is the initial core radius that removes the singularity at ${\mathrm{t}}_0$, and was set equal to 5% of the blade average chord length c in the calculations, the term α is the Oseen coefficient and is equal to 1.25643. The product δv represents the “eddy viscosity” where v is the kinematic viscosity and:

$$\delta =1+{\mathrm{a}}_1\frac{{\mathsf{\Gamma }}_{\mathrm{v}}}{\nu }$$

(18)

represents an average effective (turbulent) viscosity coefficient in which ${\mathsf{\Gamma }}_{\mathrm{v}}$ is the circulation strength of the vortex element, while the Squire’s coefficient ${\mathrm{a}}_1$ is an empirical parameter specified to vary between 0.2 and 0.0002, as indicated in Bhagwat [35]. For a small-scale rotor, like the one used for these investigations, a value of $\mathrm{O}\left({10}^{-4}\right)$ can be used. The model suggested by Donaldson & Bilanin [36] was used to take into account the decay of the circulation ${\mathsf{\Gamma }}_{\mathrm{v}}$ with time. According to this model, the circulation of the tip vortex ${\mathsf{\Gamma }}_{\mathrm{v}}\left(\mathrm{t}\right)$ is expressed as:

$${\mathsf{\Gamma }}_{\mathrm{v}}\left(\mathrm{t}\right)={\mathsf{\Gamma }}_0\exp \left(-\frac{\mathrm{bq}}{\mathrm{s}}\mathrm{t}\right)$$

(19)

being $\mathrm{b}$ a decay coefficient; $\mathrm{q}$ the ambient turbulence level and $\mathrm{s}$ the aircraft semispan. In the present calculations, the coefficient $\mathrm{bq}/\mathrm{s}$ was replaced with a single coefficient set equal to 2.5. This value was tuned, together with the empirical parameter ${\mathrm{a}}_1,$ finally set to 0.0003, to match several experimental observations [35,37,38] according to which the effective diffusion Squire/Lamb constant $\mathsf{\delta }\approx 8$ for small scale helicopter rotors.

Figure 6 illustrates the decay of the normalized vortex circulation with the wake age obtained by applying the decay model of Equation (9) and using the value of 2.5 for the decay coefficient. Despite the slope is slightly lower than the measurements reported in Ramasamy et al. [37], a good agreement with the experimental results can be observed.

Figure 7 shows the predicted growth of the vortex core radius as a function of the wake age obtained by applying Equation (17) and using Equation (18) for the diffusion parameter $\delta$ and Equation (19) for the circulation decay. The picture highlights the close agreement of the calculated parameter $\delta$ with the constant value of 8, which is typical for small scale helicopter rotors. The slight increasing deviation with the wake age from the value of 8 is produced by the application of the decay model in the vortex circulation.

4. Results

PIV measurements and BEM numerical simulations were carried out with the main purpose to investigate the structure and development of the rotor blade tip vortices.

4.1. Numerical Test Set-Up and Computational Resources

Each of the four-rotor blades was discretized by 40 chordwise panels and 23 spanwise panels. No rotor hub nor any other body (such as the motor, the controller, etc.) was modelled. Fourteen rotor revolutions and fourteen wake spirals were required to get a sufficiently long wake and a converged solution. The time discretization adopted corresponded to an azimuth step equal to 2 deg. A computational acceleration was obtained by running the code in parallel mode via the OpenMP API and using 32 cores. The run required about 8 h of elapsed computational time.

4.2. Rotor Wake Ensemble-Averaged Flow Field

The ensemble-averaged velocity field showed the shear layer region surrounding the rotor downwash wake.

The comparison between the experimental results, Figure 8a, and the numerical predictions, Figure 8b highlighted a general similarity but with two main differences:

• in the experimental PIV the origin of the shear layer is identified at about x/R = 0.96 and slightly above z/R = 0 because of the deflection produced by the blade elasticity. Instead, the numerical results show the origin of the shear layer exactly at x/R = 1 and z/R = 0, and this is because the blade was modelled as a fully rigid body;
• the diffusion produced by the viscous effects causes a marked thickening of the experimental shear layer moving downstream from the rotor disk, while the dissipation produces a reduction of the velocity magnitude which is already visible at around z/R = −0.7. These effects are less visible in the numerical results, despite diffusion and dissipation models were applied in the simulations.

A deeper investigation of the aforementioned differences was obtained by comparing the experimental and numerical z component of the induced velocity at several stations below the rotor disk in a region extending up to one radius below the rotor, Figure 9. The numerical results shown in the figure were evaluated at the several azimuthal stations of the last rotor revolution and time-averaged, whereas the velocity fluctuations, due to the flow field unsteadiness, was evaluated and represented in terms of RMS bars. The experimental PIV measurements were made in a fixed vertical plane during about 1450 rotor revolutions. The values were then ensemble-averaged and the velocity fluctuations were represented in terms of RMS bars.

The numerical predictions show a satisfactory agreement with the experiment in the radial region of the blade included from the root cut-out (16%) to the position $\mathrm{r}/\mathrm{R}=80\%$, where the maximum of the inflow is measured. In the radial region, where the tip vortex roll-up produces its greater effect, the numerical results show an upwash that is not present in the experiment. Furthermore, discrepancies between the numerical results and the experiment can also be observed in the region of the rotor hub, not modelled in the numerical simulations.

Finally, the different slope of the derivative $\partial w/\partial x$ between the experiment and the numerical simulation around 80–100% of the blade radius was further analysed by comparing the PIV measurements made with a finer resolution ($\Delta \mathrm{x}=0.47\mathrm{mm}$), with those at a resolution of ($\Delta \mathrm{x}=2.7\mathrm{mm}$) and with the numerical results evaluated on a grid having a resolution of $\Delta \mathrm{x}=\Delta \mathrm{z}=0.6\mathrm{mm}$. The results showed that the slope evaluated by the finer PIV measurements are much closer to the numerical results but this happens up to z/R = 0.4. More downstream, the slope of the finer and coarser PIVs remains almost unchanged.

4.3. Blade Tip Vortices

The application of the ${\mathsf{\Gamma }}_2$ method to the PIV and BEM instantaneous velocity fields highlighted the presence of several tip vortices in the measurement regions. Figure 10 shows an instantaneous velocity field in the ROI immediately downstream of the rotor disc. Three tip vortices can be counted and the vortex centres detected by ${\mathsf{\Gamma }}_2$ method are shown. Once the centres of the tip vortices were detected, the non-dimensional tangential velocities and circulations were calculated and represented versus the non-dimensional vortex radius (Figure 11).

The results have shown a decrement of the swirl velocity that follows a second-order trend as the distance from the rotor increases. The circulation decreases as the distance from the disk increases.

The PIV measurements were not phase-locked with the rotor so that a direct comparison between instantaneous velocity fields was not possible. The BEM/PIV comparison was carried out on the tip vortex characteristics versus the distance from the rotor disk. The growth of the normalized vortex core radius with the distance from the rotor disk is shown in Figure 12. The comparison between the experimental and numerical results show that in the latter case the growth is slightly over-estimated.

The distribution of the experimental vortex centres, Figure 13, showed the highest concentration in the proximity to the rotor disk and that their location was enclosed in the shear layer region. Moving downstream, the data scattering increased distributing the centres of the vortical structures both outside and inside the ensemble-averaged downwash, whereas the concentration of vortex decreased due to the fading and/or merging of vortices. The same representation for the numerical vortex centres showed an extremely narrow region of the shear layer and an almost full match between the boundaries of the shear layer region and the path of the vortex centres.

Additional characterization of the instantaneous tip vortices was obtained by evaluating the composition of the vorticity orthogonal to the PIV plane $\mathsf{\omega }={\left(\nabla \times \mathrm{V}\right)}_{\perp }$ and comparing the experimental results with the numerical simulations. Figure 14 illustrates the results of this comparison. A cut-off of the vorticity module was set at 1500 $1/\mathrm{s}$ to remove as much as possible the flow field small-scale turbulence and to have a more concentrated representation of the tip vortices. PIV results generally show stronger vorticity. The blade flexibility in the experiment generates the tip vortex at a higher and more inboard position with respect to the BEM result ($\mathrm{x}/\mathrm{R}\approx 0.96;\mathrm{z}/\mathrm{R}\approx 0\mathrm{vs}.\mathrm{x}/\mathrm{R}\approx 1;\mathrm{z}/\mathrm{R}\approx -0.02)$. The presence of widespread vorticity along the blade span around x/R = 0.8 and z/R = 0 can be observed in the PIV due to the blade passage.

The numerical simulations allowed the visualization of the entire wake structure generated by the rotor. In particular, it was possible to visualize the tip vortices and to associate them with the generating blades. Figure 15a shows the numerical tip vortex structure up to a distance of z = −1.5 R below the rotor disk. Each of the four colours corresponds to a vortex filament generated by the relative blade. The image illustrates that the vortex structure keeps a geometrical regularity until a distance between z/R = −0.5 and −0.6, after which the wake starts showing more chaotic behaviour. A second aspect that is highlighted in the figure is the pairing phenomenon according to which the mutual positions of the vortex filaments tend to interchange after a first rotor revolution. More specifically, looking at the sequence of the tip vortices at about x/R = + 1.0, it is: cyan-blue-green-red during the first revolution but becomes cyan-green-blue-red during the second revolution. This means that the vortices blue and green begin to roll up with each other until when they have completely interchanged their position after one revolution. Analogously, similar behaviour can be observed at about x/R = −1.0, for the vortex filaments red and cyan. This mechanism contributes to increasing the flow turbulence from the third revolution onward.

Figure 15b shows the association between the numerical vortex filaments and the generated vorticity in the plane corresponding to the PIV measurements. The pairing phenomenon is clearly visible after the second revolution between the green and blue vortex filaments, highlighted in the figure by the black circle, with the vorticity intensity of the first one being smaller than that of the second one. The dashed black circle shows that during the third revolution the green and blue vortex filaments tend to coalesce.

Figure 16 shows a comparison between the experiment and the numerical simulation in terms of the variation of the non-dimensional net circulation of all the identified tip vortices versus the distance from the rotor disk z/R. The net circulation was determined at a distance of 0.25 c from the vortex centre, and by assuming axisymmetric flow in the reference system moving with the vortex core, following the specification in Ramasamy et al. [37]. The experimental data presents three different regions characterised by different slopes, a near-field region comprised between the rotor disk to z/R = −0.23, a mid-region with z/R ranging between −0.23 and −0.75 and a far-field region from z/R = −0.75 to −1.1. In the near-field region, the experimental/numerical comparison shows a reasonable agreement in terms of slope, while an intensity difference of about 17% arises. In the mid-region, where the pairing phenomena start to occur (z/R = −0.5), the experimental slope is almost doubled with respect to the near-field region, and the intensity dissipation is larger than the numerical data. Further downstream, the experimental slope realigns to the numerical trend but with smaller values.

The change in slope of the experimental data in the mid-region and their subsequent scattering at distances z/R < −0.25 might be explained by an increase in dissipation due to the growing of turbulence and by the process of pairing and coalescence of vortex filaments, respectively. This process was also mentioned in the explanation of the BEM tip vortices path of Figure 15a. The numerical results resemble the trends illustrated in Figure 6. The merging of co-rotating vortices, which produce a higher circulation, likely causes the presence of higher values in the region $\mathrm{z}/\mathrm{R}\in \left[-0.8;-1.0\right]$.

Finally, Figure 17 shows a comparison between the experiment and the numerical simulation in terms of the maximum values of the non-dimensional swirl velocity of all the identified tip vortices versus the distance z/R. In this case, the experimental vortices also show three distinct zones with a net difference in terms of slopes and intensities. The reason for this behaviour has the same explanation as for the net circulation trend discussed in Figure 16. The numerical results show a typical hyperbolical decay of the velocity intensity according to the model of Equation (16) combined with Equation (17). An interesting match with the experiment can be observed for distances from the rotor disk lower than z/R = −0.5. The presence of higher values in the region z/R ∈ [−0.8; −1.0] is likely caused by the merging of co-rotating vortices which produce a higher velocity swirl, as also mentioned for the net circulation.

5. Conclusions

Research activity was carried out to characterize experimentally and numerically the blade tip vortices of a small scale four-bladed isolated rotor in hover flight. The investigation of the vortex decay process during the downstream convection of the wake was addressed. 2C-2D PIV measurements were carried out below the rotor disk down to a distance of one rotor radius. The numerical simulations were aimed at assessing the modelling capabilities and the accuracy of a free-wake Boundary Element Methodology.

The Γ₂ vortex centre detection criterion was applied both to the experimental and numerical results. Once detected in the centres, tip vortices were characterised in terms of vorticity, circulation, swirl velocity, core radius and trajectory.

The rotor wake mean velocity field and the instantaneous vortex characteristics were investigated. The experimental/numerical comparisons showed a reasonable agreement in the estimation of the mean velocity inside the rotor wake, whereas the BEM simulations predicted and under-estimated the effect of the diffusion thus generating a smaller shear layer region with respect to the experiment.

The numerical results provide a clear picture of the filament vortex trajectory interested in complex interaction starting at about a distance of z/R = −0.5.

The time evolution of the tip vortices was investigated in terms of net circulation and swirl velocity. The PIV results showed similar behaviour for both quantities. They showed a linear mild decay up to the region interested by vortex pairing and coalescence, where a sudden decrease, characterised by a large data scattering, occurred. The numerical modelling predicted a hyperbolic decay of the swirl velocity down to z/R = −0.4 followed by an almost constant decay. Conversely, the calculated net circulation showed a gradual decrease throughout the whole wake development.

The comparisons showed discrepancies in the region immediately downstream the rotor disk but significant similarities beyond z/R = −0.5.