Aeroacoustic prediction based on large‑eddy simulation and the Ffowcs Williams–Hawkings equation

A hybrid noise computation method is presented in this paper. Large-eddy simulation with wall-model equation is proposed to compute the flow field. With a stress-bal-anced wall-model equation, the near-wall computation cost of large eddy simulation was effectively reduced. The instantaneous flow variables obtained by the large-eddy simulation were used to compute the noise source terms of the Ffowcs Williams– Hawkings equation. The present method was investigated with two test cases: a single cylinder at Re = 10,000 and a rod-airfoil at Re = 480,000. The flow quantities and aeroacoustic characteristics were compared with the reference data. The mean velocity profiles and spectra of the flow fluctuations were consistent with data from the litera-ture. When compared with the reference data, the noise computation error was less than 3 dB. The computation results demonstrate the present wall-modeled large eddy simulation is efficient for the noise computation of complex vortex shedding flows.

turbulence structures and sound wave propagation. The sound field and flow field are two-way coupled in the direct numerical simulation. However, the computational cost of the DNS is expensive for high Reynolds number (Re) flows. According to Choi and Moin [8], the DNS grid requirement is proportional to Re 37/14 , which is not feasible for engineering flows. The hybrid noise computation method computes the near-field sound source and the far-field sound propagation separately. The sound source is computed by an unsteady CFD method, such as large eddy simulation (LES) [9] and detached eddy simulation (DES) [10]. The sound propagation is computed by an acoustic analogy equation, such as the Ffowcs Williams-Hawkings (FW-H) equation [3,11]. In the hybrid method, the flow field and sound field are one-way coupled; that is, the sound field cannot influence the flow field. The grid requirement of a wall resolved LES is proportional to Re 13/7 [8], which is too expensive for realistic aircraft flows. The grid requirements of the wall-modeled large eddy simulation (WMLES) [12] and detached eddy simulation (DES) methods are considerably lower than those of the wall-resolved LES method, which is commonly applied in aeroacoustic prediction.
The FW-H equation [13] is a sound analogy method that is widely used for engineering flows. The sound source computed by the unsteady CFD method is integrated as the source term of the wave equation. The source terms of the FW-H equation contain three types: monopoles, dipoles, and quadrupoles [14]. The monopole source can be considered a pulsating sphere that scales as the 4th power of the flow velocity [15] and can result from unsteady volumetric flow addition. The dipole source is caused by unsteady momentum fluxes that scale as the 6th power of the flow velocity [15]. The unsteady pressure load caused by vortex shedding can be thought of a dipole source. The quadrupole source is the most common type of turbulent flow and is caused by collisions of fluid elements. The quadrupole source scales as the 8th power of the flow velocity [15], which is less efficient than the monopole and dipole in low-speed flows.
In this paper, aeroacoustic predictions were carried out with LES and the FW-H equation. Two test cases were used to test the effectiveness of the methods. The first case was a flow around a circular cylinder at Re = 10,000 [16], which was computed by a wallresolved LES. The second case was an interaction flow between a circular cylinder and an airfoil [17], which is a benchmark problem for airframe noise computation (BANC) of the American Institute of Aeronautics and Astronautics (AIAA). The Reynolds number was 480,000 based on the airfoil chord length. This case was computed by a wall-modeled LES method. The flow field and sound emission of the two cases were compared with data from the literature. The results demonstrate that the present hybrid method with wall-modeled LES and FW-H equation has high reliability and broad engineering applicability.

Wall-modeled large eddy simulation
The present study is based on a compressible finite volume solver [18]. The convection and viscous terms of the Navier-Stokes equations were both discretized with a second-order numerical scheme. A central scheme and an upwind scheme were combined to form the spatial scheme. The blending parameter was determined by the grid quality [18], which makes the scheme less dissipative in the good grid quality regions and more dissipative in the poor grid quality regions. The third-order Runge-Kutta method was used for time marching. Vreman's [19] model was applied for subgrid stress closure.
The wall-stress model of the boundary layer, which is an important feature for reducing the computational cost, was implemented, as shown in the sketch in Fig. 1. The LES grid was quite coarse near the wall (usually the first grid layer Δy+ > 10). A wall-model velocity profile was embedded at the bottom of the LES grid. The apex velocity of the wall-model profile was obtained from the LES grid at the interpolation point. Kawai et al. [20] suggested that the interpolation point should be in the log-layer, and there should be more than two LES cells under the interpolation point. A one-dimensional wall-model grid with high grid density was automatically generated to solve the wall-model equation.
In the wall model framework, the viscous effect was only considered in the y-direction. In the very near-wall region, the pressure was assumed to be independent of the inner layer, and a quasi-steady state was assumed in the bottom of the boundary layer; consequently, the boundary layer equation becomes a "stress balance model" [21] or an equilibrium model [12].
As shown in Eq. (1), the equilibrium model is a one-dimensional ordinary differential equation (ODE). In the framework of a wall-modeled LES, the shear stress τ w , as shown in Eq. (2), can be computed by the wall model and used as the wall boundary The wall boundary condition for Eq. (1) is a nonslip condition, and the apex velocity is interpolated from the instantaneous flow field of the LES. To solve the ODE, approximately 40 to 60 grid points with a stretching ratio of less than 1.05 are typically automatically generated. The Δy+ value of the first wall-model cell is less than 1.0.

