Numerical method to compute acoustic scattering effect of a moving source

In this paper, the aerodynamic characteristic of a ducted tail rotor in hover has been numerically studied using CFD method. An analytical time domain formulation based on Ffowcs Williams–Hawkings (FW–H) equation is derived for the prediction of the acoustic velocity field and used as Neumann boundary condition on a rigid scattering surface. In order to predict the aerodynamic noise, a hybrid method combing computational aeroacoustics with an acoustic thin-body boundary element method has been proposed. The aerodynamic results and the calculated sound pressure levels (SPLs) are compared with the known method for validation. Simulation results show that the duct can change the value of SPLs and the sound directivity. Compared with the isolate tail rotor, the SPLs of the ducted tail rotor are smaller at certain azimuth.

ducted tail rotors (Rajagopalan and Keys 1997;Bandoh et al. 1998;Vuillet and Morelli 1986). In recent years, the main method which investigates the aerodynamic performance of ducted tail rotor is CFD (Computational Fluid Dynamic) method. Among those CFD methods, 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 are added to the N-S equations by UDF procedure. Cao and Yu (2005) used it to solve the numerical simulation of turbulent flow around ducted tail rotor. The method not only can predict the flow field but also investigate the performance of the ducted tail rotors. Song et al. (2009) applied it to analyze the aerodynamic characteristics of ducted tail rotor. To obtain the aerodynamic characteristics of the ducted tail rotor, the author solved the N-S equation added by the source. The numerical solutions were verified by the experimental data, and some aerodynamic characteristics of the ducted tail rotor were also discussed. 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 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 method 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 methods for noise prediction are based on acoustic analogy Curle equation, FW-H (Ffowcs Williams-Hawkings) equation and the generalized treatment of Goldstein (Williams and Hawkings 1969;Farassat andKenneth 1987, 1997). 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 Yu and Cao 2006), the investigation of the aerodynamic noise is relatively few. The prediction of rotating blade aerodynamic noise can be acquired easily by solving FW-H 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 filed in the exterior and interior closed domain (Ih and Lee 1997;Wu 2000). Wu and Wan (1992) 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. The author indicated the scattering mechanism of the thin body boundary, and proposed the Helmholtz boundary integral equation to analyze the scattering effect.
In the present study, the analytical formulation of the acoustic velocity computation is also derived for sources in arbitrary motion. In order to consider the scattering effect of the thin duct, the newly developed FW-H-Helmholtz boundary elements method is introduced. The derived velocity analytical formulation is used as the Neumann boundary condition for the thin-body BEM. This method can predict the far-field aerodynamic noise of the ducted tail rotor.

Governing equations
Given the constant rotary speed of the ducted tail rotor, the flow field in stable state is only taken into consideration to study its noise problems. 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 (Li and Xiong 2007) where V is control body, ∂V is the boundary area, n is the normal vector. W is conservation vector,and W =[ρ, ρu 0 , ρv 0 , ρw 0 , ρE] T , ρ is the density of the fluid, u 0 , v 0 , w 0 is projection of absolute velocity in relative coordinate, E is the internal energy of the unit fluid. H W and H V are the flux tensors. S V is source term.
The flow field of the ducted tail rotor is acquired via solving the numerical solutions of Eq. (1). The standard k − ɛ 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 non-reflecting boundary condition. The detailed procedure for solving the governing equations can be found in Cao et al. (2003).

Ffowcs Williams-Hawkings equation
In 1969, Williams and Hawkings (1969) used the generalized function theory to derive the sound equation of the control plane in arbitrary motion in static fluid, that is, the famous FW-H equation. The FW-H equation is given by is wave operator, u i is velocity, f denotes a moving Kirchhoff surface, p' is acoustic pressure, v n is normal component of surface velocity. P ij denotes the compressive stress tensor, T ij = −P ′ ij + ρu i u j − c 2 ρ ′ δ ij denotes a component of the Lighthill tensor, δ(f) is Dirac delta function, H(f) is Heaviside function and satisfied As f (x i , t) = 0, according to non penetration condition, u n = v n , FW-H equation can be reduced as below ∂[ρ 0 v n δ(f )] ∂t, −∂ P ij n j δ(f ) ∂x i and ∂ 2 [T ij H (f )] ∂x i ∂x j are monopole source, dipole source and quadruple source, respectively. According to Hanson and Fink's theory, Eq. (3) can be simplified to (4), the formulation of Farassat 1A, is expressed as follows (Farassat and Kenneth 1987) with p ′ T (x, t) is thickness noise and p ′ L (x, t) is loading noise. The subscript ret indicates that all of the values have to be taken at the retarded time. The dot over the quantity denotes the differentiation of this magnitude with respect to the emission time.

