VPM/CFD-Based Research on Rotor Performance and Loads of Individual Blade Control Rotor System

This paper aims to explore the effect of individual blade control (IBC) on aerodynamic performance of helicopter rotor and explain its formation mechanism. For this purpose, the vortex particle method (VPM)-computational fluid dynamics (CFD) coupling method was proposed to calculate rotor aerodynamic performance under open-loop IBC active control. Specifically, the near-blade flow field was calculated by the CFD method, while the far-field flow field was solved by the VPM method. In this way, the entire flow field was computed through the information interaction between the two calculated fields. Then, the UH-60A rotor was selected as an example to verify the established VPM/CFD method. First, the proposed method was proved valid; then, the effect of control frequency and phase on the helicopter performance was analysed under different forward flight conditions; finally, the mechanism of IBC control was examined by comparing the lift coefficient distribution and the induced inflows of the optimal control and the worst control. The results showed that proper IBC control parameters can lower the required power of the rotor to some extent, but the optimal control parameters vary with flight states. Comparatively, the lift distribution is more even and the induced flows are less fluctuating under optimal control than under worst control.


INTRODUCTION
For a helicopter, it is more difficult to maintain stability in forward flight than in hovering. The reasons are as follows. During the forward flight, the blades are in a periodically changing aerodynamic environment. In the near-blade region, there exist prominent nonlinear flow phenomena, such as separate flow and shock wave. Meanwhile, the strong tip vortex in the flow field seriously disturbs the movement and induces a highly nonlinear velocity field. When the blades move close to the vortex, the blade-vortex interaction (BVI) is formed, producing nonlinear oscillations of aerodynamic loads. The ensuing vibration and noise make it difficult to optimize the aerodynamic performance of the helicopter. This calls for the reduction of the nonlinear fluctuations of rotor aerodynamic loads.
In the relevant research, the initial focus was placed on the airfoil, the overall design and tip configuration of the blades. Wilby (2001), Leishman (2006), Dadone (1978) and Yamauchi and Johnson (1983) present new, multi-section airfoils capable of elevating the lift-to-drag ratio and reducing the impact of stalling and shock [1][2][3][4]. Leishman (2006) shows that a reasonable nonlinear twist distribution design will result in the blade span-wise aerodynamic distribution is more uniform [2], and thus enhances aerodynamic efficiency. Leoni (2007) and Desopper et al. (1988) point out that the shape of blade tip affects the tip vortex projectile and the shock wave effect at the blade tip [5,6]. With the development in helicopter design, the above methods are no longer suitable for improving helicopter performance.
Recent years have witnessed the popularity of the active control method for enhancing rotor aerodynamic performance. For example, the independent blade control (IBC), initially a tool to reduce rotor vibration load [7], was found capable of lowering the required power and noise of the rotor under proper pitch frequency, amplitude and phase [8,9]. Hence, the IBC was widely recognised as the future research direction of helicopter control [10,11].
Much research has been done on the effect of the IBC on rotor load. Norman et al. (2009) carried out an IBC control experiment on UH-60A rotor system in the National Full-Scale Aerodynamics Complex (NFAC) wind tunnel, proving that the IBC can effectively reduce the required power of the rotor under certain parameter settings [9]. During the numerical analysis of the IBC control system, Yeo et al. (2011) integrated the free-wake method (FWM) and the FWM/ computational fluid dynamics (CFD) coupling method into the numerical model for flow field calculation, and compared the calculated results with the experimental data [12,13]. The comparison shows that the FWM/CFD coupling method is more accurate than the FWM in most cases. Of course, the FWM/CFD falls short of expectations in accuracy and versatility.
Recently, the high-precision vortex method has been coupled with the CFD for rotor calculation. Based on viscous vortex particle method, He et al. (2009) and Adesanya et al. (2017) establish the vortex particle method (VPM)/CFD coupling calculation method, and achieve good results with this coupling method in rotor calculations [14,15].
In view of above, this paper attempts to apply the VPM/CFD coupling method to calculate rotor aerodynamic performance under open-loop IBC active control, aiming to disclose the exact effect of IBC active control on rotor aerodynamic performance of helicopters and its formation mechanism.

CDF-VPM COUPLING METHOD
In this research, the rotor flow field is divided into the near-blade region and the far-field region. The two regions were solved by different calculation methods. In near-blade region, the various nonlinear features of the flow field were computed by the CFD method based on Eulerian description. In the far-field region, the flow field is dominated by the wakes. Hence, the influence of the wakes was solved by the VPM based on Lagrangian description. This forms the CFD-VPM coupling method.

