Aeroaocustic Numerical Analysis of the Vehicle Model

: Understanding local phenomena connected with airﬂow around road vehicles allows to reduce the negative impact of transportation on the environment. This paper presents using numerical tools for Computational Fluid Dynamics (CFD) and Computational AeroAcoustic (CAA) calculation. As a model for simulation, simpliﬁed car geometry is used, which is known in the research community as an Ahmed body. The study is divided into two main parts: a validation process and a CAA analysis using the Ffowcs Williams–Hawkings (FW-H) analogy. Research is performed using k − ω Shear Stress Transport (SST) and the Large Eddy Simulation (LES) turbulence model. To compare results with other authors’ studies, three different comparison criteria are introduced: a drag coefﬁcient for different velocities, characteristic ﬂow structure, and velocity proﬁles. The CAA analysis is presented using colormaps and Fast Fourier Transformation (FFT). The methods used in this work allow visualizing the acoustic ﬁeld around reference geometry and determining the frequency range for which the A-weighted sound pressure level is the highest.


Introduction
Many countries around the world have liberalized regulations on the limit values of equivalent noise level A regarding traffic noise. However, the growing interest in electric cars is going to reverse this trend. More and more often combustion engines are replaced with electric ones, which eliminates noise from the fuel combustion process. At the same time, research is conducted to minimize the tire noise. New tire designs and silent road surfaces are created. In this situation aerodynamic noise begins to play an important role in the overall noise generated by vehicles. Based on Helfer's research [1], it can be observed that the aerodynamic noise begins to play a similar role to the tire noise at 120 km/h for passenger cars and 80 km/h for van vehicles. For this reason, it is important to develop techniques to study these types of phenomena.
The numerical calculations are generally applicable in Mechanical Engineering. The results are received without the need to prepare high-cost measurements, where in many situations professional equipment is needed. One example is the study of aerodynamic and aeroacoustic parameters of vehicles. In that case a wind tunnel with measuring equipment is necessary. Additionally, it is required to isolate the object of research from unwanted interaction with the measuring environment, such as noise from fans or sound reflection.
Computational Fluid Dynamics (CFD) responds to all of these requirements. CFD is a part of fluid mechanics that uses numerical methods to find solution for partial differential equations and approximate flow fields of fluids even for complex physical systems. For this purpose the finite-volume method is used. The issue is to prepare an appropriate mathematical and numerical model that reflects the phenomena of such a case. Many works in the automotive area present investigations into aerodynamic parameters like drag and lift coefficients or pressure distribution. This is particularly important due to the amount of fuel consumed and the generated aerodynamic downforce. More complex models are related to multi-physical phenomena. Shi et al. [2] developed a numerical model of exhaust gas propagation for two passenger cars. Their paper presents the effect of turbulence caused by vehicles on the concentration of pollutants. Another example of multi-physical phenomena is presented byŻółtowski et al. [3]. Their research focused on the analysis of structure vibrations under the influence of air flow load (a truck driving under a footbridge at high speed). Often, research is also carried out on coupled phenomena, such as flow fields with temperature exchange. The study of the cooling efficiency of a radiator by air was presented by Ahmet et al. [4]. Mo et al. [5] went a step further. Their paper presents an experimental and numerical study on noise generated by axial cooling fans. The research focused on noise caused by airflow. For this purpose, the Large Eddy Simulation (LES) turbulence model and the Ffowcs Williams-Hawkings (FW-H) analogy were used. Bergamini et al. [6] conducted a Computational AeroAcoustic (CAA) study using Lighthill's acoustic analogy. A bluff body was used as a geometric model in that work. A bluff body has a trapezoidal prism shape with a flat front and a 50 degrees slanted rear. Results received from a CFD simulation were used in the Ffwocs Williams-Hawkings equation to predict external aerodynamic noise. In the validation process, a CFD simulation of a cylindrical object was performed in two-dimensional space and the results were compared with available data in the literature. As the second stage of methodology validation, an experiment was performed in a wind tunnel where a bluff body was used. Tsai et al. [7] performed an aeroacoustic analysis of a car with a rear spoiler. In that work, an unstructructed mesh was applied. A simulation of airflow around the vehicle was carried out in two stages. At the beginning, calculations were made for the steady state using the Renormalization Group (RNG) k − turbulence model. Then the results were used as the initial condition for calculations in the transient time. The aim of the work was to find the most advantageous rear spoiler configuration due to the generated noise and aerodynamic coefficients. In order to study pressure fluctuation around an Ahmed body, Duncan et al. [8] used the Lattice-Boltzmann method implemented in PowerFLOW software. The transient simulation was performed using Very Large Eddy Simulation (VLES) as turbulence model. The model mesh was automatically generated using the Cartesian volumetric grid method (cut-cell method) with a special area of refinement. The authors used spectral analysis to test pressure fluctuation in low, medium and high frequency bands.
The purpose of this work is to create an original, verified, numerical model that allows determining the acoustic field around the simplified vehicle geometry called an Ahmed body. Simulations are based on using the Large Eddy Simulation (LES) turbulence model and the Ffwocs Williams-Hawkings analogy. The frequency range in this work is 25 k Hz. The CFD calculation is performed in the commercial software ANSYS Fluent. Special attention is given to the discretization process, which is realized in the ANSYS ICEM tool. The developed model will be used for further research on acoustic and aerodynamic phenomena associated with moving road vehicles.