Ffowcs Williams-Hawkings equation
The far-field noise was assessed by the Ffowcs Williams-Hawkings (FW-H) equation [13]. Given a stationary sound source in a fluid moving at a constant speed U = (U 1 , U 2 , U 3 ), the convective FW-H equation for a permeable surface S is given by Eq. (5) [11,22]. Q n , F i and T ij are shown in Eq. (6). The subscript 0 denotes the free stream quantity, and c 0 is the sound speed. The superscript ' represents the perturbation relative to the free stream; for example, ρ ′ = ρ − ρ 0 . According to the proposal of Spalart and Shur [3], the sound perturbation density ρ ′ can be expressed as p ′ /c 2 0 . H(S) and δ(S) are the Heaviside function and Dirac delta function, respectively; ⌢ n = ⌢ n 1 , ⌢ n 2 , ⌢ n 3 is the unit normal vector pointing outward from the surface S; u i are the local flow velocities on S; δ ij is the Kronecker delta; and σ ij is the viscous stress tensor, which is often negligible in far-field sound computation.
By defining an observer location x = (x 1 , x 2 , x 3 ) and a source location y = (y 1 , y 2 , y 3 ), the three-dimensional free-space Green's function for a convective wave equation can be expressed as Eq. (9). The definitions for R * , R and x i , i = 1, 2, 3 are shown in Eqs. (10), (11), and (12), respectively. M 0 = U 2 1 + U 2 2 + U 2 3 /c 0 is the Mach number. α and φ are the angle of attack and the sideslip angle, respectively, which are defined by With Green's function, a solution for Eq. (8) in subsonic flow (M 0 < 1) can be obtained, as shown in Eq. (13).
If S is a solid surface, then the source terms Q n and F i become Eq. (14). T ij is the quadrupole term. If a permeable integration surface S contains a turbulent flow region, the quadrupole source must be considered. However, the quadrupole source is negligible for low-speed flow, as implied by the solid surface integration.
3 Flow around a circular cylinder at Re = 10,000

Computational grid
In this section, the numerical methods were validated by a circular cylinder. The freestream Mach number of the cylinder was 0.20, and the Reynolds number was 10,000. A wall-resolved LES simulation was carried out because the Reynolds number was low and the boundary layer was laminar. FW-H integration based on a solid wall was used to calculate the far-field noise. The spanwise length of πD was computed for this case, where D is the cylinder diameter. A 12.6 million point mesh with 49 grid points in the spanwise direction and 258,000 points on each x-y plane was used. The computation domain was Figure 2 shows the computational grid in the x-y plane. As shown in Fig. 2a, the grid has a C-H-O topology. An O-grid block was placed around the cylinder to ensure good grid orthogonality near the wall. An H-grid block was placed after the cylinder to better capture the turbulent wake. A C-grid block surrounded the O-grid and H-grid blocks. Figure 2b shows the grid near the cylinder. A total of 353 grid points were located along the circumferential direction of the cylinder. The grid space of the first layer normal to the wall was 5.0D × 10 − 4 , which ensured that Δy+ < 1.
To test the grid resolution, Δx+ and Δy+ were collected and averaged after the LES computation. Figure 3 shows the Δx+ and Δy+ of the first grid layer normal to the wall. The maximum Δx+ was less than 50, and the maximum Δy+ was less than 1.0, satisfying the requirements of the wall-resolved large eddy simulation [23] in both the x and y directions.

