Computational aeroacoustics of the EAA benchmark case of an axial fan

This contribution benchmarks the aeroacoustic workflow of the perturbed convective wave equation and the Ffowcs Williams and Hawkings analogy in Farassat’s 1A version for a low-pressure axial fan. Thereby, we focus on the turbulence modeling of the flow simulation and mesh convergence concerning the complete aeroacoustic workflow. During the validation, good agreement has been found with the efficiency, the wall pressure sensor signals, and the mean velocity profiles in the duct. The analysis of the source term structures shows a strong correlation to the sound pressure spectrum. Finally, both acoustic sound propagation models are compared to the measured sound field data.


Introduction
Modern sound design and sound optimization require robust prediction of aerodynamics and aeroacoustics using numerical methods. Therefore, benchmarking computational schemes is of outermost importance and requires sound data [1,2]. The used generic fan benchmark case has been recently issued by the European Acoustics Association under the name "A Benchmark Case for Aerodynamics and Aeroacoustics of a Low-Pressure Axial Fan" [3]. Experimental data is available on the platform, including the aerodynamic performance of the fan, fluid dynamic quantities, and sound measurements. Furthermore, the geometry of the fan and the main dimensions of the experimental setup are available as an IGS (initial graphics specification) as well as a Parasolid file.
Considering low pressure axial fans, main contributions to the current state of the art have been made by several authors [4][5][6][7][8] but are not limited to them. In this contribution, we aim to present the whole aeroacoustic simulation process for a successful computation of the flow and acoustic signature of an axial fan: (1) meshing approach for the computational fluid dynamics (CFD) simulation; (2) turbulence modeling; (3) CFD convergence study; (4) evaluation of the most significant flow results for a subsequent acoustic simulation; (5) computational domain and meshing for the simulation of the acoustic field; (6) acoustic source term evaluation, possible truncation and interpolation [9] from the CFD to the acoustic grid; (7) acoustic field computation and the analysis of its most relevant physical quantities. For the aeroacoustic modeling, we investigate two sound propagation models. The integral equation according to Ffowcs Williams and Hawkings analogy (FWH) in Farassat's 1A version [10] and the finite element formulation of the perturbed convective wave equation (PCWE) [11].
We begin with a short review of the experimental results that are available for comparison. Then, we explain the CFD simulation setup and the obtained CFD results. Next, we discuss the used aeroacoustic formulations and show the results of the acoustic simulation. Finally, we give concluding remarks concerning the results and possible improvement of the benchmark.

Fan test case and experimental results
The experimental setup for the application is presented in [12]. The paper provides a fan in a short duct with an extensive amount of measurement data including aerodynamic performance (volume flow rate, pressure rise and efficiency), wall pressure fluctuations in the duct (see Fig. 1), fluid mechanical quantities (velocity in three spatial direction and turbulent kinetic energy) on the fan suction and pressure side (see Fig. 1). The wall pressure fluctuation measurements consist of data from 15 transducers that were installed 15 mm after the end of the nozzle, with a spacing of 10 mm. Additionally, the acoustic spectra at different microphone positions upstream of the fan is provided.
The fan was designed with the blade element theory for low solidity fans. In terms of size and operating conditions, it is a typical fan to be used in industrial applications. The *Corresponding author: stefan.schoder@tugraz.at reproduction of the geometry has been a main goal. Therefore, the fan was designed with zero blade skew and the design was not optimized for fluid dynamic or acoustic behavior. The blades consist of NACA 4510 profiles [13]. Table 1 shows the design parameters. The circumferential velocity at the blade tip corresponds to a Mach number of Ma % 0.113. Thus, the flow can be considered incompressible. The Reynolds number based on the chord length is almost constant over the span-wise direction of the blade and indicates a turbulent flow [14].
The measurements were made in a standardized anechoic inlet test chamber (see Fig. 2) according to ISO 5801 [15].
The volumetric flow _ V was adjusted by butterfly dampers and an auxiliary fan in the inlet section. The flow field was rectified in the first half of the inlet chamber by a flow straightener. The duct was installed in the chamber wall, with the suction side facing inwards and the pressure side facing outwards. An electrical motor outside of the measurement chamber drove the fan. For more information about measurement setup, we refer to [16]. The measured pressure rise of the fan at the design point is Dp = 126.5 Pa. However, the design pressure difference is not completely reached due to unconsidered losses like tip flow. At the design point the efficiency according to, with M T being the torque of the shaft and n the rotational speed, is g = 53%.
The available measurement results contain microphone signals with a sampling frequency of 48 kHz, and measurement length of T = 30 s. Inside the test chamber, seven (N = 7) microphones were installed horizontally in a halfcircle with a radius of 1 m in front of the nozzle of the duct.
The sound power level was computed according to [17], with the time averaged sound pressure level " L P , the hull of the measurement area S 1 = 6.28 m 2 and S 0 = 1 m 2 . The time averaged sound power level for all microphones is computed as, with the reference pressure p 0 = 20 lPa. For a frequency range of 100 Hz-10 kHz, the measured sound power level was L W = 87.3 dB. The spectrum of the sound power level is shown in Figure 3. The first blade passing frequency (BPF) can be seen as a sharp peak at 223 Hz, the first higher harmonic of the BPF at 446 Hz and the second higher harmonic of the BPF at 675 Hz. Broader peaks are around 340 Hz and 500 Hz that exceed the ones from the BPF. This first and second visible sub-harmonic peak are expected to result from the interaction of the tip flow with the blades [18,19]. Above 800 Hz, the spectrum consists of broadband noise.

