Computation of Aerodynamic Noise Radiated from Ducted Tail Rotor Using Boundary Element Method

A detailed aerodynamic performance of a ducted tail rotor in hover has been numerically studied using CFD technique.The general governing equations of turbulent flow around ducted tail rotor are given and directly solved by using finite volume discretization and Runge-Kutta time integration. The calculations of the lift characteristics of the ducted tail rotor can be obtained. In order to predict the aerodynamic noise, a hybrid method combining computational aeroacoustic with boundary element method (BEM) has been proposed. The computational steps include the following: firstly, the unsteady flow around rotor is calculated using the CFDmethod to get the noise source information; secondly, the radiate sound pressure is calculated using the acoustic analogy Curle equation in the frequency domain; lastly, the scattering effect of the duct wall on the propagation of the sound wave is presented using an acoustic thin-body BEM.The aerodynamic results and the calculated sound pressure levels are compared with the known technique for validation.The sound pressure directivity and scattering effect are shown to demonstrate the validity and applicability of the method.


Introduction
With wide applications of helicopter, its noise level has become one of substantial indexes of helicopter performance and acceptability [1][2][3][4][5].Helicopter tail rotor is a major noise source; therefore, investigation of tail rotor noise characteristic becomes distinguishably valuable and important to control the helicopter noise level.Tail rotor is an important aerodynamic assembly of conventional configuration helicopter, and its aerodynamic characteristics could have a significant impact on performance and flying qualities of the helicopter.For this reason the analysis of the tail rotor has been one of the research priorities in the field of helicopter technology [6][7][8].
Due to the effect of the duct, the slipstream pattern and aerodynamic characteristic of ducted tail rotor are different from the conventional open-type rotor.It can reduce the risk of component damage and enhance the operational safety.In addition, both the duct and rotor of the ducted tail rotor can produce thrust.Thrust which is produced by the duct is due to negative pressure on the duct inlet, which contributes maximum 50% to the total.Furthermore, the noise of ducted tail rotor is lower than the conventional open-type rotor [9,10].Consequently, the analysis of the aerodynamic noise of the ducted tail rotor is the main topic for this paper.
In order to predict the aerodynamic noise, aerodynamic characteristics of the ducted tail rotor should be obtained firstly.Along with the development of the ducted tail rotor, many researchers have studied the aerodynamic characteristics of various types of ducted tail rotors [11][12][13].In recent years, the main method which investigates the aerodynamic performance of ducted tail rotor is CFD (Computational Fluid Dynamic) technique.Among those CFD techniques, momentum source method is one of the commonly used methods.The momentum source with functional relationship to the local flow conditions represents the spanning rotor in flow field and is added to the N-S equations by UDF procedure.Cao and Yu [14] used it to solve the numerical simulation of turbulent flow around ducted tail rotor.Song et al. [15] applied it to analyze the aerodynamic characteristics of ducted tail rotor.Even though momentum source method is quite successful for the prediction of the aerodynamic performance of ducted tail rotor, detailed three-dimensional flow features cannot be accurately predicted.Because of the Mathematical Problems in Engineering approximation of the blade shape parameters, the accuracy of the momentum source method is still limited.In addition, it is difficult to obtain the accurate noise source information by using the momentum source method.In this work, a CFD technique which is considered the blade shape parameters has been proposed to investigate the aerodynamic characteristics of ducted tail rotor, and the method is based on moving reference frame (MRF) model instead of momentum source method.
Calculation of the acoustic field radiated by rotating sources is a meaningful problem in the prediction of noise of aircraft rotors.The main techniques for noise prediction are based on acoustic analogy Curle equation, FW-H (Ffowcs Williams-Hawkings) equation, and the generalized treatment of Goldstein [16][17][18].The inhomogeneous acoustic wave equation divides the aeroacoustic source into three types: monopole source, dipole source, and quadrupole source.
Although the aerodynamic characteristics of ducted tail rotor have been discussed by many researchers [19,20], the investigations of the aerodynamic noise are relatively few.The prediction of rotating blade aerodynamic noise can be acquired easily by solving acoustic analogy Curle equation, but the scattering effect of the duct wall on the propagation of the sound wave is difficult to discuss and the scattering mechanism is rather complex.The acoustic boundary element method (BEM) is usually applied to predict the sound radiating and scattering field in the exterior and interior closed domain [21,22].Hu et al. [23] proposed a thin-body BEM, in which an imaginary interface is constructed to divide the domain into interior and exterior subdomains, and the imaginary surface is not discretized in the numerical implementation.In the present study, the thin-body boundary element method is also used to consider the scattering effect of the solid wall and to predict the far-field aerodynamic noise of the ducted tail rotor.