Flow field results
The total computation included approximately 400 time units (tU ∞ /D). The last 200 time units were used for statistical averaging. The nondimensional timestep ΔtU ∞ /D was 1.82 × 10 − 4 . Approximately 2.2 million iteration steps were used in the total computation. Figure 4 shows the lift and drag coefficient histories after statistical collection. The time histories of the lift and drag coefficients demonstrate periodic oscillation with variational amplitudes. The statistical force coefficients are compared with DNS data [16] in Table 1. The present results, including the average drag coefficient and the root-meansquares of the lift and drag coefficients, were all within the error ranges of the reference compressible DNS data. The Strouhal number St of the vortex shedding was obtained based on the periodic oscillations of the force coefficients and is presented in Table 1. This result was also similar to the reference data.
The instantaneous flow fields of the cylinder at four typical timesteps are presented in Fig. 5. The curves in the figures correspond to the lift coefficients, and the symbols        [16]. The reference compressible data are a direct numerical simulation computed by a 6th order solver that were used to directly assess the far-field sound [16]. Both the averaged velocities and the RMS of the present computation were consistent with the reference data, even at far downstream locations (x/D = 10). A slight discrepancy appeared in the v RMS near the y = 0 location at x/D = 5 and 10. The grid topology and numerical scheme were quite different from those of Khalighi [16], which might result in different behaviors at far downstream locations. The comparison demonstrates that the present computation is reliable and satisfactory.

Sound results
The sound perturbation density ρ ′ = ρ − ρ 0 is important for sound emission. The contour of the instantaneous sound perturbation density is shown in Fig. 8. The acoustic wave around the cylinder is clearly captured. Additionally, some small waves are associated with downstream vortex shedding.
The accuracy of the turbulence fluctuation is important for aeroacoustic computations. To determine the computational accuracy of the turbulence fluctuation, two probes were used for comparison with the reference data [16]. The probes were located at y = 0, with one probe at x = 5D and the other at x = 15D. The power spectral density of the velocity v was compared with the compressible and incompressible DNS results of Khalighi [16], as shown in Fig. 9. The incompressible result of Khalighi [16] has lower energy at the high-frequency part than the compressible result because the 6th order compressible solver used by Khalighi has a higher resolution than the 2nd order incompressible solver [16]. The present result was consistent with the reference compressible DNS data, demonstrating the good temporal-spatial resolution of the present computation method. The main peak also matched well; however, the second peak of fD/u 0 St = 3 was slightly overpredicted. Figure 10 shows the sound spectra of the present computation and Khalighi [16] at an observer location of x = − 1.2D, y = 16.2D. The sound computation was based on the solid surface integration of the FW-H equation, which means the integration surface was the wall of the cylinder; consequently, only the dipole source was considered in Eq. (13). The frequency in the figure was normalized by the vortex shedding Strouhal number. The power spectrum density (PSD) was nondimensionalized. The main peak and the spectrum slope of the present computation agree well with those of Khalighi [16]. A minor flaw in the present computation is that the third harmonic was slightly overestimated, which corresponds with the near-field spectra in Fig. 9a. The result demonstrates that the primary flow patterns and sound propagation mechanism are well captured by the present computation.

Case description and computational grid
The rod-airfoil interaction case is a benchmark noise case of the BANC workshop of the AIAA. Abundant experimental data were provided by Jacob et al. [17]. The case includes an upstream circular cylinder and a downstream airfoil, as shown in Fig. 11a. The cylinder has a diameter of 0.1C (C is the chord length of the airfoil), with the center located at x = 0 and y = − 0.02C. The airfoil was a NACA0012, and the leading edge was located at x = 1.05C and y = 0. The incoming flow Mach number was 0.21, and the Reynolds number based on the airfoil chord was 4.8 ×  Figure 11b and c show the grid near the cylinder and the airfoil. An O-type grid was generated near the cylinder and airfoil to ensure grid orthogonality.
The height of the first grid layer was 1.5C × 10 − 4 . In this computation, the wall shear stress was computed by the wall-model equation. Consequently, the first grid layer was not required to be Δy+ < 1. The Δy+ of the first layer was determined after the computation. Figure 12 shows the Δy+ of the cylinder and airfoil. The maximum Δy+ on the cylinder was approximately 7, while it was approximately 4 on the airfoil. Because the y-coordinate of the upstream cylinder was slightly lower than that of the airfoil, the cylinder wake did not impact the airfoil symmetrically; consequently, the Δy+ of the upper and lower surfaces of the airfoil were not equal.