Flow simulation
The CFD domain was derived from the computer aided design (CAD) model of the measurement setup. The fan installed in the duct is shown in Figure 4a. The geometry of the fan used in the CFD simulations is shown in Figure 4b. The geometry was simplified to adapt the CAD model to the need of the CFD simulation. The holes in the inlet nozzle as well as the gaps and the front face   Figure 2. Measurement chamber according to [12].
of the hub are slightly modified. The electrical motor outside the measurement chamber is replaced by a cylinder with a similar extent, to account for the blocking in the wake. For the same reason, the struts behind the rotor are preserved in shape. The blade geometry is not altered. The whole simulation domain is shown in Figure 5, with the flow direction from left to right.
The inlet domain has an extent of 2.32 m Â 2.4 m Â 2.4 m which is about the size of the test chamber after the flow straightener. This domain is stationary and contains the nozzle of the duct. The rotating domain is connected with non-conforming interfaces and contains the straight section of the duct, including the fan. Downstream the rotating region, a second stationary region is connected with a further non-conforming interface. Directly after the interface, the diffuser guides the flow into the outlet region. The outlet domain has an extent of 2.0 m Â 2.4 m Â 2.4 m which is the same cross-section as the inlet domain and has a length of 4D with D denoting the duct diameter.
The volume flow rate of the design point results in a mean inlet velocity of 0.24 m/s. Therefore, the influence of the inlet velocity profile was considered negligible and a volume flow rate boundary condition was used. Due to the prescribed flow rate, errors in the numerical simulations result in an error of the pressure rise of the fan. On the outlet surface perpendicular to the main flow direction, a pressure outlet with zero pressure was applied. All other boundaries were modeled as no-slip walls.
Star-CCM+ v.12.06 [20] was used for meshing and the numerical flow simulation. Besides, we performed acoustic computations with the implemented FWH integral.

Turbulence modeling
An incompressible detached eddy simulation (DES), in the IDDES formulation [21], was used to model the turbulent, wall modeled flows. The DES simulation blends between an unsteady RANS (URANS) simulation near the wall and a large eddy simulation (LES), with Smagorinsky's sub-grid scale model, outside the boundary layers. The blending function f DES is shown in Figure 6, where f DES = 0 stands for a complete LES simulation and f DES = 1 for a complete URANS simulation. It can be seen that the URANS model is used just directly at the walls of the duct and the fan and the rest of the domain is treated by LES. Among other things, the blending depends on the mesh size as can be seen at the top, the bottom of Figure 6, as well at the walls, and at the right part of the shaft. For too coarse meshes, the blending function can treat too many areas as URANS simulations. Therefore, the blending has to be checked depending on the simulation.
For the simulation of wall-bounded flows, boundary layers have to be resolved adequately. This yields strong   demands on a mesh to resolve the boundary layer in all regions adequately. The all y + -treatment offers automatic switching between high y + -and low y + -treatment. This option was used to cover regions where the y + 1 criterion is not fulfilled. The y + values of the CFD simulation are shown in Figure 7. The highest y + values occur at the suction side of the blades at the outer radius and where the flow interacts with the duct, with the maximum value of y + = 14.
A further aspect of the simulation of the boundary layers is the used turbulence model. For the URANS part of the DES, two different models are investigated. First, the Spalart-Allmaras (SA) model [22] and second, the k À x shear stress transport (SST) model [23]. The comparison of the time-averaged velocity in the axial direction on the suction side is shown in Figure 8.
On the suction side, the velocity profiles are almost identical. On the pressure side shown in Figure 9, small deviations occur at r/r duct = 0.95, but overall the results are similar. The largest deviations compared to the measurement can be found at 0.65 < r/r duct < 0.9, where the turbulence model does not have a strong effect since this region is modeled by the LES. Due to the higher numerical effort, the k À x SST model is about 10% slower and therefore, the SA model is used.

