Abstract
In this paper, we explore the use of a second-order unstructured-grid, finite-volume code for direct noise prediction. We consider a Mach 1.5 jet impinging on a perpendicular flat plate. Hybrid LES-RANS simulations are used to calculate directly both the flow field and the radiated sound. The ANSYS Fluent commercial code is utilized for the calculations. The acoustic field is obtained directly from the simulations and is compared with the integral approach of Ffowcs-Williams–Hawkings. Results indicate the existence of a preferred radiation angle. The spectrum obtained is in good agreement with observations. This points out to the possibility of handling the effects of complicated geometries on noise radiation by using unstructured second-order codes.
Introduction
The objective of this work is to propose a methodology based on the commercial ANSYS Fluent code for the study of supersonic impinging jets. Impinging jets are important for short take-off and vertical landing of military aircraft, as well as for rocket’s lift-off. High-speed impingement leads to significant unsteady loads on the vehicle surface and it elevates the noise levels. In addition, the proposed methodology can be used for other applications such as the prediction of the noise associated with free-jet plumes.
Impinging jets are dominated by a well-known feedback loop.1–8 The feedback loop starts as instability waves in the initial jet shear layer, which amplifies and becomes vortical structures. The vortical structures are reflected upon impingement, create large pressure fluctuations and travel back upstream as acoustic waves, and when they reach the nozzle exit, they excite the shear layer instabilities completing the feedback loop, as demonstrated in Figure 1. The experimental investigations of this problem have shown that strong tones are generated by the resonance/feedback phenomenon described above.