Mathematical Model
The purpose of this work is to study the acoustic field around an Ahmed body using available CFD tools. For this reason solving the differential equation is required. Motion of fluid is described by continuity and Navier-Stokes Equations (1)- (4). In order to reduce computational time a few assumptions are introduced to the present model. Fluid is viscous, incompressible and Newtonian, and gravitation is negligible and in this work is neglected. Use of the Reynolds Averaged Navier-Stokes (RANS) turbulent model in the first stage significantly improves model convergence. The results from k − ω Shear Stress Transport (SST) are used as an initial condition for transient simulation.
Continuity equation: ∂u x ∂x + ∂u y ∂y Navier-Stokes equation: where: u x , u y , u z -air velocity vector components p-pressure ρ p -air density µ p -dynamic viscosity Received pressure and velocity distributions are used to calculate the acting drag force F d . To compare results with other studies, the non-dimensional drag coefficient C d is calculated according to Equation (5).
where: F d -acting drag force A-frontal area of the vehicle u-relative air velocity. The initial calculations are conducted using the k − ω SST turbulence model. Transport equations for this model are presented in (6) and (7).
Transport equation for the turbulence kinetic energy k: Transport equation for the specific dissipation rate ω: where: µ t -turbulent viscosity S-modulus of the mean rate-of-strain tensor Further CFD calculations were performed using the LES turbulence model. The governing equations were obtained by filtered continuity and Navier-Stokes (8)-(11) equations. The Smagorinsky-Lilly model was used to model the eddy-viscosity µ s . Filtered continuity equation for LES: Filtered Navier-Stokes equation for LES: where: U x , U y , U z -filter velocity vector components P-filter pressure µ s -eddy-viscosity |S|-invariant of the filtered-field deformation tensor L s -subgrid-scale characteristic mixing length κ-Karman constant d-distance to the closest wall C s -Smagorinsky constant ∆-local grid scale V-volume of the computational cell.
In order to predict the acoustic field around an Ahmed body the extended form of Lighthill's analogy (known as the Ffowcs Williams-Hawkings equation) is used according to Equation (12). The searched parameters are the pressure fluctuations p , which are hidden behind the density fluctuation ρ .
Ffowcs Williams-Hawkings equation: where: Because the flow is incompressible it is possible to assume p = c 2 0 (ρ − ρ 0 ) = c 0 ρ . Furthermore, the velocity components normal to the surface are u n = 0 i v n = 0. Using the above assumptions, the Ffowcs Williams-Hawkings equations simplify to Equations (13)- (15). 1 1 As a result from the calculations, the pressure fluctuation p is received as a function of time. The pressure fluctuation is used to calculate the sound pressure level L sp according to Equation (16). 10 p p re f [dB] (16) where:

Numerical Model
In order to study the flow field around vehicles and phenomena of flow separations, in 1984 Ahmed et al. proposed a simplified geometry of a vehicle [9]. The model is characterized by a rectangular shaped body with a rounded front and a slanted top part of the rear. The most common model also has four cylindrical pillars under the geometry ( Figure 1). To this day the model is in use as a reference geometry for scientist all over the world, and it is known as the Ahmed body.
The Ahmed body has many applications in the scientific community. Many scientific works are focused on the numerical investigation of aerodynamic parameters depending on the turbulence model used. Guilmineau et al. [10] compared k − ω SST, the Explicit Algebraic Stress Model (EARSM), Detached Eddy Simulation (DES) and Improved Delay Detached Eddy Simulation (IDDES) models. The last two are called hybrid models due to joining the advantages of the RANS model with the complexity of transient models. Investigation has been conducted for Ahmed body geometry with two different slant angles, 25 and 35 degrees. The best results were received from the IDDES approach. Serre et al. [11] presented studies on a reference structure with a rear slant angle of 25 degrees. They performed four simulation processes using LES and DES models, achieving overestimates of the drag coefficient between 6-16%. Additionally, the mesh used in the simulation depends on the used turbulence model. Banga et al. [12] used a k − realizable model with a tetrahedral mesh. The studies concentrate on drag and lift coefficients depending on the rear slant angle with an inlet velocity of 40 m/s. Ten different angles were tested ranging from 0 to 40 degrees. The minimum drag coefficient of 0.235 was obtained for the slant angle of 7.5 degrees. Meile et al. [13] presented studies on aerodynamic parameters for slant angles of 25 and 35 degrees. Numerical simulations used the Reynolds Stress Model (RSM) of turbulence for a velocity range of 10 to 40 m/s. The cut-cell mesh was implemented. Results were compared with experimental data. Thomas et al. [14] compared five different turbulence models to compute airflow around an Ahmed body with a rear slant angle of 25 degrees. Turbulence models employed in these work were: k − standard, k − realizable, k − ω SST, Spalart-Allmaras (SA) and Wray-Agarwal (WA). The WA model is a one-equation turbulence model implemented in the ANSYS Fluent software as a user-defined function [14,15]. The advantage of this turbulence model is less simulation time required per iteration.
In this paper, an Ahmed body is used in CFD and CAA analyses. In the first stage of this work the numerical model is validated. Experimental data from the literature are compared with the results received from CFD simulations. It is an important step in this work because the accuracy of the calculated acoustic field depends on the quality of the CFD solution (Dykas et al. [16]). The next and main step is the calculation of the acoustic field around the reference structure and on this basis the prediction of the acoustic field around real road vehicles in the future.

Geometry Model
The geometry model of Ahmed body consists of three parts: a fore body with rounded edges, a middle section with a rectangular cross section, and a slanted rear end with an angle of 25 degrees. The body is mounted on four cylindrical pillars. The Ahmed body was prepared as a 3D computer model using the part design tool in CATIA software. The dimensions of the reference structure are shown in Figure 1.

Discretization
The influences of computational domain discretization were investigated by Sosnowski [17] and Mansour et al. [18]. Comparing three types of grid elements, namely tetrahedral, hexahedral, and polyhedral meshes, a hexahedral mesh was found to be the least numerical diffusive. Furthermore, using this type of element helps with shorter convergence time. For this reason the computational domain was discretized using a hexahedron mesh that was manually implemented using ANSYS ICEM software (Figure 2). The number of elements depending on air flow velocity ranged from 0.5 to 5.1 million. The boundary layer consisted of 15 elements in height, where the height of the first element as a nondimentional coefficient y+ for 97% of elements was below 1 (17). where: The cell Courant number calculated from Equation (18) for 99.9997% of elements was below 1. The minimum orthogonal quality was 0.323. Additionally, between the boundary layer and the domain a transitional layer was applied. The transitional layer consisted of 15 elements in height. The growth ratio between the elements in both layers was 1.15, while the aspect ratio did not exceed 60.
where: u ∞ -free stream velocity ∆t-time step ∆x-cell size.