Discretization and grid convergence study
The CFD simulations were initialized with a RANS simulation. The Reynolds averaged Navier-Stokes (RANS) simulation was run for 10 k iterations. Within this initial calculation, the pressure increase and momentum (and therefore efficiency) of the fan converged. The resulting pressure and velocity fields were used as the initial conditions for the transient simulations. The transient simulations used a time step size of Dt = 10 ls to resolve a frequency up to 5 kHz within the acoustics simulation [11]. This time step size leads to a domain rotation of Da = 0.089°per step. At the outer radius of the rotating domain, this corresponds to a movement of 0.39 mm, which is between 15% and 40% of the used mesh sizes at this location. A second-order implicit time-stepping was used, and each time step was solved with 10 inner iterations. For most regions of the CFD simulations CFL < 1 was fulfilled and the maximum value was CFL max = 150 (with all cells located in the URANS region of the DES).
The incompressible DES is discretized by a cut-cell mesh (trimmed mesh) approach [20]. It uses mainly hexahedral elements, which are subdivided for mesh refinement resulting in hanging nodes. The geometry is resolved by a non-orthogonal cutting of the elements, and leads to partly polyhedral elements. Prism layers are applied to resolve the    boundary layer of the fan and duct. The CFD mesh and the prism layer at the blades are shown in Figure 10.
The rotating domain is discretized with a constant mesh size, which is referred to as the base size h. The nozzle and the diffusor are also discretized with the base size. Therefore, the non-matching interfaces have the same cell size on both sides. Downstream of the nozzle, the mesh size is 4h for 1.8D to resolve the turbulent wake of the fan. Upstream of the nozzle, the same mesh size is applied for 0.5D. Outside of this refinement regions, the mesh gradually increases to a maximum cell size of 50h. For the duct, four prism layers with a total height of 0.5h and a stretching ratio of 1.5 were used; for a base size of h = 2 mm, this leads to an initial height of 6 Â 10 À5 m for the wall-closest element. At the fan blades, five prism layers with a total height of 1/3h and a stretching ratio of 1.6 were used; for a base size of h = 2 mm, this leads to an initial height of 2 Â 10 À5 m. The target surface mesh size on the fan blades was 0.5h. According to the prism layer at the blades and the wall, the tip gap is resolved by twelve elements.
We performed a grid independence study with four different mesh sizes, and use a Richardson extrapolation [24] to estimate the extrapolated solution for the efficiency. From these four simulations, we investigate the mesh convergence. Therefore, a curve-fitting algorithm was applied to extrapolate the simulation results. The cell numbers and the base mesh sizes are displayed in Table 2. Except the prism layers height, the other regions are specified relatively to the base mesh sizes. So the refinement is mostly uniform in space.
The results of the convergence investigation are shown for the initial RANS simulations, as well as for the DES simulations. The results of the DES simulations were time-averaged and are displayed over the cell number n À 2 3 , which is commonly used as a non-dimensional cell size to the power of two [25]. Figure 11 shows the total-to-static efficiency g according to (1), using the simulated pressure difference and the simulated shaft power. The convergence is monotone and the extrapolated result of the RANS simulation underestimates the efficiency by 8.3% and the DES simulation overestimates the efficiency by 3.4%. The deviation of the efficiency is caused by the different inflow conditions, which may increase the losses and therefore explain the overestimation of the DES simulation. A rigorous comparison of the efficiency would also need a distribution of the efficiency measurement. Furthermore, in [26] convergence of the time-averaged axial velocity profile between the hub and the duct was investigated. The coarse simulation deviates from the other simulation results directly at the duct and for r/r duct > 0.8. The rest of the simulations deviate only slightly directly at the hub and the duct wall. At both sides, the measurement result is not completely met, but the influence of the mesh seems very small between the fine and finest mesh. Therefore, no improvement is expected from further mesh refinement.
This common practice of using time-averaged signals for assessing grid convergence is relevant for the simulation setup in the first place. Additionally, fluctuating quantities that are connected to aeroacoustic sources should be investigated. If assessing convergence of computational aeroacoustics simulations, it is essential to study timedependent flow results and the acoustic propagation. Using the final acoustic propagation to evaluate grid convergence, the whole CAA workflow for at least three meshes has to be done. At the current point of simulation time, this is infeasible in a real-world application. Therefore, flow quantities connected to the aeroacoustic source terms are evaluated spectrally and used as a grid convergence criterion of the simulation to reduce the computational effort. Regarding the aeroacoustic models used, the fluctuating pressure values are the primary flow quantity computing the acoustic propagation. For low-pressure axial fans and using FWH, the surface pressure fluctuations at the blades are known to be important aeroacoustic source terms; using PCWE, the CAA source term is the substantial derivative of the   Figure 11. Convergence of the total-to-static efficiency g for RANS simulations (blue) and DES simulations (orange).
incompressible pressure. The pressure measurements close to the leading edge tip gap at transducer 7 was used to assess grid convergence. At this location close to the blade tip, high pressure fluctuations and large acoustic source terms are expected. The comparison of the PSD (power spectral density) from the different simulations is shown in Figure 12. The frequency resolution is coarse due to the short evaluation time of 0.02 s. Until 1 kHz all simulations are similar and represent the BPF and its first harmonic. The coarse simulation overestimates the PSD in the high frequency range, which is partly attributed to the coarse mesh and the modeled turbulence of the filtering function of the LES part. This trend can also be observed for the "middle" mesh, whereas the tendency of the increase is shifted towards higher frequencies. The fine and the very fine simulations reproduce the decline of the spectrum between 4 and 10 kHz quite well. Between 1 and 4 kHz the simulations underestimate the PSD. To conclude, the slope and the peaks of the measurement were reproduced by the fine mesh and the fine mesh is chosen to be a compromise between accuracy and computational cost. Finally, we assess the ratio of modeled to resolved turbulent kinetic energy inside the LES region, regarding the fine mesh. The flow off the walls is modeled by LES (for f DES = 0, see Fig. 6). For too coarse meshes, too less eddies are resolved and too much turbulence is modeled by the sub grid model. The relation of resolved turbulent kinetic energy k res to the modeled turbulent kinetic energy of the sub grid model k sgs is a criterion to estimate the quality of the resolution used in the form, and corresponds to Pope's criterion M = 1 À LES res , estimating the unresolved turbulent kinetic energy [27]. For LES res , this gives a criterion between 0 and 1. For higher values (commonly used value is 0.8 or higher) most turbulent kinetic energy is resolved. Since k sgs is not provided directly in DES simulations, a criterion according to [28] was used. This criterion uses a relation of viscosity l and turbulent viscosity l turb in the form, with a = 0.05 and n = 0.53 providing similar values as (4). For the DES simulation, this criterion is just valid in regions where the LES model is used. Therefore, the values in the boundary layer should be neglected. The result of (5) is shown in Figure 13. Upstream of the fan, over 85% of the turbulent kinetic energy is resolved. In the vortex core resulting from the tip flow, it is about 70%. In the diffusor, it is between 80% and 95%. In most regions of the wake, it is between 70% and 80%. And in the outer shear layer of the wake, especially towards the coarse cells at the outlet, the values drop to a range between 57% and 70%. To conclude, in the region of interest, which is near the fan blades, the fine mesh resolves the turbulent kinetic energy regarding the 80% criteria. As a consequence, we use the fine mesh for all further investigations. Considering the fine mesh, the simulation was run for 0.1 s, or 2.48 revolutions of the fan (respectively 10 k time steps) to obtain a stationary operating point in the solution. After that, every second time step was exported; a total of 0.4 s or 9.9 revolutions (respectively 20 k time steps) were exported for the aeroacoustic simulations.
The computations were performed on the Vienna Scientific Cluster (VSC-3) [29]. The statistics for the CFD simulation are summarized in Table 3, using the final computational setup that is discussed in this section.
The CFD simulation to export the aeroacoustic sources took 225.5 h and 32.1 TB of data was exported.
The results of the CFD simulation are validated by measurements [12]. The geometry of the fan and duct and the evaluation areas are displayed in Figure 1. The two simulation evaluation planes of the flow velocity are displayed in green. They are directly before and after the root of the fan blades. The transducer locations are used as pressure probe locations to validate the CFD simulation. As referenced later, the pressure probe location 7 coincides with the location of transducer 7. Transducer 7 is 75 mm behind the nozzle and directly in front of the blade tip. Major sound sources are expected at this location.