In the past, both free and impinging jet noise studies relied heavily on experiments.10,11 However, with the tremendous improvements of computing power, there are several attempts to numerically compute the generated noise. Direct Numerical Simulations (DNS) resolves all the relevant length scales of turbulence without any turbulence models. Due to the wide range of time and length scales in high Reynolds number turbulent flows, the application of DNS to realistic jets remains prohibitively expensive. Large-Eddy Simulations (LES) is a cheaper alternative because it computes the dynamics of the large scales directly and models the effect of the small (subgrid) scales. However, including the nozzle and the impinging surface (if any) require a very fine grid near the solid boundaries. The Hybrid LES-RANS (HLR) consists of unsteady Reynolds average Navier-Stokes (RANS) calculations near the solid boundaries and LES elsewhere. The method is less computational intensive compared to LES. The first use of LES for supersonic jet noise prediction was done by Mankbadi et al. 12 Many researchers have also used LES for supersonic jet noise prediction.13–16 RANS is usually used in the nozzle.17,18 However, including the nozzles in LES can become very demanding and require hundreds of millions of grid points.19–23 For impinging jets, there also have been several attempts to simulate the flow field. For example, Erwin and Sinha 24 using a hybrid RANS-LES implementation, Brown and Frendi 25 used a detached-eddy simulation (DES) method, and Khalighi et al., 26 Nonomura et al., 27 and Uzun et al. 28 used LES.
Once the near field has been computed by one of the approaches described above, the farfield noise can be calculated using integral acoustic methods (Lyrintzis 29 ), such as the Ffowcs-Williams–Hawkings (FWH) method, 30 Kirchhoff’s method (Lyrintzis 31 ), or Lighthill’s acoustic analogy (Lighthill 32 ). The first two methods are preferable because they involve only surface integrals. Early applications of integral methods for jet aeroacoustics are given by Mankbadi 33 (acoustic analogy), Lyrintzis and Mankbadi 34 (Kirchhoff method), Mankbadi et al. 35 (Integral Methods), and Uzun et al. 36 (FWH method).
While there has been considerable progress in computational aeroacoustics over the last two decades, LES is still restricted to simple geometries. 37 For instance, it is still difficult to study the effect of the internal geometry on the radiated sound field. In the present work, we use the ANSYS Fluent commercial code that has the ability to handle complicated geometries that is needed to represent the launch pad and the details of the nozzle geometry. It has been demonstrated that the code can predict the time-averaged flow of a supersonic impinging jet. However, we explore herein its ability to predict directly the near sound field using HLR. We start by considering an impinging perpendicular jet with focus on accurately capturing the generated pressure fluctuations. We focus this work on comparing the numerical results with the experiment of Krothapalli et al. 6 to ensure the validity of the numerical approach.
We present in next section the numerical approach in which we discuss the Hybrid LES-RANS approach, the grid, the FWH approach, and further numerical details. In the following section, we present the flow and the acoustic results. This followed by the conclusions.
Numerical approach
Hybrid LES-RANS simulations
Algebraic wall-Modeled LES Model (WMLES) with the S-Ω formulation
38
was applied in transient simulations. It is based on the concept of covering the inner portion of the boundary layer by a RANS and the outer portion by a LES model. The WMLES is an attractive alternative to the classical LES approach which lessens stringent near-wall mesh resolution constrains imposed by the traditional LES methodology. In WMLES, RANS and LES formulations are combined to have RANS covering the very near-wall layer (SST
39
used in this work), and then switching to the LES when mesh resolution is sufficient to resolve local scales. The WMLES implementation in ANSYS Fluent is based on the formulation of Shur et al.
40
The eddy viscosity,
One of the deficiencies of the Shur’s WMLES model 40 is it does not produce zero eddy viscosity for flows with constant shear when a modified Smagorinsky model is applied in the LES zone. It does not allow the prediction of transitional effects and can produce excessive eddy viscosity in separating shear layers. An enhancement to the WMLES overcoming these deficiencies is based on computing the LES portion of the model using the difference of the absolute value of S – Ω instead of S, where Ω is the vorticity magnitude. This enhancement is referred to as the WMLES S-Ω formulation, and it is used in the current work.
Numerical details
The simulations were carried out using ANSYS Fluent R14.5. The Fluent code is a full Navier-Stokes solver which implements the finite-volume approach. Details about the implementation of the solver can be found in the fluent documentation. 38 The low Mach number flow far away from the impinging regions slows down convergence of the solution. Therefore, a pressure-based coupled solver is used in this study, which avoids the need for preconditioning otherwise require by the density-based approach. The coupled option was chosen for the solver since the inter-equation coupling in this system is strong. The governing equations for the conservation of mass, momentum, and energy are discretized using a control volume-based technique. Face values required for computing the convection terms in the continuity and energy equations are interpolated from the cell centers using the second-order upwind and the QUICK type scheme, 38 respectively. Face values of pressure are reconstructed using a second-order scheme in a manner similar to a multidimensional linear reconstruction approach. 41 The implementation is based on the normalized variable diagram (NVD) approach 42 together with the convection roundedness criterion (CBC). The bounded central differencing scheme is a composite NVD scheme that consists of a pure central differencing, a blended scheme of the central differencing, and the second-order upwind scheme, and the first-order upwind scheme. It should be noted that the first-order scheme is used only when the CBC is violated. To reconstruct the face values, it is necessary to calculate the gradients of the flow variables. In addition, the velocity and temperature gradients are required to compute the viscous flux. The nodal-based Green-Gauss approach is adopted in our simulations.
A second-order fully implicit time scheme is used to march the solution in time. The scheme is unconditionally stable. The time step was selected to ensure that the important frequencies are appropriately resolved. A constant time step 2.0 × 10−6 s is used. When non-dimensionalized by the nozzle throat diameter and exit velocity, this time step corresponds to 0.034, which gives a maximum CFL around 5 near the shear layer regions. The solution initialization for the transient case is done in two steps. Firstly, a full multigrid initialization is carried out to provide an approximate steady-state solution. Next, a more accurate RANS simulation is carried out to provide a better estimate of the flow field. The HLR is carried out after this step.
At the nozzle inflow boundary, the total conditions are specified to calculate the static thermodynamic properties. This is done by using pressure inlet boundary condition with the input total pressure and total temperature. The nozzle interior, lift plate, and the impingement walls are specified as adiabatic no slip wall boundaries. Since there is no freestream, the farfield boundary is specified by a given static pressure, which is the ambient pressure. This is done by using pressure outlet boundary condition. In order to avoid the wave reflections from the farfield boundaries to contaminate the interior flow fields, the characteristics-based non-reflecting formulation43,44 is used.
Computational mesh
The computational mesh is a full three-dimensional mesh, which contains about 26.61 million cells. The entire domain contains only hexahedral type cell, and most of the grid points cluster in the jet core, shear layer, and the impingement plate. The normal direction grid space on the impingement wall is 1.5 × 10−4Dt. This corresponds to y+ is less than 0.5 in the flow-attached regions. The flows are detached in most of the impinging regions, and the y+ in there is generally less than 2.0. There is no attempt to resolve the wall boundary on the lift plate; therefore, the normal grid spacing on the lift plate is about 0.029Dt. Inside the nozzle, the boundary layer is assumed laminar, and the grid spacing normal to the wall is about 0.00787Dt.
Figure 2(a) shows the grid distributions near the nozzle lip and the edge of the lift plate. The grid points are clustered in both nozzle center and shear layer regions. Even though only the hexahedral element is used, this is not a structured grid. Therefore, there is block interface existing inside the cylindrical type geometry. The grid distribution at the nozzle exit and perpendicular to the jet centerline is shown Figure 2(b). Most of the grids are clustered near the shear layer regions and the same distributions extend to the impinging wall. About at 1.18Dt away from the centerline, the grid space in radial direction is close to a constant, which is about 0.06Dt.
(a) Mesh distributions near the nozzle lip region. (b) A quarter-mesh distribution at the nozzle exit and perpendicular to the jet centerline.
The FWH acoustics integral surface
To reduce mesh resolution, acoustics in the farfield are obtained using the integral technique. The FWH equation is an inhomogeneous wave equation derived by manipulating the continuity equation and the Navier-Stokes equations. The complete solution consists of surface integrals and volume integrals.29,30 For a solid control surface, the surface integrals represent the contributions from monopole (thickness) (
If we assume that the control surface contains all acoustic sources, the volume integrals outside this surface can be dropped, and the solution in the time domain becomes (for a moving permeable surface):
The kernels of the integrals are computed at the corresponding retarded times, τ, defined as equation (10). Given the receiver time, t, and the distance to the receiver, r, the retarded time can be computed by
More detailed information can be found in the ANSYS Inc. 38
In this simulation, a cylindrical surface concentric with the jet core with a radius of 5 throat diameters was set as the FWH surface. More details are described in section 4. During the computation, the near filed variables computed by the LES are stored at the FWH surface every certain time steps. The pressure fluctuations at the observer’s locations are then computed by FWH surface integral acoustic method.
Benchmark experiment
To demonstrate the validity of the code in simulating the impingement of supersonic jets, a benchmark experiment of Krothapalli et al. 6 is used so that simulation results can be validated. The experimental setup contains a converging diverging axisymmetric nozzle with a throat diameter of 0.0254 m and an exit diameter of 0.0275 m. The converging section of the nozzle is a third-degree polynomial with a contraction ratio of approximately 5. The diverging section of the nozzle is straight walled with a divergence angle of 3°. A schematic of the entire experimental setup is shown in Figure 1. The nozzle is mounted flush with a solid circular plate called the lift plate. The impingement surface is four throat diameters, 4Dt, away from the nozzle exit. The jet exiting the nozzle is ideally expanded with a nozzle pressure ratio of 3.7. This nozzle pressure ratio is the ratio of the stagnation pressure in the nozzle to the ambient pressure. The jet is unheated with a stagnation temperature of 20℃. The microphone is placed 0.25 m away in plane with the nozzle exit. The Reynolds number based on the nozzle exit conditions and the exit diameter is about 1.56 × 106.
Results and discussions
All the simulations are carried out on Rigel at Embry-Riddle Aeronautical University, Daytona campus. This computational cluster has 22 nodes and total 280 computational cores. The processor type is Intel Xeon X5670 with the speed of 2.93 GHz. The total memory size is 982 GB. Each of our computational case uses 120 cores at one time. For the current simulations, 1000 time steps take about 18 h wall-clock time. At the beginning, the flow field was initialized by the far-field condition. The time step is constant at 2.0 × 10−7 s. In order to reduce the number of samples and saving the computational time for acoustic analysis, we then increased the time step to 2.0 × 10−6 s after 20,000 time steps. This time step gives the maximum CFL less than 5 in the entire flow filed. The maximum CFL occurs near the shear layer regions. The flow statistics were collected over 25,000 time steps to get the time-averaged results. This number of time steps corresponds to 0.05 s of flow time, which equals to an ambient sound wave travels about 675 times the impinging distance. The time-averaged results have been compared every 5000 time steps to ensure the solution changes are less than 0.2% near the impinging regions and converges properly. The computed results were examined from time-averaged flow field, and then followed by studying the acoustic field.
The flow near field
Turbulence spectra
To examine the range of the resolved scales in our LES, we show in Figure 3 the one-dimensional spectrum of the streamwise velocity fluctuations at 3.5Dt away from the nozzle exit on the jet centerline. The temporal spectrum of the streamwise velocity fluctuations at the given location is coupled with Taylor’s hypothesis of frozen turbulence to compute the one-dimensional spectrum. We collected about 25,000 time steps (840Dt/Uj) to average the results. The grid space in the streamwise direction at this region is about 0.0291Dt. The grid cut-off wave number, which shown in dashed line, is computed by assuming 20 points per wavelength, which is on the conservative side, and is about 425 m−1. Before the grid cut-off wave number, the spectrum exhibits a decay which is quite similar to Kolmogorov’s −5/3 decay rate for the inertial range of turbulence. At the same distance away from the nozzle exit, Figure 4 shows the one-dimensional spectrum of the streamwise velocity fluctuations at the jet lip line. Similar to the Figure 3, the spectrum has a −5/3 decay rate within the grid cut-off wave number. This implies that our grid resolution is fine enough to resolve a portion of the inertial range of wave numbers at the core and the shear layer locations.
One-dimensional spectrum of the streamwise velocity fluctuations at the x = 3.5Dt location on the jet centerline. One-dimensional spectrum of the streamwise velocity fluctuations at the x = 3.5Dt location on the jet lip line.

There is a peak at wave number of around 105 m−1, and the frequency is about 5920 Hz. This frequency corresponds to a single-tone noise in the pressure spectra as shown in Figure 13. Since there is a standing shock oscillating at this probe location, this 5920 Hz single tone in the pressure spectra can probably be related to the interactions of the shock and the shear layer.
Time-averaged velocities
Figure 5 shows the time-averaged centerline streamwise velocity from the nozzle exit to the impinging wall. This velocity is non-dimensionalized by nozzle exit velocity Uj, which has a value of 427.41 m/s. This exit velocity is computed from the isentropic relation based on the given total conditions. In this figure, the horizontal axis is non-dimensionalized by the nozzle throat diameter. The experimental data are taken from Krothapalli et al.
6
It can be seen from the figure that the time-averaged velocity has some oscillations before it starts to decay at about x = 3.3Dt. The experiments have the mean centerline velocity remaining normally constant up to about x = 3.0Dt. However, this difference is within an acceptable range. Figure 6 shows the time-averaged Mach number distributions along the centerline. An examination of this figure shows that the sonic condition occurs at about x = 3.58Dt. The experiment reports the sonic condition at x = 3.4Dt, where there is a change in axial velocity gradient. It is interesting to note that the experiment shows that the location of the sonic condition to shift to upstream if the nozzle NPR is increased, and that there is a standoff shock near the x = 3.5Dt regions. This standoff shock can be seen clearly from the time-averaged streamwise velocity contour as shown in Figure 7 or the instantaneous Mach number contour as shown in Figure 8.
Time-averaged streamwise velocity along jet centerline. Time-averaged Mach number along jet centerline. Time-averaged streamwise velocity contour. Snapshots of the instantaneous flow field computed using HLR simulation. (a) Pressure; (b) Mach number; (c) temperature; (d) turbulent viscosity.



Snapshots
Unlike RANS, the current approach relies on computing the time-dependent unsteady structure through the HLR techniques. We present here snapshots of instantaneous contours in Figure 8. To fully understand the dynamics of the solution, an animation of the Mach number is important to consider. We show in Figure 9 snapshots of the time histories of Mach number. A few dynamic flow features can be seen clearly from these plots. First, the end of the potential core can be seen to be fluctuating. Next, the dynamics of the wall jet shows a large circulation and high-vorticity region. This vortical structure seen in the wall jet is the main reason a high-grid resolution was chosen for this part of the computational mesh. In addition, there is a standing shock, which has a Mach number around 2.0 oscillating in front of the impinging regions.
Instantaneous Mach number contours. Each successive figure represents the next twenty time steps.
The acoustic field
Snapshots of the small pressure fluctuations in the flow field are shown in Figure 10. These small pressure fluctuations are the difference between the local mean and the local instantaneous pressure. Small pressure fluctuations in the flow field are seen as ripples in these pressure contours. This implies that these ripples are representative of the acoustics field. A closer look into Figure 10 shows that these pressure waves propagate away from the impingement region at a specific angle. This angle is a first indication of the directivity of the noise produced. One of the major things that can be inferred from these pressure contours is the absence of reflections from the boundaries of the domain. These reflections can be a major concern for most acoustics simulations.
Instantaneous acoustic pressure contours. Each successive figure represents the next 20 time steps.
The FWH surface
In order to investigate the noise generated from impingement jet, both the direct noise calculation and the FWH approach were adopted. Figure 11 shows the schematic of the computational geometries and the FWH surfaces. There are two different FWH surfaces used in the noise calculation. The first FWH surface is a ring (or cylindrical) surface concentric with the jet centerline, and is shown in blue color. Its radius equals to the lift plate’s radius, which is 5Dt. It should be noted that this ring surface does not include the top or bottom cap. The second surface is the bottom cap of the ring surface, which is the impingement wall and is shown in green color. The strategy of the placement of the ring FWH surface and the use of the bottom cap surface has also been adopted in literature.9,24 The circumferential location of the ring FWH surface is far away from the shear layer in order to minimize the non-linear effects. However, from the contour of the instantaneous vorticity magnitude shown in Figure 12, there might be additional vortical structures outside the ring FWH surface, which may cause underpredication of the noise levels by using the FWH approach.
Schematic of ring and wall Ffowcs-Williams–Hawkings surfaces. Instantaneous vorticity magnitudes. Near-field narrowband spectra.


Unlike the experiments, the noise signal generated by CFD depends on the grid resolutions. At the location of the FWH surface, the maximum grid space in both axial and radial directions is around 0.059Dt. This is also the grid space around the microphone location. Basing on this grid space and assuming 20 points per wavelength, the maximum grid resolvable frequency is about 11.33 kHz. We collected 21,000 samples (706 Dt/Uj) to analyze the acoustic field.
The spectra
For the FWH sampling, we use the same sampling time step as our LES time step. The pressure data can be converted into spectrum by using fast Fourier transform (FFT). With this sampling time step, the highest FFT resolvable frequency is 250 kHz. For long noise signals, the computed Fourier spectrum may display spurious fluctuations of amplitudes between the neighboring modes. In order to obtain a smooth spectrum, a single period sample is divided into several segments. We collected a total of 21,000 samples and then divided this single period samples into several segments. In our case, each segment contains 1040 samples. The FFT is then applied on each segment and the resulting spectra are averaged. This procedure can significantly suppress the spurious fluctuations of the spectrum. Based on the number of samples in each segment, the lowest FFT resolvable frequency is about 480 Hz.
The frequencies can be non-dimensionalized to Strouhal number by using the jet exit velocity and the nozzle throat diameter. When studying the acoustic spectra, it should be kept in mind that the high-frequency portion of the spectra is just numerical noise. The sound pressure levels (SPLs) computed by direct and FWH approach at the microphone location are shown in Figure 13. In this figure, both frequency and the Strouhal number are presented. The experimental SPL has two major tone signals at around 4610 Hz and 5920 Hz, and several small amplitude tones between 7.3 kHz and 20 kHz. The experiments also report that the tone at frequency around 4610 Hz is the nozzle screech tone, and the 5920 Hz tone is generated by the impingement. The 5920 Hz tone does not exist in the free-jet case at the same operation conditions.
All the computed spectra capture the correct tone noise at frequencies of 5920 Hz and 11.968 kHz. The later one can be regarded as the second harmonic of the impingement tone. The experimental impinging tone has maximum amplitude of 143 dB, but the calculated amplitudes are about 132.5 dB and 130.5 dB for the direct and ring FWH method, respectively. The direct method even captures the high-frequency tone at 18 kHz. The trend of the board band noise is similar to the experiment; however, the amplitudes are slightly lower than the experiments.
The estimated grid cut-off frequency is shown by a dashed line. It should be mentioned that assuming 20 points for a wavelength is a little bit conservative, and the computed grid cut-off frequency is around 20 kHz as shown in Figure 13. The ring FWH surface predicts lower SPL, about 3 to 5 dB, than the direct method. However, if we include the wall FWH surface, these losses can be recovered and the spectra become similar to the direct method. The same case was run by Brown and Frendi 25 using the Overflow code with direct sound calculation. In their report, the tones at frequency of 4610 Hz, 3325 Hz, and 5920 Hz are designated as L1, L2, and L3 tones, respectively, and shown in the same figure. Even though the amplitudes of these three tones are similar to the experimental measurements, the entire spectra are overpredicted. It should be noted that the amount of their overprediction appears to be higher than the amount of our underprediction. Some underprediction is shown also in the LES calculations of Uzun et al. 28 (for impinging distance of five diameters), although approximately 200 million points were used.
The experiments show the overall sound pressure levels (OASPL) for the impinging distances of 3Dt and 5Dt are 156.2 dB and 155.6 dB, respectively. Our impinging distance is 4Dt, and the OASPL predicted by the direct method is 153.7 dB, which is about 2 dB lower than the experimental data. The impingement wall surface can reflect the sound waves. We expect these reflected sound waves theoretically can be captured by the ring surface. Therefore, the underprediction of the SPL by the ring FWH surface may probably due to the grid resolution inside the FWH surface. The grid may be too coarse to preserve the reflected waves or the reflected waves have been dissipated before they reach the ring FWH surface. Another possible reason is that the cylindrical surface does not include all the noise sources.
Directivity
The directivity of the sound can be calculated by studying the OASPL at mics located at different angles from the jet axis. The OASPL is calculated by integrating the SPL spectrum over all frequencies, and are shown in Figure 14 for different angles. The SPL spectra for each mic were calculated by using the ring FWH surface. Even though the FWH method under predicts the amplitude of sound, it still provides good trends which can be used in estimating OASPL directivity. These microphones are distributed in a circular arc with a radius of 0.25 m. It is for this reason that the directivities are only calculated for angles greater than 70°, since anything below that lies past the impingement plate and is outside the computational domain. The angle for the directivity is calculated such that the axial direction corresponds to 0° and the radial direction corresponds to 90°. It is noticed in Figure 14 that the amplitude increases from 70 to 100° and then starts to drop off. This shows that the noise is loudest at about 100° to the jet direction.
Directivity of the sound produced estimated by ring FWH surface.
Conclusions
A supersonic Mach 1.5 fully expanded jet impinging on a flat plate is studied by a hybrid RANS/LES approach with second-order numerical scheme. The fluctuating flow field associated with sound generation and propagation was analyzed. Pressure contours were studied to view the propagation of acoustic waves throughout the domain. No artificial reflections from the boundaries were observed in the simulation.
The turbulence spectra decay as the Kolmogorov’s −5/3 decay rate in the inertial range of wave numbers. The time-averaged results are presented, and the standoff shock near the impinging wall is captured at around 3.3Dt away from the nozzle exit. For the near field acoustics, both direct and FWH methods are used to calculate the noise at the experimental microphone location. Both methods capture the tone noise at around 5920 Hz created by the impingement, and other smaller tones at higher frequencies within our grid resolvable range. The broadband noise also has a similar trend to the experiments. The spectra computed by the ring FWH surface is 2 to 4 dB lower than that predicted directly. The amplitudes of tone and broadband noise are slightly lower than the experiments. This could probably be due to the number of grid points used in our case.
The reasonable agreement with observations is encouraging enough to further explore the use of computer codes designed for handling complicated geometries for directly predicting the associated acoustic field.
Footnotes
Acknowledgments
The support of the Florida Center for Advanced Aero-propulsion (FCAAP) is greatly appreciated.
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
This work was funded by a contract from United Launch Alliances.
