Studying the wake contraction of the flow-field of a rotor in hover

This article presents a comparative study of three numerical algorithms for computation of the flow-field of a rotor in hover and the net static thrust that is produced. The flow-field is induced by a series of vortex rings, modelling the near wake of the hovering rotor and a single semi-infinite vortex cylinder, accounting for the velocity deficit in the far wake. All three models are based on the vortex theory and differ in the choice of the numerical scheme for the estimation of the exact position of the vortex rings, emitted at the tips of the blades of the rotor. Thus, the numerical models perform a real-time simulation of the propagation of the vortex rings in the downwash. The first model uses an Euler-predictor scheme, while the second and third models use respectively first and second order predictor-corrector schemes. The aim of the study is to assess the rapidity and accuracy of each algorithm. For that purpose, the numerical results are compared with the experimental data, obtained from a wind tunnel test of the model rotor. The best results in terms of computational speed and accuracy are obtained with the use of Adams-Bashforth predictor-corrector scheme of second order.


Introduction
Rotors of drones rarely perform long lasting hovers in zero wind conditions with constant performance parameters. This is due to atmospheric disturbances, which require continues thrust adjustments, in order to maintain a stable hover, while compensating for gusts and constant winds. Thus, of great interest is rotor performance during transient or unsteady operational regimes.
The study aims to develop a numerical model, capable of rapid and accurate computation of the thrust of a hovering rotor. Such a tool will make it possible to study rotor performance in unsteady flow conditions, resulting from the transition from one operational regime to another or from phenomenons, such as the tip vortex aperiodicity of a hovering rotor [1].
Previous research [2] performed by the authors, showed that the fastest numerical computation of the performance of a hovering rotor is achieved with the use of the vortex theory and more specifically with the model presented in section 3 "Numerical model". In this previous study [2], the numerical results showed good agreement when compared with the experimental data of a model rotor in hover. The main difference between the previous studies and the current study is that the latter features a free wake approach, while in the former the numerical models had a prescribed arrangement of the vortex elements, which modelled the flow-field of the rotor in hover. The vortex elements are allowed to move freely under the forces applied on them by the resulting induced velocity field, until the vortex system reaches a state of equilibrium. Thus, a greater accuracy is achieved, although the required computational time is significantly increased.
The Miller model is coupled with a BEM code, which allows for rapid determination of the forces acting on the blades of the rotor. This allows for rapid evaluation of the overall rotor performance by providing estimations for the produced net thrust and the required torque for stable hover. Those two can be later compared with experimental data measurements for the model rotor in hover.

State-of-the-art vortex models for a rotor in hover
According to [3], Joukovsky performed a series of theoretical studies on propellers, which set the mathematical foundations for the vortex theory such as applied to rotors. Those works were followed by numerous flow visualization experimental studies on small scale and full scale helicopter rotors in both forward flight and hover. They added more clarity on the structure of the flow-field around a helicopter rotor. These studies also provided the first semi-empirical numerical models for hovering rotors, derived from the obtained experimental performance data and the geometry of the trailing wake.
Jenny et al. [4] and Landgrebe [5] proposed vortex models for predicting the hovering performance of helicopter rotors. These models consist of a helical wake, composed of helical vortex sheets, trailing from each one of the rotor blades. Simplifying assumptions are introduced, which allow for the representation of the helical wake with a finite number of discrete vortex filaments, namely straight-line segmentation of the tip trailing helical vortex. The induced velocity field is computed with the use of the Biot-Savart law. These numerical models are capable of predicting adequately the performance of a hovering rotor, although the simplifying assumption for a rigid helical wake introduces a significant error when the rotor operates in too light or too heavy load conditions, corresponding to values of the coefficients of thrust CT smaller than 0,0025 and bigger than 0,01. Another inconvenience is in order to have a higher degree of accuracy a greater number of straightline segments per complete helical turn must be considered and a greater number of complete helical turns must be taken into account, which significantly increases the volume of the required numerical computations.
Landgrebe [6], Miller et al. [7,8] and Reddy [9] proposed an alternative approach for modelling the hovering rotor by decomposing the trailing wake into several parts, namely: near, intermediate and far. Each of those parts is modelled with a different type and amount of vortex elements with the aim to achieve greater accuracy, while offering significant numerical simplification. The aforementioned papers demonstrate the possibility to obtain a significant degree of accuracy for the evaluation of the performance of hovering rotors with vortex systems, which comprise less than 30 vortex elements in total.
Karpatne et al. [1] demonstrated the capability of a vortex model, in which the wake of the hovering rotor is modelled entirely with a system of vortex rings. In their model the vortex rings are emitted with a specific timestep, depending on the angular velocity of the rotor and are allowed to move freely under the induced velocities resulting from the interaction with the remainder of the previously emitted vortex rings. Reference [1] showed that vortex models are capable of adequately simulating unsteady flow phenomena, thus allowing to study them.

Numerical model
Modelling the wake of a hovering rotor with a combination of different types of vortex elements for different segments of the wake is demonstrated in [1,8,9]. The present study is performed with the vortex model proposed by Miller et al. in [8]. In the rest of the article it will be referred to as the Miller model.
The wake of the Miller model consists of two parts: a near wake composed of 20 circular vortex rings arranged in series downstream of the rotor; and a far wake modelled by a single semi-infinite vortex cylinder, placed behind the last vortex ring of the near wake. Thus the model consists of only 21 vortex elements to be accounted for in the computation of the entire flow-field around the rotor.
The semi-infinite vortex cylinder allows the model to compensate for the velocity deficit formed downstream of a wake when modelled with a series of vortex rings. It allows to stabilize the numerical computation and ensures that the induced velocities in the far wake double those induced in the rotor plane, such as per the theory for the hovering rotor.
The crucial thing for the degree of accuracy and the adequate operation of the model is the proper setup of parameters d0 and d 1 , shown on Fig. 1. The distance d 0 , at which the first vortex ring is placed behind the rotor, plays a key role in the distribution of the induced velocities along the blade. It is referred in [8] as the first blade-vortex encounter position, which represents the distance between the vortex ring emitted by the preceding blade and the following blade. This parameter can be varied, in order to study its influence on the induced velocity distribution on the blades and in the near wake downstream of the rotor.
Distance d1 is the spacing between the last vortex ring of the near wake and the semi-infinite vortex cylinder, which models the far wake. The optimal values of parameters d 0 and d 1 are estimated for a comparison between numerical and experimental data for multiple operational regimes of the rotor.
The Miller model offers the potential for rapid numerical computation when compared to CFD and BEM models due to the necessity to compute for significantly smaller amount of vortex elements. An additional reduction of the volume of numerical computations is possible for hover and vertical flight regimes due to the presence of axis-symmetric flow condition for the rotor. For the studied rotor with two blades, only half of the flow-field can be computed in order to visualize the entire flow around the rotor.

Numerical study
The performed numerical study aims to evaluate the accuracy and computational rapidity of the free wake Miller model, coupled with three different algorithms for the precise estimation of the exact arrangement of the vortex rings in the wake.
Computations are made of the mutually induced and self-induced velocities between the vortex elements modelling the wake. The induced velocities on each vortex element are then used to compute its new position in the wake for the next computational time step Δt. Thus, the vortex rings propagate in the numerical domain and take the necessary geometry arrangement, in order to adapt to the latest flow conditions.
Once the solution has converged and the equilibrium arrangement of the vortex elements is obtained, then the distribution of the induced velocity Vz along the blades is computed. Thus, by knowing the angular velocity of the rotor Ω and the pitch angle of the blades θ, it is possible to determine the aerodynamic forces acting on the blades. Those forces are used for the computation of the net trust T produced by the rotor and the required mechanical power P for the rotor to perform hover. This computational cycle repeats itself for each time step Δt of the simulated time-period. The described algorithm is shown in Fig. 2. By keeping the values of both the angular velocity of rotor Ω and the pitch angle of blades θ constant, the numerical model evaluates the rotor parameters in hover flight. However, the presented program module is also capable of studying the rotor performance in smooth transition from one operational regime to another, as well as for an abrupt step-change of the pitch angle of the blades θ.
For this present study assumed is a stable hover, which is a hover performed in zero wind conditions. The results from this study are presented in section 5 "Results".
The model rotor, which is used as reference geometry in those studies, has two rectangular untwisted blades of NACA0012 profile. The tip radius R of the blades is 288 mm, the root radius at the hub mounting r0 is 65 mm and the chord length c is 25 mm. The pitch angle of the blades θ is adjustable and thus can be varied between wind tunnel experiments. Fig. 3 shows the velocity triangle in a typical crosssection of the blade. The velocity of the axial flow through the plane of rotation, resulting from the work of the rotor, is designated by Vz. In the literature it is referred to as the axial induced velocity and is a function of the operational parameters of the rotor, namely angular velocity Ω and the pitch angle of blade θ. The speed of rotation of the rotor section in the plane is is refferred to as U = Ωy, where y is the blade station radius, which varies from r 0 to R. For the tip of the blade y = R and U = ΩR = V tip . The relative air speed for the cross-section of the blade W is computed according to (1): . (1)

Numerical approach
According to [1], by assuming uniform bound circulation and a single tip-trailing vortex per blade, the vortex strength Γ can be calculated as follows: , where T is the total thrust produced by the model rotor at a fixed angular velocity Ω and pitch angle of the blades θ. R is the tip radius of the blades, V tip = ΩR is the velocity at the tip of the blades, N b is the number of blades of the rotor and ρ is the density of the air. For every change in either the angular velocity Ω or the pitch angle of the blades θ, the total thrust of the rotor T and the bound circulation Г are also changed. Equation (2) is used for the computation of the bound circulation Г.
In [10] the coefficient of thrust C T is defined as: , where T is the total thrust produced by the rotor; ρ is the density of the air; Ω is the angular velocity of the rotor; R is the radius of the rotor and A = πR 2 is the area of the rotor disc. Equation (4) is used for the computation of the coefficient of thrust C T . From Fig. 3 it can be seen that with an increase of the axial induced velocity V z there is a decrease of the angle of attack α. The smaller the angle of attack of the crosssectional airfoil is, the smaller the local thrust increment ΔF z is. Thus, the total thrust of the rotor T decreases and with it decreases the coefficient of thrust C T . Fig. 4 shows the key dimensions, used in the computation of the induced velocity at a point of the flowfield. The y m , y n , z m and z n dimensions are universal for both the cylindrical and contracting wake models. Those are the coordinates of two points, namely M and N. It is considered that the induction happens in point M by a vortex element, situated at point N. Lewis [11] introduces the dimensionless parameters z and y with (4) and (5): Тhose are used for deriving the equations for the axial V z and radial V y induced velocities. in the case of the semiinfinite vortex cylinder, the axial velocity is computed with (8) and in the case of y = 1, (9) is used: In the case of a semi-infinite vortex cylinder, the radial induced velocity is computed with (11): In the case of the system of vortex rings, the axial and radial induced velocities are respectively computed with (14) and (15): � ������ � ����� � � � �� � ��� � ������ � �� ; (14) Equations (8), (9), (11), (14) and (15) contain complete elliptic integrals of the first E(k 2 ), second K(k 2 ) and third kind П(n,k 2 ). The elliptic parameters n and k are introduced with (16) and (17): Equations (18) and (19) are used for the application of first-order Euler-predictor: Equations (20) and (21) are used for the application of second-order Adams-Bashforth predictor: Equations (22) and (23) are used in the corrector step: The first model only updates coordinates y and z for each vortex element by applying only a first-order Eulerpredictor without performing any correction step.
The second model updates coordinates y and z with a first-order Euler-predictor step, followed by a corrector step. This solution scheme is known as the Trapezoidal predictor-corrector scheme of first order.
The third model updates coordinates y and z with an Adams-Bashforth predictor step, followed by a corrector step. It is referred to as the Adams-Bashforth predictorcorrector scheme of second order.
Numerical instability issues are observed if the vortex rings in the near wake are not being initially slightly contacted. Thus, a simple arrangement with linear contraction higher than 10% is found to be enough to ensure stable computation.  The axial induced velocity V z /V tip in the far wake at z/R > 10 doubles its value, when compared with the induced velocity in the plane of the rotor at z/R = 0, only when the semi-infinite vortex cylinder is positioned at a specific distance behind the last vortex ring from the near wake. This distance was found to be equal to half of the averaged distance between the vortex rings in the near wake. Thus d1 = 0,5p, where p is the average spacing between the vortex rings in the near wake. Table 1 compares the coefficient of thrust C T , computed with all three models, against the wind tunnel experiment for 2000 RPM and θ = 6 deg. It presents also the percentage error |Δ| and the required computational time t comp for the simulations. The spacing parameters for both computations are set as follows: d 0 = 0,25p and d 1 = 0,5p.

Conclusion and Future work
The accuracy and rapidity of the numerical approach for the computation of the performance of a model rotor in hover was evaluated with three free wake vortex based models. On the basis of the results from the performed studies, it can be concluded that the free wake Miller model with the second-order predictor-corrector is only slightly slower than the model using the first-order predictor-corrector in terms of computational speed, while providing significantly higher degree of accuracy. The existing percentage of error between the computational results and the experimental data is acceptable. Therefore, the proposed model can be used for preliminary studies and simulations of the expected thrust and overall performance of newly designed rotors.
The results from the conducted studies, presented in this paper, encouraged the authors to use a free wake Miller model with a predictor-corrector scheme of second order, in order to study the thrust of small rotors in unsteady operational condition. The intent is to conduct a series of studies on rotor performance during transient operational regimes, produced by both big and small rates of change in the angular velocity of the rotor Ω; and for both positive and negative rates of change.
It is of special interest to study the changes in the thrust, resulting from a rapid change of the pitch angle of the blades θ for a constant rotor angular velocity of Ω. Such a study will provide greater insights into rotor performance during the transition period, in which the flow-field adjusts to the new operational conditions of the rotor.