Flow field result
The nondimensional timestep ΔtU ∞ /C has a value of 3.15 × 10 − 5 . Approximately 40 time units (tU ∞ /C) were computed, and the last 25 time units were used to collect the averaged flow field quantities. Figure 13 shows the three-dimensional and two-dimensional instantaneous flow fields of the rod-airfoil case. It is clear that the Karman vortex-street after the circular cylinder is well captured and impacts the airfoil. The Karman vortexstreet is highly three-dimensional because the Reynolds number based on the cylinder diameter was in the subcritical region, and the flow became turbulent in the cylinder wake. Several points are marked in Fig. 13b for comparison with the experimental data. The coordinates of the points are listed in Table 2.
The flow fields were time-averaged and collected for comparison with the experimental data [17]. Figure 14 shows the time-averaged velocity profiles at several streamwise locations. The x/C = 0.01 location is above the cylinder boundary layer, before the  The velocity profile at this location matched well with the experiment. The x/C = 0.18 location is in the cylinder wake. The velocity deficit peak at this location was higher than the experimental data, implying that the turbulence mixing was underpredicted in the computation. The velocity deficit peak matched better at x/C = 0.795, indicating that the turbulence mixing was well captured at this location. The velocity profiles above the airfoil (x/C = 1.3 and 1.9) and in the airfoil wake (x/C = 2.15) matched well with the experimental data, which illustrates that the present computation is reliable.
The velocity spectra of the sampling points marked in Fig. 13b were compared with the experimental data [17] in Fig. 15. The first three points (A, B, and C) were in the circular cylinder wake, and the remaining three points (D, E, and F) were near the airfoil. The velocity spectra profiles all matched well with the experimental data. The vortex shedding frequency was clearly visible at the first five sampling points, and the peak became indistinct at the trailing edge of the airfoil (point F). The computed vortex shedding frequency, which was normalized by the airfoil chord length, was approximately fC/U ∞ = 1.85. The vortex shedding frequency in the experiment was slightly higher than that in the computation. The relative error was less than 5%. The amplitudes of the velocity spectra, even at low frequencies, were well captured by the present computation, demonstrating that the flow fluctuations were accurately simulated by the numerical method and providing satisfactory input for the aeroacoustic computations.

Aeroacoustic result
The FW-H equation was adopted to compute the aeroacoustic characteristics of the rod-airfoil test case. Both the permeable integration surface and the solid-wall surface integration were compared in this computation. Figure 16 shows the computational surfaces for this case. The solid-wall surfaces of the cylinder and the airfoil can be integrated by the FW-H equation together or separately. The purple rectangular surface was used as the permeable integration surface. Closing disk averaging [25] was used on the planes on the right side of the integration surface to reduce the integration error caused by the interaction between the turbulent structures and the integration surface. A far-field location 18.5C away from the leading edge of the airfoil was adopted to test the aeroacoustic characteristics. The experimental data from Jacob  [17] were used for comparison. The spanwise length in the experiment [17] was 3.0C, while the length in the present computation was 0.1πC. The computed sound pressure energy was multiplied by 3.0/0.1π to match the experiment. Figure 17 shows the power spectrum density of the sound pressure at r = 18.5C and α = 45° (α = 0° is the x-positive direction). The solid-surface integration and permeable-surface integration were nearly identical. A small difference can be observed in the high-frequency region. Due to the low-speed flow, the quadrupole noise source caused by the turbulent structures was not significant in this case. Consequently, only solid-surface integration was used in the following comparison. Figure 18 shows a comparison of the computational and experimental results. The sampling position was located at r = 18.5C and α = 45°. The present computation was consistent with the experimental results. The computation errors of the tone peak and the broadband noise were less than 3 dB. The tone peak frequency was slightly lower than that in the experiment, while the broadband noise was well captured. Figure 19 shows the contribution of each component. The airfoil contributes most of the sound pressure level, while the cylinder has a small contribution. The airfoil generated strong noise because it was located in the turbulent wake of the cylinder. Figure 20 compares the sound pressure spectra at different angles along the circumferential direction. The abscissa is the circumferential angle, and α = 0° represents the x-positive direction. The ordinate is the frequency on a logarithmic scale. The computation spectra map matched well with the experimental data. Both the first peak and the second peak were visible in the spectra map. The result demonstrates that the present numerical method is satisfactory for aeroacoustic computation.

Conclusions
Aeroacoustic predictions were carried out by large eddy simulations and the Ffowcs Williams-Hawkings equation in this paper. The work can be summarized as follows.
(1) The present hybrid method for aeroacoustic computation demonstrates good computation accuracy in both the single-cylinder and rod-airfoil interaction cases. The wall-model method can efficiently reduce the computational cost of the large eddy simulation. The velocity profile of the bottom part of the boundary layer can be solved by a one-dimensional wall-model equation. The flow fields match well with the reference data. The Ffowcs Williams-Hawkings equation demonstrates good aeroacoustic postprocessing accuracy for calculating flow-induced noise. (2) For the rod-airfoil case, permeable surface integration and solid surface integration produce nearly identical results because of the low-speed flow. The airfoil is the main noise source due to its location in the cylinder wake. The vortex-shedding frequency is the main frequency tone of the airfoil.