Acoustic velocity formulation for sources in arbitrary motion
The procedure of the velocity formulation for the thickness and loading sources has been recently proposed by (Farassat 2007). Following the same procedure gives the thickness and loading acoustic velocity as follows: Ti and a ′ Li are the acoustic velocity components for the thickness and loading sources. u n and v n are the fluid and the data surface normal velocity, respectively. ρ denotes the local fluid density, and ρ 0 is the density of the undisturbed medium. The subscript ret * indicates that all of the values have to be taken at the retarded time t * .
For any function F(x, τ(x, t)), one has where Ṁ r =Ṁ iri . Then the acoustic velocity components for the thickness sources are got Next, simplifying Eq. (12) Fig. 1 to connect the relationship between sections "Governing equations", "Radiated sound field" and "Acoustic velocity formulation for sources in arbitrary motion". The flowchart also explains the whole numerical procedure of the total sound field.

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 P ′ I and P ′ S are incident and scattered sound pressure in the frequency domain, respectively. Based on Eqs. (6), (7), the incident sound pressure can be obtained and was transformed into data in frequency domain via using Fast Fourier Transform (FFT) method. The scattered sound pressure will be got by using BEM. The calculated domain is shown in Fig. 2, the surface of the duct is denoted as S, an imaginary surface sis added to close the duct and divide the computational domain into an interior subdomain D + and an exterior subdomain D − . The sound pressure on the outside of the surface S + s is denoted by P ′− and that on the inside is denoted by P ′+ . The integral equation can be used to each subdomain (Wu and Wan 1992) where n 1 and n 2 are normal unit vectors at the two sides of the wall, C + (x) and C − (x) are the two constants that depend on the position of x ∂n 2 (y) G(x, y, ω) − P ′− (y, ω) ∂G(x, y, ω) ∂n 2 (y) dS(y) Fig. 1 A flowchart of the total sound field calculation

Fig. 2 A diagram of acoustic scattering by a thin-body
When x ∊ S, C + (x) = C − (x) = 1 2 . Therefore, adding Eqs. (27) and (28) gives the thinbody boundary integral equation When x / ∈ S, adding Eqs. (27) and (28), then where ∂/∂n 1 = −∂/∂n 2 = ∂/∂n, and the continuous boundary conditions of the pressure and its partial derivation on the imaginary surface s are used The assumption of acoustic rigid boundary conditions are applied over the entire surface S Then Eqs. (31) and (32) reduce to Equation (36) is not sufficient to obtain the two unknowns P ′+ (x, ω) and P ′− (x, ω). Differentiating Eq. (35) with regard to the direction of normal vectorn(x), it can be transformed into The problem of scattering by the duct can be dealt with, by initially solving Eq. (37) to calculate the sound pressure jump P ′+ (y, ω) − P ′− (y, ω) on the surface S, and afterwards evaluating the acoustic pressure at any filed point. Then the acoustic pressure on both sides of the duct could be easily got. The value of ∂n(x) cannot be obtained easily via using Eqs. (6) or (7), but can be acquired it indirectly by using the acoustic velocity formulation. If Eq. (37) satisfies the Neumann boundary condition, then for the loading sources. Then, to solve the Eq. (37), a discretized scheme based on BEM should be used to calculate the unknown value P ′+ (y, ω) − P ′− (y, ω). The simplest constant boundary element is applied in this paper and the thin-body surface S is discretized into N elements. Each element has one node which is located in the center of the element. Then, Eq. (37) can be transformed into a system of algebraic equations.
∂n(y j )∂n(x i ) and ϕ j = P ′+ (y j , ω) − P ′− (y j , ω), x i and y j are the center of element i, j, respectively. dS j is the area of the element j. Therefore, the form of matrix of Eq. (39) is written as Equation (40) may be solved easily by using the software Mathematica. When the unknown ϕ is calculated, the acoustic pressure at any filed point can be obtained through Eqs. (35) and (36).

Numerical results and discussions
In this section, the sound pressure is expressed as dB (decibels) and the predicted SPLs (sound pressure levels) is given by the following where P e denotes predicted pressure, P r denotes the reference pressure and equals to 2 × 10 −5 Pa.

Pulsating sphere
In order to verify the algorithms, the analytical solution of a monopole source has been investigated. The monopole is identified with a pulsating sphere as the small sphere with a radius a in Fig. 3. The pressure fluctuation induced by the pulsating sphere is expressed by a harmonic spherical wave where ω and k are angular velocity and the wave number, respectively, and The velocity is given by where A = 4πa 2 U, the radius of the spherical penetrable data surface r s equals to 3.25a. The speed of sound c 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 a = 0.01 m and U = 8 m/s. The pulsating sphere is located at the center of the duct, which is shown in Fig. 4. The diameter of the duct is 0.07 m. The length of the duct is 0.5 m. The observer distance is assumed to be 1 m.
The thin-body BEM is used to perform the acoustic scattering problems of the duct. And the acoustic velocity is calculated by using Eq. (20) and used as Neumann boundary condition for this problem. Figure 5 shows the scattering performance of the pulsating sphere. The left is the incident sound pressure, the middle is the scattering effect of the duct and the right is the total sound pressure.
From the Fig. 4, it is possible to see that the directivity of the incident sound pressure is circular as the property of monopole sound source. When the scattering effect is considered, the directivity of the sum sound becomes non-circular. For the angle (210°-330° and 30°-150°), SPLs of the total field is louder due to scattering where at the angle (−30° to 30° and 150° to 210°), it is quieter. It shows that the sound pressure is strengthened in the direction of duct both ends, and the sound pressure is subdued in the other direction.