Boundary Conditions
A finite-difference approximation requires defining a solution at the domain boundaries to receive a unique solution. This case is an example of airflow around a vehicle model in a channel. In order to avoid interaction of the airflow structures caused by the Ahmed body shape with channel planes, it is necessary to ensure the appropriate size of the computational domain. For this purpose, a 1 L × 2 L × 8 L calculation area was created around the Ahmed body geometry, where L is the vehicle's characteristic length of 1.044 m. The major turbulence and swirl structures of the airflow are observed in the wake behind the vehicle. This area is defined as five vehicle lengths in the x-coordinate. Because of a high pressure area at the front of the vehicle model, the distance between inlet plane and Ahmed body surface is defined as two characteristic lengths. Similar recommendations on the size of the computational domain are used by the European Research Community on Flow, Turbulence and Combustion (ERCOFTAC).
The computational domain with surfaces names is presented in Figure 3. The no-slip condition was applied to the surface of the Ahmed body (Surface G) and bottom plane (Surface A). The front plane (Surface B), where the fluid enters the computational area, was selected as the inlet boundary condition. Depending on the simulation, air velocity was in range from 20 to 100 m/s. The pressure gradient was zero for this surface. The outlet boundary condition was employed on the plane where the air leaves the flow field (Surface E). A reference pressure of zero was assumed. Zero gradients were applied for other parameters. A symmetry condition was applied along the middle section of the body for plane y = 0 (Surface F). This results in reducing the number of elements and increasing the speed of calculations. The boundary condition with the coordinates of each surface are described in detail in Equations (19)-(44). (Surface A) The bottom surface of the channel is defined as a stationary wall with a no-slip condition: where: k w -turbulent kinetic energy in the wall cell ω w -specific turbulence dissipation in the wall cell y w -distance from wall to cell centroid. (Surface B) At the inlet of the channel the air velocity is fixed to a constant value. Depending on the case the velocity is from 20 to 100 m/s: (a) Studies of pressure and velocity distribution, streamlines in symmetry plane y = 0, comparison of vortex system (Section 3.1), streamwise velocity profiles, drag and lift coefficient as a function of time (Section 3.2), overall sound pressure level and A-weighting sound pressure level (Section 3.4): where: where: I = 1%-turbulence intensity µ t µ p = 10-turbulent viscosity ratio.
(Surface C) The side of the channel is defined as a symmetry: for − 3.132m ≤ x ≤ 5.22m, y = 1.044m, 0 ≤ z ≤ 2.088m (Surface D) The top of the channel is defined as a symmetry: (Surface E) At the outlet of the channel the pressure is fixed: where: u avg -mean flow velocity I = 5%-backflow turbulent intensity µ t µ p = 10-backflow turbulent viscosity ratio.
(Surface F) The side of the channel with an Ahmed body contour is defined as a symmetry: (Surface G) The surface of the Ahmed body is defined as a wall with a no-slip condition: ∇k · n = 0, ∇ω · n = 0, ∇p = 0 (42) where: n-unit vector normal to surface.

Initial Condition
To solve the parental equation for unsteady flow using the LES (Large Eddy Simulation) turbulence model, an initial condition was required. In the beginning for t = 0, the velocity vector at each point is equal to the velocity vector calculated in the steady-state using the k − ω SST turbulence model (45). Pressures at each surface element were calculated every two computational time steps.
The maximum frequency related with the Nyquist frequency was 25000 Hz. These sound frequencies cover the full range of human hearing.
where: i-grid index on the x axis j-grid index on the y axis k-grid index on the z axis.