Flow results
In Figure 14, a slice of the unsteady velocity field is shown. Upstream of the fan, the velocity field is relatively smooth, which is the reason why less refinement is necessary. However, this smoothness of the incoming flow field will reduce the generated sound due to the reduced interaction of the flow with the blades leading edge, and we expect no tonal component at the BPF inside the acoustic spectra. The largest velocity amplitudes occur directly at the blade tips. An increased velocity towards the lower duct wall can be seen, which results from the tip flow. At the hub, turbulent structures can be seen, which arise from the horseshoe vortices at the blade roots (see the Q-criterion [30] in Fig. 15). Downstream of the fan, the flow is characterized by large scale eddies that are convected through the duct. The turbulent structures are dissipated as the mesh size increases towards the outlet.
In Figure 16, the wall pressure fluctuations are shown from the measurement and the CFD simulation for the selected transducers (see Fig. 1). Each transducer is either mounted before the blade, inside the tip gap, or after the fan blades and is characteristic for the aeroacoustic sources of the tip gap. The CFD simulation time for the evaluation of the spectral results was 0.1 s.
Upstream of the fan at transducer position 2, the spectrum is dominated by a peak at the BPF. The BPF peak is met by the simulation, but the broadband component above 1 kHz is underestimated. For the positions in the middle of the duct, the overall level is higher and the higher harmonics of the BPF are visible. For transducer 7 the higher harmonics are met, but for transducer 9 the harmonics of order 4 and higher are overestimated. For both positions, the broadband component is slightly overestimated. The highest wall pressure levels occur at transducer 9. Downstream of the fan, the first BPF is similar in amplitude to the upstream spectrum but the broadband component is higher over the whole spectrum. This discrepancy can be explained by the smoothed measurement signal and by modeled turbulent structures of the sub-grid filtering.
In Figure 17 the velocity results in the axial direction on the two measurement planes are displayed, both the measurement results from Laser Doppler Anemometry and CFD simulation results are compared. On the suction side, the velocity is underestimated by the CFD simulation close to the wall, where the overall amplitude and shape of the velocity profile are reproduced. On the pressure side, the velocity is higher close to the wall but lower towards the middle of the duct. On the pressure side, the wake profile from the fan blades is reproduced in shape.

Acoustics -PCWE
The acoustic wave propagation was solved in the time domain using the PCWE, In (6) / a denotes the scalar acoustic potential, c 0 the speed of sound and D/Dt = o/ot + " u Á r the substantial derivative with the mean flow velocity " u. The aeroacoustic source term of this equation is the substantial derivative of the incompressible flow pressure p ic . The PCWE [31,32] is an exact reformulation of APE2 [33] equations as a convective wave equation. For rotating regions, the second term of the     substantial derivative is corrected by the rotational velocity u r of the mesh, with the mean flow " u. This approach corresponds to an Arbitrary Lagrangian Eulerian (ALE) formulation (see [34]). For low Mach numbers, the mean flow effects are small and are neglected in the acoustic simulations. The acoustic pressure can be derived from the acoustic potential as a post-processing result by PCWE has a reduced computational effort compared to the acoustic perturbation equations since it is a single scalar equation. Therefore, the computational operations to solve this equation are fewer, and the needed memory for the system matrices is smaller. Additionally, the scalar source term reduces the amount of CFD results that need to be stored. The PCWE is solved by the finite element method (FEM), which is available in the research software Coupled Field Systems (CFS++) [35]. This method was already successfully applied to different rotating systems [11,36].

Acoustic domain
For the application to fan noise prediction, the geometry of the whole setup has to be represented in the acoustic domain. This means that not only the source domain has to be matched correctly, but also scattering obstacles in the propagation domain. Therefore, the used acoustic domain contains all scattering objects like the fan, duct, nozzle, shaft, strut, and motor. The upstream anechoic chamber is modeled by a free-field condition, due to the high absorption of the walls. The downstream box is also modeled by a free-field condition, assuming negligible sound backscattering of surrounding laboratory equipment to the microphone positions inside the anechoic chamber. An overview of the acoustic domain is shown in Figure 18. The rotating region (in red) is exactly the same as in the CFD simulation. The rotating region and the stationary regions are connected by Nitsche-type mortaring [11]. The inlet region is cropped to still include all the microphones of the measurement (in black) and the outlet region is dimensioned to contain the motor.
The microphones are arranged in front of the nozzle on a half-circle with a radius of 1 m. The whole domain is surrounded by a perfectly matched layer (PML) [37] displayed in green to satisfy the free radiation property of the measurement chamber.

Mesh discretization
The accuracy of the numerical solution depends on the used mesh discretization and the used polynomial order of the finite elements. For all acoustic simulations, first-order polynomial basis functions are used. The finer the mesh becomes, the more accurate is the solution. In contrast to the CFD mesh, where near-wall regions need mesh refinement, the acoustic mesh needs a uniform mesh size to preserve the acoustic waves in the whole computational domain. To reduce errors for finite elements with first-order polynomial basis functions, the smallest acoustic wavelength of interest is discretized with 10-20 points in space and the shortest period is discretized with 10-20 steps in time [35]. From the measured spectrum in Figure 3, the tonal components of the fan can be seen under 1 kHz, where above broadband noise occurs. In this investigation, three different mesh resolutions were used in the acoustic propagation simulations. For the finest resolution of the propagation domain, a spatial resolution of 15 linear elements was chosen for a frequency up to f max = 1500 Hz. The smallest wavelength computes as, and with a standard isentropic speed of sound c 0 = 343 m/s and the maximum frequency the wavelength computes as k min = 0.228 " 6 m. Therefore, a maximum mesh size of h % 0.015 m was used. To investigate the effect of coarser meshes a second mesh with a maximum frequency of f max = 1000 Hz was chosen, which leads to a maximum mesh size of h % 0.023 m (15 elements per wavelength). A third mesh for the same maximum frequency and just 10 elements, which leads to h % 0.034 m was chosen. For harmonic simulations an estimation for the needed discretization exists. A relation between the element order q, the wave number k and the spatial resolution h is given in [38] as, The constant can be assumed as c 1 = 1. For a given element order, the necessary mesh size h can be computed, at which the numerical solution starts to monotonically converge to the true solution. According to the used mesh sizes, this inequality gives a maximum frequency for the fine mesh  Table 4. The numerical accuracy and efficiency depend not only on the resolution but also on the used element type. In the FEM, the hexahedral elements are best in terms of accuracy and efficiency. As a consequence, the large inlet and outlet regions are meshed with pure hexahedral elements. These regions are mainly for the propagation and have little to no contribution of sound-producing source terms. Inside the inlet region, it is possible to create a purely blockstructured mesh. But inside the outlet region, the complex geometry of the strut leads to highly skewed elements. Therefore, a combination of a block-structured mesh and swept mesh blocks is used. Pure quadrilateral surface element meshes were generated. These surface elements are then extruded to a three-dimensional block consisting of hexahedral elements. The acoustic meshes were generated with ANSYS ICEM CFD 18.0.
The rotating region contains the fan, which is complex geometry and difficult to mesh structured. In this region, the ability of the FEM to use different kinds of elements is beneficial, and this region is meshed with tetrahedral elements. This leads to an unstructured mesh with a larger number of volume elements compared to hexahedral volume elements. But it decreases the time used for meshing and improves the overall mesh quality. The acoustic mesh is shown in Figure 19, with the inlet region on the left, the rotating region in the middle, and the outlet region on the right.
Inside the source region, the needed mesh size is not known a priori. Since in addition to the wave propagation resolution, acoustic source terms must be resolved. As shown in [9], all sources in the dimension of the mesh size can be resolved, properly. So the first guess is to discretize this region with the same resolution as the propagation region and to refine where small scale source terms are assumed. This is mainly in the tip gap and closes to the boundary of the fan. Afterward, mesh refinement is done similarly, as it has been performed for the mesh convergence study in CFD. The results of this investigation will be discussed in Section 4.5.
A relevant aspect of the meshing process is to avoid different extents of defeaturing of the underlying geometry. The CFD mesh is created by using CAD geometry. Since meshing defeatures small details, we use the CFD surface and volume mesh geometry to mesh the acoustic domain in a second step.