VPM
In an unsteady flow field of high Reynolds number, the rotor wake, especially the tip vortex, has a compact structure. The vorticity is mainly concentrated within a small "tube", while the external flow field can be regarded as incompressible [15,16]. The incompressible and unsteady viscous flow can be discretized by the vortex particle method, that is, the vorticity is discretized into a blob-shaped vortex with a high vorticity and a certain distribution. By this method, the vorticity at any time t and any position can be expressed as: smoothing parameter designed to avoid the singularity of the Biot-Savart law equation; res h is the discrete resolution; c ε is the overlap factor. For convergence, c ε must be equal to or greater than 1 in actual numerical calculation [17].
In the smoothed discrete vorticity field, the incompressible viscous flow in the external flow field of the blade can be controlled by the following equations: The induced speed can be solved by the Biot-Savart law equation: is Green's function [14].
The governing Eqs. (2) depicts the temporal vorticity change resulted from wake distortion and air viscosity as the vortex moves with the local flow field. The first term on the right side of the equation, p ⋅∇ α u , is the stretching term, an indicator of the effect of tension and bending on local vorticity; the second term 2 p n∇ α is the viscosity term, a measure of the vortex diffusion effect caused by air viscosity [14]. During the solution of the control equations, the two terms were calculated separately in a step-by-step manner.
(1) Solution of the stretching term: where ST is short for the stretching term; [ ] ∇u is matrix containing the elements of ( ) , (2) Solution of the viscosity term: The viscosity term is solved by the particle strength exchange (PSE) method proposed by Mas-Gallic [18].
To accurately capture the tip vortex structure, a large number of vortex particles are required for rotor flow field calculations: over 30,000 for forward flight and over 40,000 for hovering. The above equation is extremely time-consuming, for the magnitude of the direct calculation of self-induced velocity and velocity gradient is O(N 2 ). Hence, the fast multipole method (FMM) [19] was introduced to speed up the computation.

CFD
Taking absolute physical quantities as variables, the Reynolds-averaged Navier-Stokes equations were adopted to solve the near-blade flow field. In the inertial system (x, y, z), the governing equation can be expressed as: where the conservation variable being the velocity vector in the inertial system, p being the density and E being the energy; V is the cell volume, S is the face area; F and G are the convective (inviscid) and viscous flux vectors, respectively . The inviscid terms were computed using a secondorder upwind Roe scheme and the viscous terms were computed using second-order central differencing. To meet the need of unsteady computation, a dual time stepping method was applied and the pseudo time step is explicit five-step Runge-Kutta method. All physical quantities are based on dimensionless parameters of hover tip speed, far-field density and blade chord length.

Coupling Algorithm
The information exchange between computing domains is the key to the successful establishment of the coupling method. The computing accuracy and efficiency hinge on the rationality of the information exchange strategy. Our coupling method covers two computing domains: the CFD and the VPM. The former transfers the vorticity information from the blade to the VPM, while the latter provides the CFD with the field boundary information. The information exchange strategy between the two domains is as follows: (1) CFD to VPM: The vorticity information is transferred from the CFD to the VPM through distributed or centralized vortex source method. The distributed method directly transfers the calculated vorticity distribution from the CFD to the VPM. Despite the definite physical meanings, this method may introduce the numerical dissipation from the CFD to the VPM, and suffer from heavy computing loads.
Therefore, the concentrated method is adopted in this research. Specifically, the CFD method was employed to compute the span-wise load distribution of the blade, and then the vorticity distribution of the bound vortex strength is obtained by the Kutta-Joukowski theorem: where l C is the lift coefficient of the blade profile; b Γ is the bound vortex strength; c is the chord length of the blade tip section; q v is the airflow rate at 1/4 chord.
The blade vortex information mainly comes from the change in bound vortex strength. After the near-blade vortex distribution was obtained, the strength of the new vortex element in the VPM can be determined as: where ω γ is the magnitude of the new vortex ring; b v is the relative airflow rate of the blade profile, including the incoming airflow speed, the rotational speed of the rotor and the waving speed of the blade.
(2) VPM to CFD: There are three common ways to correct the influence of far-field flow field from the VPM to the CFD: first, correction of the blade's local angle of attack by the vortex-induced velocity; second, imposing induced velocity on all mesh grid using an equivalent method; third, direct application of the boundary condition on the outer boundary of the CFD.
In spite of its high efficiency, the first method is low in accuracy due to the errors incurred during the correction of the attack angle. The second method applies the wake effect to all grid points, aiming to directly reflect the strong impact of tip vortex on the flow field. However, the application drags down the efficiency of the method. The third method boasts better physical meaning than the first and second methods, and efficiently computes satisfactory results without scarifying the mesh quality. In view of these, the third method is adopted for this research.
In addition to velocity, density and energy are also part of boundary condition information. Here, the compressibility from the outer boundary of the blade grid is neglected. The density at the boundary was taken as the far-field density. Hence, the energy boundary condition can be expressed as: where M is the local Mach number.

Trim control strategy and flow chart
The input to the IBC active control system can be expressed as: , , . Assuming that the difference between the calculated target state quantity and the target value is ∆Y , the fine adjustment amount of the input quantity can be obtained by the tangent method: (11) where the Jacobian matrix is:  (12) The Jacobian matrix was calculated by incremental method. In other words, the 0.75 ϕ , 1c θ and 1s θ were respectively increased, and the corresponding converged Y was processed by Eq. (12). By the conventional method, the entire balancing process requires multiple iterations and consumes too much time. For better efficiency, the trim strategy of the difference method [20] was adopted for the calculation. As shown in Fig. 1, the core idea of this method is to divide the aerodynamic calculation into two stages, and apply a simpler and faster aerodynamic model in the trim stage. To ensure the accuracy, the results of the simple aerodynamic model were corrected against the aerodynamic force acquired the accurate aerodynamic model. In this way, the computing speed in the trim stage was improved without sacrificing computing accuracy.

RESULTS AND DISCUSSION
The proposed method was applied to the UH-60A rotor. The blade design of the rotor carries the features of advanced helicopters, e.g. sweptback, nonlinear negative torsion and multi-section airfoil. The detailed parameters are listed in Tab. 1.   The following cases were selected to verify the proposed method: low-speed and medium-speed forward flight (Case 8513), high-speed medium-overload (Case 8534) [21], and high-speed forward flight (Case Norman) [9]. The parameters in each case are given in Tab. 2.

Verification of Our Method
The accurate capture of shock wave is the key and difficult point to this research. As a control group, pressure measuring points were arranged across the crosssection of one of the rotor blades in the UH-60A, with the aim to capture the load data. To verify the shock wave capture capability, Case 8534 was selected to analyse the blade surface pressure distribution at the forward side and span-wise location r/R=0.865 (Fig. 3). It can be seen that our method captures the exact pressure change induced by the shock wave, and the resultant pressure distribution curve agrees well with the experimental value. Then, the overall capability of rotor load prediction was validated, to see if our method can accurately calculate the rotor aerodynamic load distribution. Fig. 4 compares the sectional normal force distribution of Case 8513 and Case 8534 at r/R = 0.965. It is clear that the calculated results are consistent with the experimental results in terms of magnitude and azimuth change, although the blade was assumed to be rigid and only the first-order influence was considered in the wiggling dynamics model.

Rotor Performance and Loads for IBC System
Three forward flight conditions were selected to examine the effect of the IBC active control system on the aerodynamic performance of the helicopter. Note that the amplitude of the IBC high-frequency control was maintained at 1.2°. For simplicity, Case Norman was taken as the example of analysis and calculation. Figs. 5 and 6 respectively display the phase sweep results at the frequency of 2/rev and 4/rev. As shown in the figures, the trimmed rotor tension difference is so small as to be negligible, while the required power differs significantly. The required power at 2/rev and 240°is about 2% lower than the non-controlled state. The result is basically consistent with the experimental conclusion in Reference [9]. When the frequency is 4/rev, the optimal control point for required power is about 30, but the control effect is not obvious.  Fig. 7 illustrates the distribution of the optimal and the worst control points for the two frequencies under each pre-flight condition. It can be seen that, for highspeed forward flight cases like Case Norman and Case 8534, the phase distribution between optimal and worst control points is consistent, despite the difference in specific flight parameters; for the medium-speed forward flight case, the phase difference is quite different from that of high-speed forward flight conditions. Meanwhile, the medium-speed flight condition has the most significant power reduction at the optimal control point under the two frequencies, about 6% lower than the uncontrolled condition. In addition, the optimal control points of different frequencies have different phases under the same working condition. Next, the author analysed the IBC's enhancement mechanism of BVI aerodynamic performance. Figs. 8 and 9 show the variation in normal force of blade rotating under the significant control effect of moderate speed Case 8513. Both 2/dev and 4/dev examples show that the blade normal force changes more smoothly in the optimal control state. Fig. 10 compares the time-varying wake-induced inflows under the optimal control and the worst control for Case 8513. The observation point is 0° azimuth at the plane of the paddle and 0.8R radius. It is observed that the inflows are less volatile under the optimal control.

CONCLUSION
In this paper, the VPM/CFD coupling method was proposed to examine the effect of IBC active control system on aerodynamic performance of helicopter rotor and its formation mechanism. Through the analysis, the following conclusions were drawn: First, different rotor advance ratios correspond to different optimal control point phases at the same IBC control frequency; the control effect is more obvious in the medium speed state than the high-speed state. When the rotor advance ratio is 0.15, the power required by the optimal control point is 6% lower than that under the uncontrolled state. Second, the optimal control point phase varies with the control frequencies at the same flight status. Third, after reaching the optimal control point, the induced inflow at the propeller plate becomes gentler, the peak value turns smaller, and the normal lift coefficient of the blade changes more smoothly during blade rotation.