Geometric configuration and mesh generation
The model of TsAGI ducted tail rotor, as shown in Fig. 6, is used to demonstrate the effectiveness of the present method. It has been installed on the Ka-60 helicopter. Its geometric dimensions are tabulated in Table 1 as reported in the Bourtsev and Selemenev (2000).
A structured mesh is generated around the duct of the TsAGI model, and the surface quadrilateral for the all configuration is presented in Fig. 7. The number of the structured grid cell is 153 k. The tail rotor of the TsAGI model is meshed with an unstructured mesh which is shown in Fig. 8, and an adaptive encryption method is adopted to mesh the leading edge and the trailing edge of the tail rotor. During the calculation, the tail rotor is involved in a rotating area which is meshed with an unstructured grid (See Fig. 9). The number of the unstructured grid cell is 343 k. So the total number of the computational grid is 496 k. Figure 10 shows the overall mesh of the ducted 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° at the blade root. At this moment, the predicted thrust from the present calculation is 9.547 kg. The rotor thrust in the Tan et al. (2008) and Lee and Kwon (2003) are 9.54 kg and 9.69 kg respectively. The coincidence indicates that our method can investigate aerodynamic performance of ducted tail rotor effectively.
In Figs. 11, 12, spanwise distribution of the induced downwash and the circumferential induced velocity are compared with the experiment and the results obtained by vortex theory (Wu 2000). It is also shown that the induced downwash and the circumferential induced velocity are generally in good coherence with the experiment results. The results obtained from our method also compare well with that of the vortex theory. Figure 13 displays the pressure variation on the surface of the ducted tail rotor. As can be seen in Fig. 13, the pressure changes violently on the surface of the rotor and the duct lip, and others remain the same. Furthermore, the conclusion can be reached from Fig. 13: the thrust of ducted tail rotor is generated mainly by the rotors and the duct lip.
A grid sensitivity analysis is shown for the CFD simulation to illustrate the convergence of the method. The induced downwash and the circumferential induced velocity are discussed for the hybrid grids by increasing the meshing parameters. Their mean error (e mean ) and maximum error (e max ) are shown in Table 2, where D denotes the max grid cell size.   Table 2, it is possible to find that the mean and maximum errors for the flow variables are smaller and smaller with the grid cells increasing. The small errors also show a converging trend.

Aerodynamic noise of the ducted tail rotor
The unsteady flow properties and the enough pressure data of the surface of the rotor can be acquired via the analysis of aerodynamic performance of TsAGI ducted tail rotor. Then they are transformed into data in frequency domain by applying fast Fourier transform method. The far-filed sound pressure is calculated using FW-H equation in the frequency domain. The acoustic velocity is calculated by using Eqs. (20), (25) and used as Neumann boundary condition. Before using the BEM method to study the radiation and propagation of the sound source, a simplified configuration, as shown in Fig. 14, should be used to consider the scattering effect of the duct. In this section, for a low speed rotor application, it is shown that dipole source is dominant and monopole source and quadrupole source shown to be negligible.
Acoustic pressure in frequency domain is calculated by combing FW-H equation and the thin-body BEM. For different observation distance, the directivity of SPLs at 440 Hz and 880 Hz are shown in Figs. 15 and 16 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°-360°), upstream of the rotor (180°-270°) is slightly louder due to scattering where at the angle (270°-360°), it is quieter. From the Fig. 15a, it is easy to see that the sound pressure at inlet is bigger than that of the outlet. From the Figs. 15b, c and 16, it is possible to find  In order to illustrate the validation of the scheme, the results obtained in this work are compared with the results obtained by commercial software Virtual. Lab. Acoustics (VLA) (Zhan and Xu 2013). VLA is commercial simulation software which used specially to noise analysis. It can approximately analyze the aerodynamic noise and vibration noise by using boundary element method and finite element method. The comparisons are shown in 17.
In Fig. 17, the SPLs obtained by our method are in good coincidence with those obtained by using VLA.

Conclusions
The purpose of this work is to propose the FW-H/thin-body BEM method for predicting the aerodynamic noise of the ducted tail rotor in hover numerically. Based on the CFD method, aerodynamic characteristic of the TsAGI ducted tail rotor are discussed, and the results are compared with vortex theory and momentum source method for validation. The numerical method is convergent by increasing the number of the grid. The thrust of the duct is generated by duct lip.
A thin-body BEM is developed to predict the far-field aerodynamic noise of the ducted tail rotor. An analytical formulation is derived for prediction the acoustic velocity generated by moving bodies. It can be used as boundary condition for the thin-body BEM. The scattering effect of the pulsating sphere and ducted tail rotor are investigated by using the thin-body BEM. The calculated SPLs accord well with the results obtained by VLA, and show that our method effectively predict the farfield SPLs of the ducted tail rotor. Numerical examples show that the duct can change the value of SPLs and the sound directivity both the pulsating sphere and tail rotor.