Aeroacoustic sources PCWE
In numerical simulations, the computational domain has to be finite. This leads to an introduction of artificial boundaries that truncate the physical field. The truncation of the aeroacoustic source changes the acoustic solution, which increases the computational error. For an abrupt truncation of the source domain, some correction techniques have been derived [39]. In some cases, the truncation of acoustic sources is desirable to reduce the impact of nonphysical artifact sources. For example, in the far-field, no acoustic source terms should occur. For this purpose, it is common to use blending functions (see, e.g. [40,41]) to suppress un-physical sources, because a smooth blending function reduces the un-physical sources in the blending area and leads to less un-physical acoustic radiation.
In this work, we use a blending function for the aeroacoustic sources (see Fig. 20). It is designed to utilize the sources in the rotating domain and drop smoothly to cancel out the sources before the non-conforming interfaces. This should reduce the un-physical sources, which occur due to numerical errors at the interfaces.
A similar truncation of the simulation occurs in time, where the abrupt starting of the simulation can be seen as a multiplication of the sources with a Heaviside function in time. This can lead to high acoustic impulses, propagating through the computational domain. In this case, for the very first 0.5 ms (25 acoustic time steps) all source terms are suppressed to remove any numerical artifacts from the start of the simulations, e.g., time derivatives with not enough time steps for the complete computational stencil. In the second 0.5 ms, the sources are increased C 1 -smoothly to their actual value for the remaining simulation time.

Acoustic directivity of the duct
The acoustic result in the measurement chamber depends on the radiated sound of the fan and the interaction with the surrounding geometry. The geometry influences the acoustic directivity of a source by reflections and diffractions of the propagating waves. For the acoustic wave, the wavelength k is the characteristic parameter. For the geometry, the duct diameter D is a reasonable choice. As long as the wavelength is larger than the duct diameter,  Figure 19. Acoustic mesh, with pure hexahedral meshes in the inlet and outlet region (gray) and a tethrahedral mesh in the rotating region (red).
the interaction of the acoustic propagation with the geometry has a small effect. If the wavelength and the duct diameter are of similar size, modifications of the directivity are expected. The directivity was computed using a time-harmonic simulation on the fine acoustic mesh, where an artificial excitation was applied on a single node on the rotation axis in front of the fan. The source has a monopole behavior. In the case of free-field radiation, the directivity would be uniform. The amplitude was evaluated at the microphone positions of the measurement setup. The result is displayed in Figure 21 for different relations of the wavelength and the duct diameter. The directivity inside the measurement chamber is displayed in the left half of the diagram, with the rotational axis of the fan at the horizontal line. It can be seen that for lower frequencies the radiation is almost omnidirectional. For wavelengths much larger than the duct diameter, like for D = 0.15k in Figure 21, the acoustic radiation is almost not influenced by the duct. From the measurements, it is known that in this range of frequencies, the tonal components of the acoustic spectrum occur. Therefore, the influence on the directivity is supposed to be small for the tonal components. When the wavelength comes in the range of the duct diameter (700 Hz), the acoustic radiation starts to be influenced by the duct's nozzle geometry and the directivity towards the rotational axis increases. From the measurements, it is known that above this frequency the acoustic spectrum contains mostly broadband noise. Therefore, the broadband noise is supposed to be influenced by the geometry. With increasing frequency and therefore shorter wavelength side lobes start to form. For D = 1.47k (1000 Hz), this leads to a directivity, which has its main lobes 40°off the rotational axis. For D = 2.35 (1600 Hz), a strong directivity towards the rotational axis forms, where the side lobes are about 10 dB smaller than the main lobe. The nozzle has an amplifying effect on the result since it is almost formed like an exponential horn [42].

Influence of mesh discretization
The mesh discretization influences the spatial resolution of the simulation, the acoustic result, and the computational effort. In the acoustic propagation simulation, different acoustic mesh discretizations have been investigated as described in Section 4.2. The acoustic results from the different meshes at microphone position 4 are shown in Figure 22. Overall, the different acoustic meshes lead to very similar results. The results in the frequency range below 1 kHz are expected, since all meshes were designed to resolve this range well. It also means that the mesh inside the source region is accurate enough to resolve the spatial characteristics of the sources. In the higher frequency range between 2 and 5 kHz, small deviations occur in the acoustic results where the results of the fine mesh are slightly closer to the measurement. In the frequency range above 5 kHz the simulation on the coarse mesh is strongly influenced by the time-stepping method (Hilber -Hughes -Taylor time-stepping scheme, see [43]). Since this mesh does not resolve high frequency waves, the time-stepping algorithm strongly damps these waves, so that the numerical solution is not deteriorated in the low frequency range [34]. From a simulation point of view, only the fine mesh is capable of resolving the acoustic propagation in the high frequency range accurately.
The statistics for the acoustic propagation simulations with different mesh sizes are shown in Table 5. The propagation simulations were again performed on single nodes of the VSC3 with 16 CPUs. The computational effort in simulation time and data export rises with the total number of elements. For the coarse mesh, the simulation time is in the order of the interpolation time, but for larger meshes the simulation time exceeds the interpolation time and for the fine mesh the solution time is even higher than the source computation time on the CFD mesh.