Results
In the presented numerical model, a hexahedral mesh was used in a rectangular channel. For this reason, the airflow was aligned with the mesh for most of the domain area. This is important due to the reduction of the numerical discretization error and allowed the use of the first-order accuracy scheme in the first stage of calculations. The first-order accuracy scheme yields better convergence so that the required criteria could be obtained in a lower number of iterations. It should be noted that at the same time the results yielded less accuracy, and therefore this scheme was changed later. The desired convergence criteria for this stage was residuals below 10 −4 and it was reached after 411 iterations. In the second stage of calculations a second-order accuracy scheme was used. The same convergence criteria were assumed, where 821 iterations were required (two times more). The last stage is connected with a transient simulation using LES as a turbulence model. In addition to suitable low residues below 10 −4 , a condition of a steady value of drag and lift coefficient with fluctuation around the mean was assumed. Both criteria were met after approximately 7500 iterations (Figure 4). For further calculation, we assumed 10,000 iterations. The mean value for the drag and lift coefficient was calculated from a time band of 0.1 to 0.35 s. Details of the simulation settings, used methods, and convergence criteria are described in Table 1. The first part of the results is related with the validation process. The purpose of this stage is to choose the appropriate simulation parameter. Important factors in CFD calculation include the assumed model simplifications, model discretization (mesh type and density, boundary layer, cell quality), turbulence model, boundary condition, time step, and initial conditions (the last two refer only to transient simulations).

Velocity and Pressure Field in Symmetry Plane
As results of the steady-state simulation using the k − ω SST turbulence model, the pressure field ( Figure 4) and velocity field (Figures 5 and 6) were received. The results from the simulation are compared with studies by other authors. The similarity criteria are characteristic structures generated behind an Ahmed body (recirculation and edge vortices) and velocity profiles in the symmetry plane, as well as drag coefficients. Figure 7 schematically shows the vortex structures: the separation bubble above the rear part of the vehicle, two recirculation bubbles behind the geometry, and the c-pillar vortex generated on the rear side edge of the Ahmed body. Similar structures were observed from the simulation in Figures 6  and 8. Additionally in Figure 8 the results are presented in isometric view. Area A represent the separation bubble, areas B1, B2 are recirculation bubbles, and area C is the c-pillar vortex.

Steamwise Velocity Profile in Symmetry Plane
In order to compare results from the numerical simulation with the experiment, the streamwise velocity profile were used. The velocity profile shows velocity vectors for specific coordinates marked as a straight rake of black dots in Figure 9. The distance between a rake of black dots and the results is the value of the velocity vector in the x direction. Streamwise velocity profiles were compared around the reference geometry for two different turbulence models: the k − ω SST (blue rake of dots) and LES (red rake of dots). The air velocity was assumed as 40 m/s. The reference point for the simulation tests was the experimental data (black lines) available from Lienhart et al. [19].
The results received using the k − ω SST turbulence model allowed to map the velocity field very precisely, especially for the front and middle parts of the body. In the wake a deviation is observable from the experimental data.
Further calculation using the LES turbulence model allowed us to achieve a much improved map of the wake. The absolute errors in the x direction for the k − ω SST and LES turbulence models are shown in the Figures 10 and 11. To facilitate the comparison of data, a velocity colormap is used. In order to compare results from the experiment with results from the LES unsteady turbulence model, the arithmetic average of the velocity was taken. The results from the simulation for the LES turbulence model were stored every 1305 time steps, where each time step was 10 −5 s.
After 10 4 iterations, the drag and lift coefficients were stable ( Figure 12) and changing around the mean value. The streamwise velocities in that case were calculated as the average of the velocity results from the time band of 0.1 to 0.35 s.

Drag Coefficient
Simulation studies were carried out for nine different speeds in the range from 20 to 100 m/s for the k − ω SST turbulence model and two different speeds, 40 and 100 m/s, for the LES turbulence model. The drag coefficients were compared with the results of other researchers [13,20] and presented as a function of the Reynolds number ( Figure 13). The obtained relative error did not exceed 8.35% comparing to the Meile fit. For the aeroacoustic calculation the free stream velocity was 40 m/s. In that case the drag coefficient relative errors for the k − ω SST and LES turbulence models were 0.16 and 0.06%.

