Abstract
This paper presents the simulation of a two-rotor aircraft in different geometric configurations during hover flight. The analysis was performed using an implementation of the unsteady vortex-lattice method (UVLM). A description of the UVLM is presented as well as the techniques used to enhance the stability of results for rotors in hover flight. The model is validated for an isolated rotor in hover, comparing numerical results to experimental data (high-Reynolds, low-Mach conditions). Results show that an exclusion of the root vortex generates a more stable wake, without affecting results. Results for the two-rotor aircraft show an important influence of the number of blades on the vertical thrust. Furthermore, the geometric configuration has a considerable influence on the pitching moment.
1. Introduction
For many years, aeronautics and aerodynamics research communities have shown interest in vertical and short take-off and landing (V/STOL) aircraft. Their ability to reach unconventional areas make them ideal for certain applications, such as search and rescue missions. However, few configurations have been successful, as it is difficult to understand the aerodynamic interactions during hover, transition (from vertical to forward flight and vice versa), and near-ground flight [1]. Rotorcraft have been most successful in this area. Even if helicopters are predominant VTOL aircraft, their complex aerodynamics are yet to be fully understood.
Rotor wake modeling plays a fundamental role in the design of V/STOL aircraft. Accurate prediction of loads on rotor blades is necessary in assessing rotor performance. Blade-vortex interactions are associated with intrusive external noise and vibration of the blades, thus having a strong impact on the fatigue life of the aircraft and its environmental footprint [2]. Rotor-fuselage interaction affects the overall drag of the aircraft [3], but in most cases it is computationally prohibitive to do complete aircraft simulations, which is the reason why results from rotor wake modeling are often used in the study of fuselages [4, 5].
Low-order Lagrangian numerical methods based on potential flow theory, like the source/doublet panel and vortex-lattice method, are an inexpensive alternative to full Navier-Stokes or Euler equation computational fluid dynamics (CFD) approaches, where the complete domain of the fluid needs to be discretized into a grid and the corresponding partial differential equations are solved. In vortex-lattice methods, lifting surfaces are approximated as thin surfaces represented by vortex sheets of unknown circulations. Since the velocity field of the potential flow must always be tangential to the surface, a nonpenetration boundary condition is imposed at the centroids of the sheets.
Modeling wake structure is the most difficult part when simulating rotary-wing aerodynamics. Initial experiments done by Gray [6, 7] determined that the wake is formed by a high-strength tip-vortex and an inboard vortex sheet where all vorticity is confined. A common problem when using CFD codes is that vorticity at the wake diffuses too quickly due to numerical dissipation [8]. Different solutions that have been proposed, such as adaptive mesh refinement and overlapping grids, tend to increase computational cost without eliminating numerical dissipation effectively. For this reason, vortex methods have been strongly recommended for the analysis of rotary wings because of their ability to model the wake structure and predict loads that agree with experimental data [9].
A special challenge is presented when trying to model the wake of a hovering rotor. A theoretical analysis done by Gupta and Loewy [10] showed that the helical shape of the rotor wake in hover flight has a wide range of unstable modes of motion in both single and interdigitated helices. Bhagwat and Leishman [11] used eigenvalue analysis on the aerodynamic stability of helicopter rotor wakes in hover to show the intrinsic instabilities that appear due to possible unstable deformation modes. These instabilities have made it very difficult for researchers to obtain converged solutions for rotor wakes in hover, which is the reason why different approaches have been taken to circumvent this problem.
Time-marching free-vortex wake models that capture unsteady aerodynamics tend to be unstable in the case of modeling hover flight. In these methods, induced velocities at each fluid particle in the wake are calculated at each time step and the new position of the wake is obtained integrating over time. Initial attempts made by Scully [12] to model rotor hover performance using time-marching schemes failed to obtain convergent results. In order to impose periodicity of the wake, he developed a technique in which an initially prescribed wake is distorted at each time step in an iterative manner and the far wake structure is extrapolated. Miller and Bliss [13] developed a method based on enforcing a periodic boundary condition on the wake over one rotor revolution, which allowed for obtaining converged solutions for rotor wakes in hover and low advance ratio. Pseudoimplicit methods developed by Leishman et al. [9] rewrite the convection equation that describes the motion of the vortex filaments into blade-fixed coordinate systems, resulting in a partial differential equation (PDE) in terms of azimuthal position of the blade and wake age of the vortex filament. The resulting PDE is then solved using finite-difference schemes that introduce an additional term which reduces instabilities and imposes periodicity on the wake. Satisfactory results were obtained with this method for different flight regimes, including hover. To reduce the effect of this instability in time-marching codes without imposing periodicity, high-order predictor-corrector time integration has been implemented to predict performance during hover and maneuvering flight [14]. Quackenbush et al. [15] present an approach based on influence coefficients which proved to diminish wake instabilities. Chung et al. [16] implemented a method based on a slow start of the rotational movement. It reduces the nonphysical instabilities in the initial wake by eliminating the strong circulation of the root vortex generated by an impulsively started rotor.
Vortex-lattice methods based on time-marching schemes have been used successfully in simulation of loads distribution and wake structure of rotary wings. Röttgermann et al. [17] coupled the vortex-lattice method with a field panel method to compute the transonic nonlinear effects at the blade tip. Long and Fritz [18] presented an unsteady vortex-lattice method to study flapping flight as a complex flow behavior. Then, a panel method free-wake code was presented by Roura et al. [19] where an important division of the wake (near and far) was implemented in the code. The near wake is modeled by a vortex lattice that represents the full inbound vortex sheet and tip-vortex while the far wake is modeled by a single concentrated tip-vortex, thus reducing computational cost. Tan and Wang [20] present an unsteady panel time-marching method which is more efficient compared to traditional panel/vortex rings free-wake method.
The present paper discusses the results from the simulation of a “generic” two-rotor VTOL aircraft in hover flight, obtained using the General Unsteady Aerodynamics Vortex-Lattice Method (GUAVLAM) code, which is an implementation of the unsteady vortex-lattice method (UVLM) developed by Konstadinopoulos et al. [21] and Preidikman [22]. Because the method was designed for general aerodynamics, no periodicity is imposed. The purpose of this study was to analyze the influence of the orientation of the rotors relative to the fuselage on the aerodynamic performance of the aircraft in hover flight. Wake stability was improved by implementing a slow-starting method and neglecting the effect of the root vortex, producing a smoother evolution of unsteady loads. The model was validated for the case of an isolated rotor in hover flight using experimental data from Caradonna and Tung [23]. Results show pitching moment is greatly influenced by the orientation of the rotors, while the number of blades per rotor has a greater influence on vertical thrust. Aspect ratio and number of blades were also numerically studied by Steinhoff and Ramachandran [24]; their results show the same tendency measured by Caradonna and Tung.
2. Aerodynamic Model
The unsteady vortex-lattice method (UVLM) was used to model the aerodynamic loads on rotor blades. This method assumes an inviscid, irrotational, and incompressible flow around the body; therefore, it does not account for viscous effects, flow separation at the leading edge, or compressibility at high Mach numbers.
Given the previous assumptions, a velocity potential can be defined; that is, , such that the continuity equation in the inertial frame becomes the Laplace’s equation:
The basic solutions to (1) are the source () and doublet () and the general solution is found by integrating the contribution of the basic solutions distributions over the body’s surface:
In UVLM, the source term is neglected in the case of thin surfaces. Constant-strength doublet panel is equivalent to a closed vortex lattice with the same strength of circulation (). Lifting surfaces are discretized into a lattice of short, straight vortex segments that form closed triangular or quadrilateral vortex loops, as illustrated in Figure 1. Such panels represent the bound vortex sheets [21] where vorticity is confined. The no-penetration boundary condition is enforced at the geometrical centroids of the sheets, called control points.
(a)
(b)
The velocity induced by each vortex segment () is computed using the Biot-Savart law (see equation (3)). A viscous core radius () is introduced to eliminate the singularity due to the velocity field induced by the Biot-Savart law at a point on the vortex segment. A velocity profile proposed by Vatistas et al. [25] was used in the current implementation. According to Leishman et al. [26], the Vatistas vortex model (see (3)) with best fits their experimental data. Graphical representation of the B-S law is depicted in Figure 2.
Spatial conservation of circulation is enforced by assuming a constant circulation (), at each closed vortex loop of index . In order to satisfy the no-penetration condition, the sum of the velocities induced by the vortex segments of the closed loops at the control points of each panel should be tangential to the surface (nonpenetration condition):where the subscript means velocity induced by panel on panel , the subindex is the velocity of the surface at the control point of panel , and is the total number of vortex sheets, including those on the wake. This summation can be separated into velocities induced by the bound vortex sheets (LHS of (5)), velocities induced by free vortex sheets in the wake ( in RHS of (5)), and the free stream velocity . The problem is expressed as a linear system of equations with a dense matrix:where is the normal component of the velocity induced on panel by panel with u unit circulation, is the velocity induced by the wake on panel , and is the number of bound vortex sheets.
The circulation of the vortex segments is calculated as the difference between the two adjacent closed vortex loops:
In the current implementation, all vortex segments are straight and finite. More details on the UVLM can be found in [21, 22].
2.1. Unsteady Kutta Condition
The Kutta condition states that flow leaves tangentially the trailing edge; that is, velocity at the trailing edge is finite. The unsteady Kutta condition is imposed in UVLM by shedding the vorticity generated along a sharp trailing edge into the wake [21]. This condition forces the pressure jump across the trailing edge to be zero, allowing the fluid to leave the surface smoothly. The nodes on the wake are allowed to move freely with local fluid particle speed, rendering the wake force-free. Although some researchers use predictor-corrector [27] and/or high order schemes [14] to calculate the positions of the wake at each time step, the present method uses the explicit Euler method for the convection of the fluid particles of the wake.
2.2. Slow-Starting Method
A slow-starting method was used for the rotation of the blades. Chung et al. [16] showed that, for free-wake calculations on rotor blades, an impulsive start, initial rotation at the final angular velocity, causes nonphysical instability of the initial wake. Also, the strong root-vortex circulation causes vortices near the root to move upward because the velocity influence from the root is larger than the downward velocity. In some applications, such as CFD modeling, a solid body is placed in between the blades to stop upward flow [28]. The slow-starting method was proposed to overcome the nonphysical instabilities and achieve better convergence on simulations. This method consists of starting rotation at an angular velocity near zero, which is gradually incremented until the desired angular velocity is achieved. This was implemented by specifying the number of time steps () taken for the angular speed to reach the desired value . The magnitude of the angular speed is increased linearly as a function of the current time step () according to the following:Physically, the slow start of the rotational movement could be seen as the rotational inertia that every rotor would have to overcome before achieving a final, steady angular velocity.
2.3. Treatment of the Root Vortex
Miller [29] showed that, for hover rotor performance analysis, the effect of the root vortex in the wake was negligible. The current implementation allows for the inclusion or exclusion of the root vortex. When root vorticity is excluded, the circulations of the root vortices () are forced to be zero throughout the simulation. Even though this is strictly a violation of the conservation of circulation, zero vorticity at the root vortex is achieved on a hovering rotor at steady state according to Miller [29]. Therefore, this approximation should only affect the initial (transient) part of the solution and should not affect the steady-state solution.
2.4. Aerodynamic Forces and Moments
Pressure distribution can be obtained from the velocity field using the unsteady Bernoulli’s equation:
For thin lifting surfaces, the pressure jump is calculated at the control points of each panel using nondimensional quantities, so that (8) can be rewritten aswhere is the nondimensional mean velocity at the control point , is the tangential velocity jump at panel , and is the loop circulation at the panel. The pressure distribution is numerically integrated along the surface to yield the acting force () and moment () coefficients: where is the nondimensional position vector of control point from the center of rotation normalized by the characteristic length (). is the nondimensional area of panel normalized by the characteristic area .
The thrust coefficient is given by where is the thrust force, is the wingtip velocity given by , and is the area of the rotor disc. Thrust coefficient () is obtained from the force coefficient () by taking the component in the direction of thrust ( direction) as follows:
3. Results
Numerical results are shown for hover flight of an isolated helicopter rotor and for a generic two-rotor aircraft. Numerical results of the helicopter rotor are compared to experimental data in order to validate the model. The forces and moments at different geometrical configurations of the two-rotor aircraft are shown.
3.1. Helicopter Rotor in Hover Flight
The rotor geometry was taken from [23], where a rotor with two untapered and untwisted rectangular blades of aspect ratio is modelled. The blades were configured at a collective pitch of , , and , with a precone angle of . Each blade was discretized with 15 panels in the spanwise direction using a cosine spacing and 6 panels in the chordwise direction using a uniform spacing.
The current method shows satisfactory results when compared with experimental data, especially for the test case with blades rotating at 1250 RPM, which corresponds to a tip Mach number of 0.439 [14]. Figure 3(a) shows the spanwise distribution of sectional thrust. It can be seen that the numerical prediction of the implemented model is in good agreement with the experimental data except near the tip of the blade, where the thrust is overpredicted. Similar results were obtained by other authors using higher-order methods [14].
(a)
(b)
Results for the thrust coefficient of the rotor at different collective pitch angles are shown in Figure 3(b). Acceptable correspondence between numerical and experimental results is observed, with an average error of less than 6.0%. The values were averaged taking the last 3 rotor revolutions after a total of 9 revolutions.
The effect of the presence of the root vortex was also studied. It was found that the exclusion of the root vortex in the wake had negligible effect on the “steady state” result of thrust coefficient, converging to an average value similar to that achieved with the model with root vorticity. The transient part of results also seems unaffected. Furthermore, simulations without root vortex showed a more stable behavior of the thrust coefficient, as seen in Figure 4, for different collective pitch angles. A comparison of the two wake structures is shown in Figure 5. Simulations without root-vortex circulation show a smoother evolution of the wake compared to results from the model with vorticity at the root, reducing the initial upwash caused by the strong root-vortex circulation, also present in [16]. It is also clear that the no-root vorticity model achieves a stable near wake structure at lower revolutions than the root vorticity model.
(a) 3 revs. RV
(b) 6 revs. RV
(c) 9 revs. RV
(d) 3 revs.
(e) 6 revs.
(f) 9 revs.
3.2. Two-Rotor Aircraft in Hover Flight
The tested aircraft consists of two rotors with rectangular, untwisted blades of AR = 10. The fuselage consists of an elliptic body of major radius and minor radius , where is the rotor radius. Simulations were run with 2, 3, and 4 blades per rotor for 6 to 8 revolutions depending on the convergence of the forces. The collective pitch angle is . The origin of the reference frame is set between the rotors, above the fuselage, from where moments are measured. For the nondimensional analysis, the rotor radius and the wingtip velocity were chosen as characteristic length and velocity, respectively. This means that the nondimensional rotor radius and angular velocity are and . Force and moment coefficients are nondimensionalized by wingtip velocity, . The time step was chosen such that at each step the rotor rotates . The loads on the fuselage are not taken into consideration.
The positions of the rotors are moved using Euler angle rotation in a 3-1-2 sequence, that is, the first angle rotates about the -axis, the second angle rotates about the -axis, and the third angle rotates about the -axis. The initial geometric configuration is at angles . Other three configurations were tested, each rotating about a different axis, as shown in Figure 6.
(a)
(b)
(c)
The evolution of the force and moment coefficients for the configuration with 3 blades is shown in Figure 7. It can be seen that an oscillation of the aerodynamic loads is caused by the interaction of the blades with the fuselage. Also, slight pitching moment occurs due to the orientation of the rotors.
(a)
(b)
Averaged values of the force and moment coefficients were obtained for last two rotor revolutions. These values are shown in Tables 1, 2, and 3. It is clear that, in most cases, force in the direction is near zero (~10−5), which is expected. A slight force is generated in the direction for most cases (~10−4), except for the configuration, where a significant force appears in the direction (~10−2) due to the orientation of the rotors. Predominant force component is in the vertical () direction (~10−1), as expected.
The vertical force coefficients of the different test cases are shown in Figure 8(a). It can be seen that geometrical configuration has minor influence on vertical thrust, while the number of blades has major influence. Numerical results show that the base case generates a small negative moment coefficient, as well as the configuration. However, the remaining configurations generate significant positive moments about the -axis, as shown in Figure 8(b). The wake structure of configuration is illustrated in Figure 9.
(a) Vertical thrust
(b) Pitching moment
4. Conclusions
Numerical results were presented for the hover flight of an isolated rotor and a generic two-rotor aircraft. These results were obtained by the GUAVLAM code, which is an implementation of the unsteady vortex-lattice method for general aerodynamics.
The geometry of the isolated rotor in hover presented in this paper is based on the experimental setup from [23]. Aerodynamic loads from the numerical results are in good agreement with experimental data. Although the method overpredicts sectional thrust near the tip of the blade, the thrust coefficients for all three collective pitch angles tested show good correspondence with experimental data. Numerical data also shows that exclusion of root vorticity generates a smoother wake and has no influence on results at steady conditions. This suggests that exclusion of vorticity in this region is an approximation that stabilizes the solution without compromising final thrust and load distribution. It should be noted that this cannot be generalized for the case of highly twisted blades where root vorticity may be a significant feature of the flow and excluding root vortex could have a negative effect on final results.
A generic two-rotor aircraft was modeled in hover flight with 2, 3, and 4 blades per rotor at a collective pitch angle of . For each geometry, 4 different configurations were tested, changing position and orientation of the rotors. The behavior of forces and moments agree with expected results: force component in the -axis is dominant in all configurations, while other components tend to zero, except for the case of configuration where the component on the -axis is significant. Also, there is a considerable positive pitching moment () for the and configurations, and a relatively small negative for the remaining configurations (including the base case). It was found that the number of blades has a major influence on vertical thrust, while the geometrical configuration has a greater influence on the pitching moment. The current model is able to capture the blade/fuselage interaction which is observed as an oscillatory motion of unsteady loads on the rotor.
Attempts were made to obtain converged solutions without using slow start. However, this was not possible since the strong starting-vortex shed from the trailing edge on the first time step tended to cause numerical instabilities when it was cut across by the on-coming blade. This coincides with the observations made by Chung et al. [16]. The slow-starting method proved to be an essential part of modeling the rotor wake in hover using a low-order time-marching method like the one used in the present study.
In conclusion, the complex case of a rotor wake in hover was modeled successfully using a low-order method thanks to the use of simple solutions such as the slow-starting method and the exclusion of the root vortex. No higher-order time integration or relaxation methods were necessary. The test cases also show the versatility of the code for modeling rotorcraft in different configurations, with and without fuselage, being able to capture unsteady loads.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was supported by the Pablo Neruda Program and Vicerectoria de Investigacion at Universidad de los Andes. The authors gratefully acknowledge the contribution of Centro de Computacion Avanzada-MOX of the School of Engineering at Universidad de los Andes for providing all the computational resources.