Analysis of source terms
More detailed information about the actual sound sources can be gained by the sources of the PCWE equation. At an instantaneous time, iso-surfaces are shown in Figure 23. Sources are at the leading edges of the blades, constantly distributed from the root to the tip, the leading edge noise sources are responsible for the tonal components at the BPF [44]. However, due to the smooth inflow  conditions, the noise sources at the leading edge are not predicted that dominantly. Smaller sources occur at the outer radius of the duct, covering the whole circumference. They are expected to result from turbulent flow structures of the tip flow. Tiny noise sources occur at the trailing edges of the blades, which are supposed to result from vortex shedding. To discuss the frequency content and spatial extent of the sources, the aeroacoustic sources were transformed into the frequency domain. A moving reference frame was used to perform the Fourier transform inside the rotating domain. Selected results are shown in Figure 24. The isosurfaces have the same values in all figures and show instantaneous temporal animation of the time-harmonic signal. The short simulation time restricts the frequency resolution of the Fourier transformation. Therefore, the closest frequency to the BPF of 223 Hz was 220 Hz, and the Fourier transform of the sources at the first sub-harmonic peak is shown. At the BPF, sources are located at the leading edge, which are stronger than at the sub-harmonic peak frequency but do not dominate. At the outer radius, long stretched sources occur in the circumferential direction. For the first sub-harmonic peak, the acoustic sources are weaker at the leading edge and at the hub. At the outer radius, the acoustic sources are more compact, due to the higher frequency.
From the CFD simulations, information can be obtained to observe possible sound sources of the transient flow. One of the main sound sources for this fan is the interaction of the tip flow with the following blades. In Figure 25, the instantaneous streamlines through the tip gap of the blade pointing towards the observer are shown at the last time step of the data exportation. These streamlines reach the blade tip of the following blade, which leads to interactions with the blade and generates noise. Some of the streamlines pass this blade and reach even the next two following blades. This interaction of the tip flow following blades is expected to cause the sub-harmonic hump [19]. But some of the streamlines are redirected downstream in the diffusor and do not interact with the following blades. For larger flow rates, it was recognized that the more streamlines are redirected downstream and this interaction is reduced. Figure 15 shows iso-surfaces of the Q-criterion colored with the flow velocity. These vortical structures are known to be sound sources [45]. Large structures can be seen to originate from the blade tips. They coincide with the large acoustic sources of the PCWE equation. Furthermore, at the hub between the blades, vortices can be seen that arise from the secondary flow (see, e.g. [46]). They can be expected to be responsible for the sources between the blades in Figure 24. In addition to that, small horseshoe vortices at the roots of the blades are visible, but the Q-criterion gives no information about the acoustic sources at the leading edges of the blades. This indicates that the leading edge noise is not resolved using a uniform inflow condition and hence the peaks at the blade passing frequency are not dominant. Using CFD simulation, particular flow phenomena can be investigated that result in acoustic sources. However, only the assessment of aeroacoustic sources in connection with an acoustic propagation model detects relevant flow-induced sound sources.

Acoustics -FWH
The application of the original FWH analogy yields some challenges using incompressible flow data. The most significant one is that the acoustic fluctuations have to be propagated correctly until the enveloping evaluationsurface. However, for low Mach number flows, the volume terms can be neglected and a simple surface formulation can be used [18]. Farassat's formulation 1A [10,47] gives the total fluctuating pressure from the surface as a contribution from the thickness and the loading of the surface, The thickness term is given as, with the distance from source to receiver r = |x À y|. In (12) Ma r denotes the Mach number in direction of the receiver, the dot quantities are derivatives with respect to the source time and the integrands are evaluated at the retarded time. The vector U n is defined as, The loading term is given as, with the blade loading vector, which is applied in the direction of the receiver L r = L i Á r i and the direction of the surface movement L M = L i Á v i /|v|. The terms of order 1/r are considered the far-field terms and the terms of order 1/r 2 are considered the near-field terms. For transient rotating motion Farassat's formulation 1A is implemented in StarCCM+ for impermeable integration surfaces only [20]. This means, that the integration surfaces have to coincide with solid surfaces and the fluid velocity and normal velocity are equal, and therefore the equations simplify with u n = v n . Furthermore, for low Mach number flows, the flow velocity is considered incompressible. With this simplification used in the investigation, the density reduced to a constant. Theoretically, the constant density leads to an infinite speed of sound. Thus, the pressure fluctuations propagate instantly through the whole domain. This means that in the surface formulation all 1/c 0 terms vanish. Therefore, the acoustic results are only valid for compact sources on the fan blades, being compact for f < 2c 0 /(D À d) = 686/(0.495 À 0.248) Hz = 2777.32 Hz. For this example, the integration surface coincides with the surface of the fan blades and neglects the volume contributions [18].