Sound Pressure Level (SPL)
Assuming windless weather conditions, the speed of the air relative to the road is zero. For this reason in everyday road situations the vehicle is a source of aerodynamic noise. However, in the case of the studied Ahmed body geometry it is quite different. We will use a numerical model with a reference structure to compare the simulation data with a future experimental study. Therefore two surfaces were chosen as the sources of sound in this model. The first is connected with the surface of the Ahmed body, which is irregular. The second is the bottom plane under the body. The generated pressure fluctuation was measured at the mid height of the Ahmed body core for the z-coordinate of 194 mm. Around the geometry 39 receivers were defined, which are shown in Figure 14. The receiver number is located at the center of the data point. The distances between receivers are 1, 2 and 4 m, except for receivers that are located close to the Ahmed body surface. The overall sound pressure level (OASPL) was calculated according to Equation (16) to one decimal place, where p is substituted by the root mean square pressure fluctuation from the time range of 0.1 to 0.35 s. The results are also shown in Figure 14 but in a northwestern position from the data point. To facilitate analysis a colormap is used. Except for receivers located just next to the surface of the Ahmed body, the highest OASPL is observed in the wake, behind the geometry.
The A-weighting sound pressure level for three receivers, 12, 13, and 19, is presented in Figure 15. Points 12 and 13 are located behind the Ahmed body in an area of a high turbulence structure. Point 19 was chosen on the side in the distance of 1 m from the vehicle surface. Similar localization of receivers will be used in a future study, which will be concentrated on the emission of aerodynamic noise generated by a moving column of vehicles to the environment.
The highest values of A-weighting sound pressure levels were observed for low frequencies in the range from 180 to 1400 Hz in all three positions.

Discussion and Conclusions
In this work, numerical calculation was used in order to study an Ahmed body model. The purpose of the work was to build a numerical model using two turbulence models, k − ω SST and LES, as well as its verification based on experimental measurement data. The model quality was verified in two stages. The first stage was to prepare and ensure the good quality of the constructed computational mesh. The boundary layer was defined with a corresponding parameter y+, which was around 1. The growth ratio between the elements was set at 1.15. The minimum orthogonal quality parameter was 0.323, while the aspect ratio did not exceed 60. The second stage of the model verification consisted in comparing the simulation calculation with the experimental measurement data, as well as with the numerical results obtained by other authors. The main criterion adopted was the comparison of the streamwise velocity profiles in the symmetry plane and the drag coefficient C d . In addition, the steady-state calculation confirmed the characteristic structures formed behind the vehicle model. The numerical model developed in the work, using the LES turbulence model, is of good quality. The two-stage research process related to the development of the model and the verification by experimental study allows to clearly state that the developed numerical model is correct. Further CAA analysis was possible using the FW-H analogy.
FFT analysis transformed pressure versus time characteristic to frequency domain. In order to reduce acoustic noise generated by the Ahmed body geometry it is important to eliminate high values of sound pressure levels from the frequency range of 180 to 1400 Hz.
The introduced methodology of using a colormap to compare of results of velocity profiles simplified the analysis of troubling areas of the calculation grid. For further analysis using the k − ω SST turbulence model the mesh in the wake area should be improved. In case of using the LES turbulence model the absolute velocity error was partly the result of the method used to average the results every 1305 time steps.
When analyzing the drag coefficient in respect to other studies, the biggest relative errors were observed for higher velocities, especially for the LES turbulence model. This kind of problem could be related with insufficient quality and density of the mesh. One possible solution would be using adaptive mesh refinement for automatically increasing the density in high turbulence areas [21].
Computing the steady-state flow using the k − ω SST turbulence model and the setting results from this simulation as initial conditions for transient simulation was a useful method to reduce computation time. However, in order to evaluate computational resource reduction further calculation is needed.
The study on an aeroacoustic numerical model of an Ahmed body was performed as part of a calculation grant in CYFRONET on the computing cluster PROMETHEUS. The total used walltime exceeded 60,000 h, where the number of cores reached up to 24 for a single simulation. ANSYS Fluent software was used as a solver. Further calculations and visualizations of the results were made in MATLAB software. Created an aeroacoustic model of the Ahmed body will be used in a future study on aerodynamic and acoustic phenomena connected with a moving column of vehicles.  Acknowledgments: This research was conducted with support of the PL-Grid (Polish NGI) Consortium, coordinated by the ACC Cyfronet AGH, grant id: aia.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: