, ,
National Key Laboratory of Science and Technology on Rotorcraft Aeromechanics, College of Aerospace Engineering,Nanjing University of Aeronautics and Astronautics, Nanjing 210016, P.R. China
Abstract: A coupled Navier-Stokes/free-wake method is developed to predict the rotor aerodynamics and wake. The widely-used Farassat 1A formulation is adopted to predict the rotor noise. In the coupled method, the Reynolds-averaged Navier-Stokes (RANS) solver is established to simulate complex aerodynamic phenomena around blade and the tip-wake is captured by a free-wake model without numerical dissipation in the off-body wake zone. To overcome the time-consuming of the coupling strategy in previous studies, a more efficient coupling strategy is presented, by which only the induced velocity on the outer boundary grid need to be calculated. In order to obtain blade control settings, a delta trimming procedure is developed, which is more efficient than traditional trim method in the calculation of Jacobian matrix. Several flight conditions are simulated to demonstrate the validity of the coupled method. Then the rotor noise of operational load survey (OLS) is studied by the developed method as an application and the computational results are shown to be in good agreements with the available experimental data.
Key words: helicopter; rotor; coupled N-S/free-wake method; delta trim method; BVI noise
Helicopter rotor noise, especially the blade-vortex interaction (BVI) noise, has a significant impact on ground observers. The BVI noise occurs in low speed and descending flight conditions, where the generated tip wakes remain within the vicinity of the rotating blades. The close interaction between tip wakes and the blades leads to impulsive fluctuation of blade air-loadings and strong BVI noise[1]. Therefore, an accurate numerical model that is able to capture the rotor wakes and predict the fluctuation of blade air-loadings is necessary for simulation of BVI noise. Meanwhile, the trimming procedure should be incorporated in numerical schemes to match the rotor thrust and eliminate the rotor moments. However, in previous studies, the control inputs are usually obtained from experimental data and the rotor has barely been trimmed by computational fluid dynamics (CFD) method in numerical simulation in most of the previous studies. This is because the complexity and computational time comsumption will have a significant crease dramatically with the traditional trimming method. Therefore, it is also necessary to develop a more efficient rotor trim method for obtaining the accurate blade control settings.
There exist three kinds of approaches in the simulation of rotor aerodynamics and tip wakes,e.g.CFD method[2], potential/wake method[3]and the coupled CFD/wake method[4]. The CFD method is capable of capturing the wakes as part of the solution and has been achieved rapid improvements in past decades. But the application of CFD method is limited due to low computation efficiency and severe numerical dissipation[4]. On the contrary, in free-wake method, the rotor wake can be preserved long distance without numerical dissipation. But the free-wake method cannot simulate the viscous effects around blades. The coupled CFD/wake method combines the advantages of CFD method in simulating of compressibility/viscous effects and wake method in rotor wake capturing. In coupled method, the CFD (Euler or Navier-Stokes equation) method is used to calculte the flow-field around the blade and the wake method is employed to capture the tip wake in wake region. It is computationally efficient with a good numerical precision and receives more and more attention. Zhao et al.[5], Berkman et al.[6]developed an N-S/full potential/free-wake coupled method for the simulation of rotor flows, respectively, and the predicted results are in good agreement with the experimental data. However, the complexity of the method is increased by the flowfield transform among the three different computation domains. Afterwards, Sitaraman et al.[7]presented a coupled N-S/free-wake method and also established the field-velocity method for passing the velocity from free-wake to CFD domain. In the field-velocity method, the induced velocities on all the grid points are incorporated via grid movement in CFD solution. Accordingly, the geometric conservation law (GCL) is enforced to satisfy the conservative relations of areas and volumes of the control cells. Accordingly, there is the fault of large calculation quantity for the “field-velocity method”. Sugiura et al.[8]studied the coupled Euler/prescribed wake method, but the prescribed wake has low accuracy of calculation.
As pointed out in Ref.[9], inaccurate numerical results may be obtained without incorporating a trimming procedure in aerodynamic model. In previous rotor aerodynamics studies, the traditional trimming procedure[9]is directly incorporated to the numerical schemes. For example, Kim et al.[10]incorporated the trimming procedure to a rotor CFD solver. To match the aerodynamic parameters to the desired levels, 20 revolutions of CFD solver are run with their solver. The researchers developed an efficient trim method[11], in which the time-consuming calculation of Jacobian matrix is conducted by the simple rotor inflow model. This provides an idea on the establishment of efficient trim method in coupled method.
In view of this, a coupled N-S/free-wake method with an efficient trimming procedure is presented for the calculation of BVI aerodynamics in this paper. In this method, the coupling strategy called “outer boundary correction approach” is established. Furthermore, a delta trimming procedure is developed based on the “delta idea”, In the procedure, the time-consuming computation of Jacobian matrix is performed by the efficient free-wake method instead of CFD method. And the numerical cases at several flight conditions are calculated to validate the method and illustrate the improvement of the proposed method. Then, combing the coupled CFD/free-wake method with the well-known FW-H equation, a methodology for the prediction of rotor BVI noise is established. The BVI noise of the widely-used operational load survey (OLS) helicopter rotor is calculated by the developed methodology. It is shown that BVI noise energy is mainly radiated to the 120° azimuth and 45° elevation.
Fig.1 shows a schematic of the coupled model. In this model, the N-S based CFD method is employed to predict the blade air-loadings and complex aerodynamic phenomena, while a pseudo-implicit predictor-corrector free-wake method[12]is used to capture the rotor tip-wake without dissipation. For coupling algorithm between the two computation domains, an “integrated vorticity approach”[13]and an “outer boundary correction approach” are adopted. The details of the coupled method are described in Fig.1.
Fig.1 Schematic of the proposed CFD/free-wake coupled method
1.1.1 N-S based CFD solver
A three dimensional unsteady compressible N-S equation[14]is used as the governing equation of CFD
(1)
where
(2)
In the above equations,Wis the vector of conserved variables.F(W) andG(W) are the convective (inviscid) and the viscous flux vectors, respectively,qandqbthe absolute velocity and the blade moving velocity, respectively. The blade moving velocity includes blade rotation, flapping and pitching motions;p,ρa(bǔ)nderepresent pressure, density and energy, respectively.Vis the grid cell volume,Sthe grid surface aream,nthe normal vector of grid surface,Φthe viscous stresses, andτthe components of the viscous stress tensor.
The inviscid terms are computed with a second-order upwind Roe scheme[15]and the viscous terms are computed by second-order central differencing. For the time integration, a dual-time stepping algorithm is employed with the LU-SGS scheme[16]to simulate the unsteady flow phenomenon at every pseudo-time step. The Spalart-Allmaras turbulence model[17]is applied to the closure RANS equation.
1.1.2 Free-wake model
The rotor wake is composed of near wake and far wake. In the CFD/free-wake coupled method, the near-wake is simulated directly by CFD method. The far-wake comprises only a single tip vortex filament. It evolves from a roll-up of the near-wake detached from the trailing edge of blade after revolving 30°. The motion equation of a vortex filament can be written as[3]
(3)
wherer(ψ,ζ) denotes the positions of control points,ψthe azimuth at which the blade located,ζthe wake age angle, andV(r(ψ,ζ)) the velocity of control points including free-stream and induced velocities from bound vortex as well as the other vortex filaments. The induced velocities are computed by Biot-Savart law with Scully vortex core model[18]with an initial vortex core radius of 0.2C(Cis the blade chord).
1.1.3 Coupled strategy between two domains
In the coupled method, it is necessary to pass the vorticity from CFD domain to wake domain and feedback the induced velocity from wake domain back to the CFD domain. In the proposed coupled method, the “integrated vorticity approach” is adopted to pass the vorticity from the CFD domain to the free-wake solution. In this method, the radial distribution of sectional lift coefficient is calculated first by integrating the pressure distributions on blade surfaces by CFD. Then the circulation strength of each blade section is calculated using the Kutta-Joukowski theorem[19]
(4)
whereClis the sectional lift coefficient andvthe local flow velocity.
The “outer boundary correction approach” is used to reflect the influence of rotor wake on CFD domain. Different from the traditional “field-velocity approach”[20], in “outer boundary correction approach”, the variables of flow field for only the outer grid is calculated by the following equations
(5)
(6)
(7)
(8)
(9)
where (u∞,v∞,w∞),(ui,vi,wi) are the free-stream velocities and wake-induced velocities, respectively.γis the ratio of specific heat,athe local mach number and ∞ the free stream.
In the coupled method, the induced velocities need to be calculated at every time step. However, the near wake is computed as part of the CFD solution. To avoid the double counting of the near wake, the first 30° of the wake geometry of the free wake are excluded from the induced-velocity calculations.
To simulate the rotor aerodynamics accurately, it is necessary to input the precise blade control settings[12]. An efficient trimming procedure is developed in this paper. The trim means that the desired aerodynamics level (including rotor thrust coefficient, rolling moment and pitching moment coefficient) is matched by adjusting the control settings (including collective and cyclic pitch angle) continually.
The rotor thrust and aerodynamic moment coefficients can be expressed as the function of control settings.
(10)
whereCT,CMx,CMydenote the rotor thrust, the pitching moment and the rolling moment coefficients, respectively. In the motion of rotor, the blade pitch angle can be described using the first harmonic term of the Fourier series, which can be written as
θ(ψ)=θ0+θ1ssin(ψ)+θ1ccos(ψ)
(11)
To trim the rotor, an iterative technique is necessary due to the nonlinear relationship between aerodynamic parameters and control settings. The Newton-Raphson iterative method is used to obtain the correction angles in this paper. The steps obtaining the correction angles include:
(12)
Thirdly, with the Newton-Raphson iterative method, the correction angle of control settings ΔXis calculated. The calculation equations are as follows
(13)
(14)
In this research, the Jacobian matrix is computed by the efficient free-wake method using a lifting-line technique. But, the aerodynamic parameters cannot be predicted accurately by the technique. Accordingly, they are corrected with the results from the coupled method. The trimming method is called “the delta trim method”. Taking the rotor thrust coefficient as a example, the correction equations is
(15)
Using the blade surface pressure calculated by the coupled method, acoustic prediction is conducted by the widely used Ffowcs-Hawkings (FW-H) equation, which is based on the Lighthill’s acoustic analogy. The time domain solution of the equation, which is also called F1A formulation[21], is employed in BVI noise prediction
(16)
(17)
(18)
The flowchart of the calculation of rotor aerodynamics and noise by the preposed CFD/free-wake/FW-H method is given in Fig.2. The solution is conducted as follows: (1) Estimate the control settings and initialize CFD and free-wake solutions. (2) Calculate the bound circulation and to generate the vorticity of tip-wake. (3) Calculate the motion equation of tip wake and to calculate the induced velocity of out grids. (4) Trim the rotor blades until convergence and calculate the blade surface pressures. (5) Calculate rotor noise.
Fig.2 Flowchart of the proposed CFD/free-wake/FW-H method
It should be pointed out that in CFD domain the time step is set to be 0.25°, while the time step of 1° is used in free-wake solution. The rotor wake geometry between time steps in CFD domain is obtained by the fast Fourier transform interpolation method.
For coupled method, the induced velocity is important for the accurate prediction of rotor air-loadings. The induced velocity for a four-bladed rotor is calculated and compared with the available test data[22], which is measured one chord above the rotor plane. The rotor has a forward shaft tilt of -3°. The blade has a linear twist of -8° from the root to tip and a radius of 0.86 m. The wake-induced velocities are simulated for the case, where thrust coefficient is 0.008 and the advance ratio is 0.15. The comparison of induced velocity between the experimental data and the propsed results is given in Fig.3. It can be seen from Fig.3, the proposed results demonstrate agreements with the experimental data. The discrepancies of inflow in central regions may be caused by the presence of fuselage and rotor hub in experiment. The “vortex-induced” upwash effects in the front of the rotor disk is also precisely predicted with the proposed method. Besides, it can be seen that the induced velocities are bigger in the rear of rotor disk from the induced velocities contour. This is because in this flight condition, the wakes are blowed down and backwards to rotor disk.
A numerical case is simulated to discuss the effects of coupling algorithm in coupled method in Fig.4. The C-T rotor[23]in hover flight is used, where the blade tip mach number is 0.794 and the collective angle is 8°. The rotor has two blades with aspect ratio of 6 and a chord of 0.190 5 m. It can be seen that the results by “outer boundary correction approach” agree better with the experimental data. This is because in “field-velocity ap-proach”, the influence of wake-induced velocity on all grid faces in CFD domain is calculated. The result induced velocity on flow field is relatively large[4]. This leads to the blade angle of attack and thus the pressure difference to be small. But, with the increase of rotational speed from the blade inboard regions to the outboard regions, the influence of induced velocity on blade angle of attack diminishes gradually. As a result, the pressure distributions on outboard regions has a better agreement with the experimental data.
Fig.3 The calculated induced velocity distribution
Fig.4 Comparison of different method for the feed-back of wake-induced velocity
By “outer boundary correction approach”, a better accordance with experimental data is obtained at all blade span locations.
The blade surface pressure distribution are calculated for the C-T model rotor in hover case by the full CFD (embedded grid CFD) method[24]and the proposed coupled method. Considering the symmetrical characteristics of flow field in hover case, only one blade is modeled. The grid used in full CFD is shown in Fig.5.
The computation condition is set as follows:thrust coefficient of 0.004 6, hover tip mach number of 0.439. The history of collective angle and thrust coefficient during trimming procedure is shown in Fig.6. The trimmed collective angle and thrust coefficient are almost same with the experiment data[23]after 3 trim cycles.
Fig.5 The generated embedded grid for hover case
Fig.6 The history of collective angle and thrust coefficient with trim cycles
The chordwise pressure distribution for several radial stations are simulated by both the coupled method and the full CFD method. The results show good agreements with measurements for both inboard and outboard span locations as seen in Fig.7.
Under the same condition, the effects of rotor wakes are also analyzed in Fig.7. The blade pressure are calculated by the proposed coupled method, but in the simulation, the effects of rotor wakes are not coupled. It can be seen that without the rotor wake, the blade surface pressure difference (the dashed dot line in Fig.7) increases obviously due to the downward induced velocity caused by rotor wake is neglected and the blade angle of attack is enlarged.
Fig.7 Chordwise pressure distributions of Caradonna-Tung model rotor
Fig.8 Sectional lift coefficient distribution along blade span
The lift coefficient along span is calculated and compared with experiment data in Fig.8. It can be seen from Fig.8, although the simulation of aerodynamics is challenging, the calculated results agree with the experimental data well. Only a little discrepancies happen in blade tip region.
Table 1 shows the comparison of computation time between the coupled method and full CFD method under the same number of blade grid. The computations are performed at the PC with Core i7,3.6 GHz. As shown in Table 1, it takes 20 h to simulate the flow-field for full CFD and2.5 h for the coupled method. These results indicate that the coupled CFD/free-wake method reduces the computational time dramatically.
Table1Comparisonofcomputationaltimebetweenthetwomethods
Number of gridFull CFD methodCoupled CFD/free-wake methodNumber ofblade grid239×34×63239×34×63Number of back-ground grid189×148×1890Computationtime20 h2.5 h
In this section, the well-known OLS rotor in low-speed descending flight condition is used to validate the capability of coupled method in the calculation of BVI air-loadings, which is very challenging in rotor aerodynamic field.The two-bladed OLS rotor is a 1/7 scale model of AH-1 helicopter main rotor[25]. The blade has an aspect ratio of 9.22 and a linear twist of -8.2°. The control inputs are obtained from the built delta trim method and compared with the experimental data[25]and the results from Ref.[26], as shown in Table 2. Here,θ0,θ1candθ1sare the blade collective pitch, longitudinal cyclic pitch and lateral cyclic pitch angles.β0,β1candβ1sare coning, longitudinal flapping and lateral flapping angles.
Table 2 Control inputs of OLS rotor
To indicate the accuracy and efficiency of the proposed trimming method, the histories of control settings and rotor performance with trim cycles by the traditional trimming method and the proposed delta trimming method are shown in Fig.9. And the calculation time is summarized in Table 3.
From Fig.9, it can be seen that the desired aerodynamic values achieve convergence after four to five trim cycles by both methods. Because of the precise coupled model is used to calculate the Jacobian matrix in the traditional trimming method, it converges slightly faster than the delta trimming model from the aspect of the number of trim cycles. But from Table 3, it can be seen that in traditional trimming model the coupled N-S/free-wake solver runs for 16 revolutions, while the solver runs for only 5 revolutions in delta trimming model. In this case, about 68% computation time is saved.
Fig.9 Histories of control settings and aerodynamic levels with trim cycles
Table3Computationaltimeofdifferenttrimmingmethods
TrimmingmethodNumber oftrim cyclesRevolutions ofcoupled solverComputationaltime/hTraditionalDelta451654815.5
Fig.10 Comparisons of blade loadings at several spanwise locations and contour of air-loading
Fig.10 shows the comparison of the predicted blade loads by the proposed method and other three methods (AFDD, DLR, ONERA) from Ref.[27] at 0.91Rand 0.95Rspanwise locations. It is shown that the predicted value by the proposed method is intermediate compared with other results and has good agreements with them. In addition, the proposed method captures four BVI interactions on the advancing side (#1—#4) and two interactions on the retreating side (#5,#6). These results are in accordance with the theoretical analysis ones from Ref.[27]. While other methods (all belong to wake methods) only capture a single strong interaction. This is because a larger time step is necessary in wake model for the computational convergence[12]. As a result, the interactions between time steps are missed.
To illustrate the importance and contribution of this method, the same case is also simulated by the full CFD method. Although the full CFD method has a minor time step (0.25° is used in this calculation), the numerical dissipation and interpolation errors inherent to the difference schemes employed causes the rotor wake to be modelled with less intensity than that occurs in physical reality. This is why the blade-vortex interactions are also not captured.
But, it should also be pointed out that the proposed coupled method has a discrepancy on the calculation of blade-vortex interactions occurring at 270° azimuth. This is caused by the integrated vorticity approach, which induces a smaller velocity on the retreating side. In future work, a distributed vorticity approach, which is more complex one, will be included in the solution to improve prediction accuracy.
The 10 014 cases in OLS rotor far-field noise tests[25]is applied to predict the BVI noise.
The comparison of noise signals at two observer points is conducted. Both observers are 3.44R(Ris the rotor radius) away from the rotor hub, 30° below the rotor plane. Laterally, Observers No.3 and No.9 are located atψ=180° andψ=210°, respectively. Locations of the two observers can be identified in Fig.11.
The comparisons of sound pressure time histories at two different observers among predicted values, experimental data and results from Ref.[27] are shown in Fig.12. It demonstrates that the proposed method has the capability of ca-pturing two peaks of noise signals and a better agreement with the experimental data. But the method shows a little discrepancy with experimental results locally where rotor noise drops after the second peak.
Fig.11 Observation locations in radiation sphere
Besides, the comparison of noise signals by the proposed method and the full CFD method is shown in left-hand side of Fig.12. As shown in Fig.12, due to the numerical dissipation and interpolation errors in the full CFD method, the rotor wakes are not preserved well and thus the BVI noise amplitudes are not predicted precisely.
Fig.12 Comparisons of sound pressure time histories at two observer locations
Sound-pressure level (SPL) contours of the acoustic hemisphere is shown in Fig.13. The acoustic hemisphere, where it has a radius of 3.44Rand the distributed observers on the surface have an interval of 5° in azimuth and elevation direction, is plotted by the Lambert projection method[28]. As shown in Fig.13, the maximum noise level is observed at the 120° azimuth and -45° elevation, which is the preferred direction of BVI noise.
Fig.13 SPL contours of the OLS rotor
A coupled CFD and free-wake method for the prediction of rotor BVI noise is developed with an efficient delta trimming procedure. The method is validated for rotors at hover and forward flight conditions. Based on the FW-H equation, the developed method is used to calculate the rotor BVI noise. The following conclusions can be drawn:
(1) Compared with the field-velocity approach, the outer boundary correction approach is more efficient. The computation time of Jacobin matrix is reduced by the delta trimming procedure due to the efficiency of free-wake model.
(2) From the results regarding the distribution of induced inflow and rotor surface pressure distribution, it can be seen that the proposed coupled method can trim the rotor and predict the BVI aerodynamics accurately.
(3) The fluctuations of air-loadings caused by BVI interaction are fairly well captured with the current coupled method. Compared with wake method and full CFD method, more BVI interactions are captured by the proposed method and the predicted BVI noise by the proposed method agrees with the experimental data.
(4) Rotor BVI noise has strong directivity and it is mainly radiated to 120° azimuth and -45° elevation. This makes the BVI noise more audible to observers on the ground.
Transactions of Nanjing University of Aeronautics and Astronautics2018年5期