Governing Equations
Because the rotary speed is constant, in order to study the noise of the ducted tail rotor, we focus on the flow field in stable state.So the unsteady Navier-Stokes equation is applied here to calculate the flow field.In rotating coordinate system, the unsteady Navier-Stokes equation can be expressed as follows [24]: where  is control body,  is the boundary area, and  is the normal vector. is conservation vector, and  = [,  0 , V 0 ,  0 , ]  ,  is the density of the fluid,  0 , V 0 ,  0 are projection of absolute velocity in relative coordinate, and  is the internal energy of the unit fluid.  and   are the flux tensors.  is source term.
Equation ( 1) is discretized by using Jameson cell-centered finite volume method.An explicit time integration algorithm based on the multistep Runge-Kutta is used to advance the solution in time.The linear system obtained by using the above method is solved by applying Gauss-Seidel method.The standard  −  model is adopted here because of its reasonable accuracy and fewer resources.The moving reference frame is also used to solve the problem of rotating.The far-field boundary is treated as the nonreflecting boundary condition.In order to accelerate the convergence, a local time step and implicit residue smoothing have been used.The detailed procedure for solving the governing equations can be found in [25].

Free Field of Rotor Noise
Lighthill aeroacoustic analogy method is used to investigate the generation of the aerodynamic noise.The general equation is given by [16] ( where  0 is the speed of sound and   is the acoustic pressure. The terms on the right-hand side of (2) are interpreted as source terms. denotes the volume source per unit volume per unit time, and / is the monopole source.  denotes the force per unit volume, −  /  means the dipole source.  is the Lighthill stress tensor and  2   /    denotes the quadrupole source.After the CFD calculation, the aerodynamic noise at observation point is calculated using (2) by referring to the time history of the blade loading of the CFD results.Since a low rotational speed is considered for the tail rotor, the monopole source and dipole source are dominant, and quadrupole source could be negligible.A brief description will be given later to explain that the loading noise generated by dipole source is much bigger than the thickness noise generated by the monopole source for the tail rotor.Then (2) can be reduced to The free space solution of (3) is [26] where (, ,  − ) is Green function given by where  = | − | denotes the distance between the observer location  and source location .The relation between observe time  and source time  is  =  − / 0 .The square bracket means evaluation at the retarded time .
Using the equality (4) can be written as By using the Fourier transform and its inverse transform, force   , sound pressure   (, ), and (, , −) can be given in the frequency domain as where  = / 0 is called the wave number and  = √ −1 the complex unit.Substituting ( 10) into (7), we obtain Comparing ( 11) with ( 9), we have Differentiating ( 12) with regard to the direction of the normal unit vector (), we get where (, ) is the static pressure on the source surface in frequency domain.The function (, , ) and its derivations (, , )/() and  2 (, , )/()() in above equations are expressed as follows:

Thin-Body Acoustic Boundary Element Method (BEM)
In this section, a thin-body boundary integral formulation is applied to calculate the far-field sound pressure.Due to the scattering effect of the solid wall in the duct, the total sound pressure is acquired as the sum of the incident and scattered pressure where    and    are incident and scattered sound pressure in the frequency domain, respectively.The incident sound pressure can be obtained by using (12), and the scattered sound pressure will be got by using BEM.The calculated domain is shown in Figure 1; the surface of the duct is denoted by  and an imaginary surface  is added to close the duct and divide the computational domain into an interior subdomain  + and an exterior subdomain  − .The sound pressure on the outside of the surface  +  is denoted by  − and that on the Mathematical Problems in Engineering inside is denoted by  + .The integral equation can be used to each subdomain [23]  + ()  + (, ) =    (, ) where  1 and  2 are normal unit vectors at the two sides of the wall and  + () and  − () are the two constants that depend on the position of Adding (16) gives the thin-body boundary integral equation where / 1 = −/ 2 = /, and the continuous boundary conditions of the pressure and its partial derivation on the imaginary surface  are used The assumption of acoustic rigid boundary conditions are applied over the entire surface Then ( 18) reduce to Equation ( 22) is not sufficient to obtain the two unknowns  + (, ) and  − (, ).Differentiating (21) with regard to the direction of normal vector (), it can be transformed into The problem of scattering by the duct can be dealt with by initially solving (23) to calculate the sound pressure jump  + (, ) −  − (, ) on the surface , and afterwards evaluating the acoustic pressure at any field point.Then the acoustic pressure on both sides of the duct could be easily got.The integral Equations ( 21)-( 23) cannot be directly calculated; a discretized scheme based on BEM should be used to calculate the unknown value  + (, ) −  − (, ).The simplest constant boundary element is applied in this paper and the thin-body surface is discretized into  elements.Each element has one node which is located in the center of the element.Then, (23) can be transformed into a system of algebraic equations.where   =    (  , )/(  ),   =  2 (  ,   , )/ (  )(  ), and   =  + (  , ) −  − (  , );   and   are the center of element , , respectively.  is the area of the element .Therefore, the form of matrix of ( 24) is written as Equation ( 25) may be solved easily by using the software Mathematica.When the unknown  is calculated, the acoustic pressure at any field point can be obtained through (21) and (22).

Geometric Configuration and Mesh Generation.
To demonstrate the effectiveness of the present method, we consider a model of ducted tail rotor at TsAGI, which has been installed on the Ka-60 helicopter.Figure 2 shows CATIA geometric model of TsAGI ducted tail rotor.Its geometric dimensions are tabulated in Table 1 as reported in [27].A structured mesh is generated around the duct of the TsAGI model, and the surface quadrilateral for the all configuration is presented in Figure 3.The tail rotor of the TsAGI model is meshed with an unstructured mesh which is shown in Figure 4, and an adaptive encryption technique is adopted to mesh the leading edge and the trailing edge of the tail rotor.

Aerodynamic Performance of the Ducted Tail
Rotor.An initial condition is made for the ducted rail rotor at a tip Mach number of 0.22 and a collective pitch angle of 40 degree at the blade root.At this moment, the predicted thrust from the present calculation is 9.547 kg.The rotor thrust in [28] and [29] are 9.54 kg and 9.69 kg, respectively.The coincidence indicates that our method can investigate aerodynamic performance of ducted tail rotor effectively.The predicted sectional thrust distribution along the rotor in Figure 5 is also compared with the momentum source theory and the vortex theory [27], which shows a good agreement.
In Figures 6 and 7, spanwise distribution of the induced downwash and the circumferential induced velocity are compared with the experiment and the results obtained by vortex theory [27].It is also shown that the induced downwash and the circumferential induced velocity are generally in good coincidence with the experiment results.The results obtained by our method also compare well with those of the vortex theory.In addition, the present results indicate better comparison than those of the vortex theory when approaching toward the tip.
Figure 8 shows the distribution of the pressure on the surface of the ducted tail rotor in hover.It is obvious that the pressure changes violently on the inlet of the duct while few on the other zones remain the same.
A grid sensitivity analysis is given for a steady-state CFD simulation.The induced downwash, the circumferential induced velocity, and the induced attack angle are analyzed for the hybrid grids which are obtained by steadily increasing the meshing parameters.Their mean error ( mean ) and maximum error ( max ) with respect to the reference grid are established in Table 2.  denotes the max grid cell size.
From Table 2, we can see that the errors for all three flow variables are smaller and smaller when the grid cells increase.The small errors also display a converging trend.
To verify the effectiveness and the accuracy of the CFD method, the analytical solution of a dipole source has been   investigated.A dipole has an axis, which is the line connecting the centers of the two monopoles (pulsating sphere with radius of ) in Figure 9.The two monopoles are positioned in the  direction at a distance of  = 2.5 from the center ( is the distance between the two monopole points).The pressure fluctuation induced by the dipole source for an observer positioned at (, , ) is    where  0 = tan −1 (1/),  1 = √ 2 +  2 + ( + 0.5) 2 , and To perform the flow field characteristic of the dipole source, we put the dipole source into the center of a sphere (the radius   of the sphere to be 3.25).The speed of sound  0 is 340 m/s.The density for medium is 1.2 kg/m 3 .The angular velocity of the source is 1020 rad/s.The other parameters are  = 4 2 ,  = 0.01 m, and  = 8 m/s.The value of the spherical pressure at any point can be acquired by using the CFD method and (26), respectively.
The calculation is performed for this case, and the spherical pressure at any point is calculated by using CFD method and analytic method.The CFD results are compared  with those of the analytical solution in Figure 10.From Figure 10, we can see that the CFD results are in perfect coincidence with the analytical solution.

Aerodynamic Noise of the Ducted Tail
Rotor.The sound pressure is expressed as dB (decibels) and the predicted SPLs (sound pressure levels) are given by the following: where   denotes predicted pressure and   denotes the reference pressure and equals 2 × 10 −5 Pa.Before investigating the aerodynamic noise of the TsAGI ducted tail rotor, we firstly consider the radiation sound field generated by its isolate tail rotor.Figures 11 and 12  noise, respectively.They can be obtained by using Farassat 1A formulation [26].From the two figures, we can find that the thickness noise generated by the monopole source is much smaller than the loading noise generated by the dipole source for the isolate tail rotor.For the ducted tail rotor, the loading noise is dominant in rotation noise.The main reason is that the number of tail rotors is too high, and the size of tail rotors is too small.Therefore, the loading noise generated by the dipole source is considered only for the following discussions.
Through the analysis of aerodynamic performance of TsAGI ducted tail rotor, we can acquire the unsteady flow properties and the enough pressure data of the surface of the rotor.Then they are transformed into data in frequency domain by applying fast Fourier transform method.The farfield sound pressure is calculated using Curle's analogy in the frequency domain.Before using the BEM method to study the radiation and propagation of the sound source, we should use a simplified configuration which is shown in Figure 13 to consider the scattering effect of the duct.From Figure 13, we know that the surface is constructed by quadrilateral mesh.The maximum size of the mesh is 0.023 m and the number of elements is 2584.In this section, we use a point force method which regards the rotors as a series of points to investigate the scattering effect of the duct.The method can be applied in the prediction of the noise when the force fluctuations of the rotors are obtained.In the present article, we just focus on rotation noise of the ducted tail rotor.The dominant noise in rotation noise is blade passing frequency (BPF) noise and its harmonic noise [30].And the BPF of the ducted tail rotor is 2400 × 11/60 = 440 Hz.
Acoustic pressure in frequency domain is calculated by combining Curle equation and the thin-body BEM.For different observation distance, the directivity of SPLs at the BPF and 2BPF are shown in Figures 14 and 15, respectively.Free field is also rotor alone.The more interesting result is the comparison of the total field and the free field.This provides information about the scattering effect of the duct.For the angle (180 to 360 degrees), upstream of the rotor (180 to 270 degrees) is slightly louder due to scattering.And at the angle (270 to 360 degrees), it is quieter.From Figure 14(a), we can see clearly that the sound pressure at inlet is bigger than that of the outlet.From Figures 14(b), 14(c), and 15, we can also find the directivity of the incident sound pressure is symmetrical as the properties of the dipole sound source.These characteristics are consistent with the fact.
For 30 m/s forward flight velocity, the directivity of SPLs at the BPF is shown in Figure 16.The polar plot is made at a distance of ten times the rotor radius (10R) away from the center of the duct.Compared with the SPLs in Figure 14, the SPLs for forward flight are bigger than the SPLs for hover, and the noise directivity is also changed.
The directivity of SPLs at the BPF for rotational speed 2100 rpm and 2400 rpm are displayed in Figure 17.The SPLs of rotational speed 2400 rpm are bigger than that of 2100 rpm; this can be noticed by taking a closer look at Figure 17.The directivity of SPLs at the BPF for various observation distances is presented in Figure 18.From Figure 18, we can find that the SPLs become more and more small when the observation distance increases.In order to illustrate the validation of the scheme, the results obtained in this paper are compared with the results obtained by commercial software Virtual.Lab Acoustics (VLA) [30].The comparisons are shown in Figure 19.
In Figure 19, at the majority of azimuth of the ducted tail rotor, the SPLs obtained by our method are in good agreement with those obtained by using VLA.

Conclusions
The aim of this paper is to develop the Curle/thin-body BEM method for predicting the aerodynamic noise of the ducted tail rotor numerically.A CFD method based on MRF model is used to study the aerodynamic performance of ducted tail rotor.Calculations are made for TsAGI ducted tail rotor, and the results are compared with the known results and the analytical solution of a dipole source for validation.
The unsteady blade force is also obtained and transferred to the Curle equation to calculate the sound source without the duct.In order to consider the scattering effect of the duct, a thin-body BEM is proposed to predict the far-field aerodynamic noise.The calculated SPLs at BPF accord well with the results obtained by VLA and show that our method can effectively predict the far-field SPLs of the ducted tail rotor.The numerical results also show that the scattering effect of the duct should be included in the prediction of the noise.The control of the inlet noise is advisable for the ducted tail rotor noise reduction.Distance between observer and source :

Nomenclature
W a v en u m b e r : Duct surface : Imaginary surface (): N o r m a lu n i tv e c t o ro nt h eo b s e r v e r surface (): Normal unit vector on the source surface   : Predicted pressure   : Reference pressure.

Figure 1 :
Figure 1: A diagram of acoustic scattering by a thin-body.

Figure 3 :
Figure 3: Surface quadrilateral for the duct.

Figure 4 :
Figure 4: Surface quadrilateral for the tail rotor.

Figure 5 :Figure 6 :Figure 7 :Figure 8 :
Figure 5: Comparison of the spanwise distribution of the blade loading.

Figure 9 :
Figure 9: Dipole source and spherical are shown.

Figure 10 :Figure 11 :
Figure 10: The spherical pressure at any point for different method.

Figure 12 :
Figure 12: Directivity of SPLs of the loading noise at different distance.

Figure 13 :
Figure 13: Mesh of duct used in BEM.

Figure 18 :Figure 19 :
Figure 18: Directivity of SPLs of the ducted tail rotor at different distance.

Table 1 :
Geometric dimensions of the ducted tail rotor.

Table 2 :
Mean and maximum values of the induced downwash, velocity, and attack angle for different grid densities.