Sound propagation results
In this section, the results of the PCWE equation and FWH integral solution of the complete simulation are shown at the different microphone positions (see Fig. 18). The transient CFD results were used after 5 revolutions of the rotor and data applied to the acoustic computations cover 7.5 revolutions (respectively 0.3 s). For the PCWE equation, the acoustic wave propagation was computed on the coarse mesh. Since the microphone positions 5-7 are arranged symmetrically to the positions 3-1, only the results at the locations 1-4 are displayed in Figure 26. The spectrum of the result at microphone position 1 is shown in Figure 26a, which is the outermost position. The low frequency range is underestimated by both prediction methods. This difference can be explained partly by the lower frequency limit and a box mode of the test rig close to 60 Hz. Whereas, we assume free field radiation for the simulations. Additionally, the PCWE results have a slightly lower amplitude at the lower frequencies as the FWH. This can be explained by the used aeroacoustic sources, which are the time derivative of the incompressible pressure. The missing peaks at the blade passing frequency are   connected to the smooth inflow field used for the simulation. As already discussed, no vortical structures are present inside the undisturbed inflow. This leads to fewer interactions at the leading edge and a reduced peak at the blade passing frequency [48]. The first sub-harmonic peak is underestimated by PCWE, whereas the FWH analogy meets the measurement quite well. The second subharmonic peak is met by both approaches. From 1 kHz to 2.5 kHz the two spectra are similar with an additional peak at 1.7 kHz only for PCWE. Above 2.5 kHz, the PCWE result drops sharply, where the FWH result stays at the same level but gets noisier. The sharp drop of the sound signal can be explained by the time-stepping scheme, which damps unresolved frequencies. The noisy characteristics of the FWH results may be due to the violation of the compactness assumption. Similar behavior, but at a higher amplitude, can be observed for microphone position 2 ( Fig. 26b) and 3 (Fig. 26c), which are closer to the axis of rotation. For these two positions, the low frequencies are predicted better by the FWH analogy. Whereas the PCWE is a bit closer to the measurements in the higher frequency range. Microphone position 2 shows the same peak at 1.7 kHz as position 1 and an additional peak at 4.2 kHz. For position 3, these peaks are not visible. Microphone position 4 is directly on the axis of rotation (Fig. 26d). The PCWE result is slightly overestimating the frequency range between 500 Hz and 5 kHz and the FWH analogy is underestimating up to a frequency range of 3 kHz after which the spectrum is dominated by noise. At this location, the largest differences between the two aeroacoustic methods occur. This difference is explained by the different assumptions regarding the two aeroacoustic formulations. Firstly, the FWH underestimates the signal due to the free-field radiation assumptions and neglecting the duct directivity at the frequency above 1 kHz. Secondly, the amplitudes of the PSD for the PCWE results are stronger at the higher frequencies because it considers the directivity of the duct. Summarizing it can be said that the FWH analogy yields better results for the first sub-harmonic peak and gives reasonable results up to a frequency of 2-3 kHz. Above that, the results are partly distorted because of the simplifications of the used setup of the surface formulation.  In the higher frequency range, this analogy is underestimating the PSD. The largest difference to the measurement occurs directly at the axis of rotation. This might be due to the influence of the surrounding geometry. At low frequencies, this influence is small, but at high frequencies, the amplitudes directly in front of the nozzle increase. The PCWE results are closer to the measurements in the high frequency range, but underestimating the first subharmonic peak. Since this is estimated better with the FWH analogy, it can be assumed that it is not a problem of the underlying CFD simulation. Furthermore, due to the high resolution, it is unlikely that it is a problem of the acoustic propagation simulation. For the first subharmonic peak, the wavelength is k % 1 m and even with the coarse discretization of h = 0.034 m this means that the smallest wavelength is discretized with 29 linear elements, which results in a high resolution. Since the FWH analogy is incorporated in the flow solver, it needs little additional effort to obtain these results compared to PCWE. However, it should be mentioned that the used FWH solver does not take the surrounding geometry into account and does not enable an investigation of the acoustic source mechanisms. The outermost microphone positions are hardest to predict in the high frequency range for both analogies.
In addition to the spectral analysis of both prediction methods, the time signal can be used to compute the over all sound power level using (2). The sound power level gives an insight how good the prediction agrees with the measurement in an energetic sense (see Tab. 6). The PCWE method has a deviation of only 0.6 dB compared to the measurements and is therefore closer to the measurement than the FWH analogy. The deviation of 2.1 dB might be a result from to the used simplifications described in Section 5.

Conclusions
The perturbed convective wave equation (PCWE) and the Ffowcs Williams and Hawkings (FWH) analogy have been successfully applied to the EAA benchmark case of a low-pressure axial fan. Computational fluid dynamics and computational acoustic mesh convergence studies have been performed for the underlying meshes by global and local quantities. Flow quantities, such as the pressure and velocity field adjacent to rotating fan and blades show good agreement to measured data. Based on the test rig description, an aeroacoustic setup has been developed and the hybrid workflow using PCWE has been applied. The duct, where the fan is installed, interacts with the radiated acoustic field using the PCWE. At frequencies lower than 700 Hz, the characteristic is nearly uniform, but with increasing frequency, the directivity gets more focused. This directivity behavior explains the discrepancy of the acoustic sound radiation considering the aeroacoustic simulations at the microphone position 4. The hybrid workflow applying PCWE has been compared to an FWH far-field prediction method resulting in similar trends and results. Concerning the benchmark case, it should be considered to provide additional data regarding the inflow velocity profile and the sampling distributions for global performance parameters (like the efficiency).

Conflict of interest
Author declared no conflict of interests.