Abstract
The three-dimensional turbulent mean flow and acoustic field of a supersonic jet impinging on a solid plate is studied computationally using the general purpose CFD code Ansys Fluent. A pressure-based coupled solver formulation with the second order weighted central-upwind spatial discretization is applied to compute transient solutions. Cold and hot jet thermal conditions are considered. Mean flow characteristics are investigated by a steady-state modeling approach. Acoustic radiation of impingement tones is simulated using a transient time-domain formulation. The effects of turbulence in steady-state are modeled by the SST k-ω turbulence model. The Wall-Modeled Large-Eddy Simulation (WMLES) model is applied to compute transient solutions. The near-wall mesh on the impingement plate is fine enough to resolve the viscosity-affected near-wall region all the way to the laminar sublayer. Nozzle-to-plate distance is parameterized in the model for automatic re-generation of the mesh and results. Steady-state predictions of hover lift loss and mean jet velocity distributions are compared with experimental data, and favorable agreement is reported. The transient solution reproduces the mechanism of impingement tone generation by the interaction of large scale vortical structures with the impingement plate. The acoustic near-field is directly resolved by Computational Aeroacoustics (CAA) to accurately propagate impingement tone waves to near-field microphone locations. Calculated impingement tone frequencies and sound pressure levels agree with experimental values.
Keywords
Introduction
The accurate numerical prediction of impinging jet flows can be a valuable tool in the analysis of a short takeoff and vertical landing (STOVL) aircraft, rocket take off and aircraft exhaust. Broad-spectrum, easy-to-implement numerical models are being sought to augment engineering analysis of applications involving impinging supersonic jets. The physics of supersonic jet impingement on a surface is complicated by very strong impingement acoustic tones generated by the interaction of large turbulent structures of the jet with the impingement surface, which may affect performance of vehicle components exposed to the impingement tones and potentially even induce undesired resonant vibrations of the structure thus leading to vehicle control instabilities. It is essential for the analysis to accurately reproduce the mechanism of impingement tone generation and propagation into the near- and far-field, and its impact on vehicle structures.
A series of experimental tests1–3 provide verification and validation base for the present numerical investigation. These experiments introduced a lift plate flushed with the nozzle to represent an underbody of a STOVL vehicle thus allowing for a more realistic reproduction of flow and acoustic conditions between the impingement surface and vehicle structure induced by a supersonic jet in hover. By considering both cold and hot jet conditions, the experiments also gave insight into the thermal effects on the mean flow and impinging acoustic fields, thus providing a unique comprehensive set of test data for benchmarking simulation approaches.
This study presents an extension of the previous works4–7 on numerical simulation of mean flow and acoustic radiation of an impinging supersonic jet. The focus of investigations4,5 was on the validation of numerical approaches and establishing best practices to analyze supersonic jet impingement noise in the framework of the general-purpose CFD solver Ansys Fluent. 8 The effects of inclination angle on mean flow and noise radiation from an impinging jet were presented in publication. 6 A simulation approach 7 to the problem of mitigating impinging jet noise by aqueous injections presented a novel multiphase modeling methodology capable of including water jet injections into the supersonic flow.
In this work, an accurate and robust simulation approach for predicting mean flow effects and acoustics field of a supersonic impinging jet, which incorporates recent numerical and modeling advances in Ansys Fluent, is presented. Mean flow properties and acoustic characteristics of cold and hot supersonic impingement jets are studied using the pressure-based coupled algorithm. The problem formulation, geometry and flow conditions are described in the problem description section, an overview of the solver algorithm is given in the numerical model section, and finally the numerical predictions and their comparisons with test data are presented in the numerical results and comparison with test data section.
Problem description
The model geometry (Figure 1) and flow conditions correspond to those used in the experimental tests.1–3 A supersonic jet is discharged from a convergent-divergent axisymmetric nozzle, with the throat and exit diameters, D and D
e
, of 25.4 mm and 27.5 mm, respectively. A circular plate with a diameter of 10D and a thickness of 0.787D (20 mm), referred to as the lift plate, is flush mounted with the nozzle exit. The impingement plate is located at a distance h from the nozzle exit. H is varied from h/D = 4.0 to h/D = 8.0 in steady-state simulations focusing on predicting hover lift losses. H is fixed at h/D = 4.0 in all the unsteady simulations, where impingement acoustic tone predictions are of interest. The effects of varying h/D on impinging jet tones will be a subject of a future investigation. Model of the nozzle, lift plate and impingement surface.
The nozzle has a design Mach number of 1.5 and is simulated at a nozzle pressure ratio NPR = 3.7 corresponding to an ideally expanded jet. Two jet temperature ratios, TR = 1.0 (cold jet) and TR = 1.4 (hot jet), are considered in the study. The Reynolds numbers based on D and fully expanded jet velocity V j are 1.4e+06 and 9.3e+05 for cold and hot jet cases, respectively.
Numerical model
Most numerical approaches to high-speed jet flows, e. g.9–12 employ density-based coupled formulations where the governing equations of continuity, momentum, energy and (where appropriate) species transport are solved simultaneously as a set, or vector, of equations. In this approach, density is used as a primary variable found from the continuity equation, and then pressure is obtained using an equation of state. Density-based techniques are found to be efficient when used for high subsonic, transonic, or supersonic flows, however they require modifications, such as preconditioning,13–15 in low Mach number flow regions, e. g. outside the jet, to overcome the problem of the system matrix becoming singular in the incompressible limit.
As an alternative to the density-based approach, several coupled pressure-based methods have been proposed16–19 to extend applicability of pressure-based segregated techniques to problems where the inter-equation coupling is strong. Unlike a segregated algorithm, in which the momentum equations and pressure correction equations are solved one after another in a decoupled manner, a pressure-based coupled algorithm solves a coupled system of equations comprising the momentum equations and the pressure equation. Since the momentum and pressure equations are solved in a closely coupled manner, the rate of solution convergence significantly improves when compared to the segregated solver. The coupling also makes pressure-based coupled algorithms applicable to supersonic and even hypersonic problems that cannot be tackled by a segregated approach. For example, the study 20 confirmed the ability of the pressure-based coupled solver 8 to obtain a solution to a Mach 3.5 flow over a pod exhausting a sonic counterflowing jet which favorably compared with the density-based 8 solution and with the experimental data.
The pressure-based coupled double-precision solver implemented in Ansys Fluent 8 is employed in this study.
Pressure-based coupled solver
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 momentum (steady-state) and energy equations are interpolated from the cell centers using a QUICK-type scheme8,21 based on a weighted average of second-order upwind and second-order central differencing of the variable. The implementation in Ansys Fluent 8 uses a variable, solution-dependent value of the weight factor, chosen to avoid introducing new solution extrema. Face values of pressure are reconstructed using a second-order scheme in a manner similar to a multidimensional linear reconstruction approach. 22 In this approach, higher-order accuracy is achieved at cell faces through a Taylor series expansion of the cell-centered solution about the cell centroid.
Bounded central differencing (BCD) scheme is used to interpolate face values in the momentum equations in transient runs. The implementation 8 is based on the normalized variable diagram (NVD) approach 23 together with the convection boundedness 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. The boundedness strength of BCD can be varied by changing a control parameter from 0 (pure central differencing) to 1 (the strict BCD) to balance between minimizing numerical dissipation and maintaining solution stability. The results in the present work were calculated using the control parameter set to 0.25 which has been found to be optimal. This differentiates current results from those reported in earlier works4–7 which were calculated using the strict BCD as this was the only BCD formulation available in the code at that time.
Mass flux vector is evaluated by a momentum-coefficient-weighted high-order velocity interpolation with a Rhie-Chow correction for the pressure gradient difference. 8
Gradients needed for constructing values of a scalar at the cell faces and for computing secondary diffusion terms and velocity derivatives are calculated using the least squares cell-based gradient evaluation 8 which preserves a second-order spatial accuracy.
An implicit discretization of the pressure gradient terms in the momentum equations, and an implicit discretization of the face mass flux, including the Rhie-Chow pressure dissipation terms, provide fully implicit coupling between the momentum and continuity equations. This discretization yields a system of algebraic equations whose matrix depends on the discretization coefficients of the momentum equations, 8 and it is then solved using the coupled algebraic multigrid (AMG) scheme.8,24 An Incomplete Lower Upper (ILU) smoother is applied to smooth the residuals between levels of the AMG. The ILU smoother is more expensive than standard Gauss-Seidel, but has better smoothing properties, especially for block-coupled systems solved by the coupled AMG, which permits more aggressive coarsening of AMG levels.
Gradient limiters
The solver formulations take advantage of gradient (or slope) limiters used on higher-order schemes to prevent spurious oscillations, which would otherwise appear in the solution flow field near shocks, discontinuities, or near rapid local changes in the flow field. The gradient limiter attempts to invoke and enforce the monotonicity principle by prohibiting the linearly reconstructed field variable on the cell faces to exceed the maximum or minimum values of the neighboring cells. A non-differentiable limiter 22 based on the Minmod function (Minimum Modulus) is utilized in this study to limit and clip the reconstructed solution overshoots and undershoots. Cell to face limiting direction is chosen, where the limited value of the reconstruction gradient is determined at cell face centers.
Time-marching scheme
A second order fully implicit time scheme is used to march the solution in time. 8 Implicit equations are solved iteratively at each time level before moving to the next time step. The advantage of the fully implicit scheme is that it is unconditionally stable with respect to a time step size. The time step size is chosen to maintain flow and acoustic CFL numbers on the order of unity to balance between minimizing temporal diffusion and the overall number of time steps required for the calculation to progress over a desired time interval. The time step value used in all transient calculations is Δt = 1e−06 sec, or in a non-dimensional form Δt = 0.004 D/V j , based on fully expanded cold jet velocity V j . Five sub-iterations are found to be sufficient for converging the implicit equations at each time level before moving to the next time step.
Physical models and boundary conditions
Air is modeled as a single-species calorically perfect gas. Air viscosity is defined as a function of temperature by Sutherland's viscosity law.
The pressure inlet boundary condition at the nozzle inlet specifies total pressure and total temperature which correspond to the test conditions.1–3 A uniform (constant) profiles of total pressure and total temperature are prescribed at the inlet which corresponds to a situation of no upstream boundary profile development takes place, e. g. when a nozzle inlet opens directly to a larger pressurized plenum chamber. No synthetic turbulence was generated at the nozzle inlet in transient runs. The inner wall of the nozzle is a no-slip surface, and the turbulent boundary layer develops along the wall starting at zero thickness at the nozzle inlet and reaches it naturally developed thickness at the nozzle exit. The effects of the boundary layer thickness at the nozzle on jet impinging tones are outside of the scope of this investigation.
The far-field boundary (Figure 2) in steady-state cases is treated as a pressure boundary which uses specified static pressure and extrapolates all other flow variables from the interior of the domain if the flow is locally subsonic. In supersonic regions, all flow variables including static pressure are extrapolated from the interior. An acoustically non-reflective pressure boundary condition is applied at the far-field boundary in transient runs. It is based on characteristic wave relations derived from Euler equations.25,26 To obtain the primitive flow quantities (p, u, v, w, T), reformulated Euler equations are solved on non-reflecting boundaries, along with the interior governing flow equations, using similar time stepping algorithms. Computational domain for ¼ model.
Computational meshes
A solid model of the computational domain and a three-dimensional conformal hexahedral computational mesh (Figure 3) were generated using Ansys Pre-Processing.
8
Computational mesh for ¼ model.
Computational domain has a cylindrical shape; its far-field non-reflective boundaries are positioned 17.3D from the axis and (h + 2.4D) away from the impingement plane (Figure 2). An exact thermal condition on the impingement plate was not stated in the experimental work, 2 and the impingement plane and inner nozzle walls are modeled as adiabatic no-slip walls. The outer nozzle and lift plate surfaces are adiabatic slip walls, since there is no flow along these walls and viscous effects are negligible when generated acoustic waves reflect off the lift plate.
The mesh size distribution was determined by this work’s objective which focuses on both the impingement tone generation mechanism and near field propagation of the tones, and resolution of the impingement jet flow over the impingement plate. The mesh size at the nozzle exit is equal to 8.0e-03D in the radial direction, and 3e-02D in the axial direction towards the impingement plate. Mesh size is smoothly increased away from the jet mixing region into the far field to 0.06D and maintained at this resolution all the way to the far-field boundary. The mesh was built for ¼ model, and ¼ mesh with the symmetry planes is used in the steady-state runs. The mesh was copied and reflected about two symmetry planes to construct a mesh for the full 360 deg. model used in transient simulations. The total size of the full 360 deg. mesh is 26.7 million cells.
Turbulence models
RANS-based shear-stress transport (SST) k-ω 27 model was applied in the steady-state runs. The SST model effectively blends the robust and accurate formulation of the classic Wilcox’s k-ω 28 model in the near-wall region with the freestream independence of the k-ε model in the far field. In addition, the SST model accounts for the transport of the turbulence shear stress in the definition of the turbulent viscosity, thus making it an accurate and reliable choice for a wider class of flows, including those with adverse pressure gradient and shock waves.
The Algebraic Wall-Modeled LES (WMLES) model with the S-Omega formulation 8 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 modified LES formulation,29–31 which avoids the very high-resolution requirements of LES in the inner wall layer. The WMLES formulation combines a mixing length model with a modified Smagorinsky model and with a wall-damping function. 32 The S-Omega enhancement to the WMLES overcomes the artifact of having non-zero eddy viscosity for flows with constant shear when a modified Smagorinsky model is used in the LES zone. This is achieved by computing the LES portion of the model using the absolute difference |S - Ω | instead of S, where S is the strain rate and Ω is the vorticity magnitude. This enhancement enables the WMLES to compute transitional effects and accurately produce eddy viscosities in separating shear layers like those from the nozzle lip in the present work.
Numerical results and comparison with test data
The numerical runs were carried out in parallel on six 20-core nodes of a Linux cluster (total of 120 CPUs). To provide an initial flow field to steady-state simulations, the numerical solution is initialized with ambient conditions and the nozzle inlet value of pressure is patched inside the nozzle. Then, the full multigrid (FMG) initialization
8
is utilized to obtain the initial solution. FMG initialization is based on the full-approximation storage (FAS) multigrid technology.8,33 FMG procedure constructs several grid levels to combine groups of cells on the finer grid to form coarse grid cells. FAS multigrid cycle is applied on each level until a given order of residual reduction is obtained, then the solution is interpolated to the next finer grid level, and the FAS cycle is repeated from the current level all the way down to the coarsest level. This process is continued until the finest grid level is reached. FMG initialization is relatively inexpensive since most of computational work is done on coarse levels, which allows to obtain a good initial solution which already recovers some flow physics, as illustrated in Figure 4, thus reducing the number of iterations required to converge the solution to its steady sate. Mach number field of the cold jet for h/d = 4.0, TR = 1.0 (a) initial field after the FMG initialization, and (b) converged steady-state solution.
Converged steady-state RANS solution is used as an initial flow field for the transient simulation. The solution is then time-marched until it locks into a self-sustained impingement tone feedback loop.
Steady-state results – hover lift loss
When a STOVL aircraft hovers close to the ground, the entrainment flow induced by the jet can create a significant download resulting in a reduction in lift force. This effect is commonly referred to as hover lift loss.
34
Figure 5 qualitative illustrates the effect of hover lift loss by plotting calculated pressure field under and above the lift plate at TR = 1.4. Pressure under the plate is lower than pressure above suggesting an overall downforce on the plate. Velocity vector plots in Figure 6 show entrainment flow pattern and flow recirculation under the plate. As the ambient air is entrained under the plate by the jet flow, it separates off the corner of the lift plate and forms a recirculation zone characterized by lower pressure which contributed to the overall downlift force. Normalized static pressure, p/pa - 1.0, under and above the lift plate. H/D = 4.0, TR = 1.4. Velocity vectors, colored by normalized velocity magnitude |U|/V
j
showing entrainment flow pattern and recirculation under the plate. Vectors are scaled to the same length for better visualization.

Hover lift loss is quantified in Figure 7 by comparing calculated lift-loss variation as a function of h/D, which is varied from 4 to 8 in the numerical study, with experimental data
2
for two temperature ratios TR = 1.0 and 1.4. The lift force in Figure 7 is normalized by the isentropic thrust of the jet. There is a higher lift loss for smaller values of h/D. Lift loss diminishes and becomes negligible as the nozzle and plate are moved away from the ground. There is a slightly higher lift loss induced by the hot jet, as shown in Figure 8. Numerically predicted values of lift loss compare well with the experimental data. Comparison of calculated hover lift loss as a function of h/D with experimental data. Comparison of calculated hover lift loss generated by cold and hot jets as a function of h/D.

Steady-state results – mean jet velocity
Mean velocity profiles of the jet flow calculated with the SST model are compared with experimental data in Figure 9 for h/D = 5.0 at TR = 1.0 and 1.4. The profiles in Figure 9 are plotted at two downstream locations from the nozzle exit, z/D = 0.5 and 4.0; they are normalized using the fully expanded jet velocity V
j
. There is a favorable agreement of the numerical profiles with the test distributions.
2
Spreading of the jet downstream of the nozzle is also well-captured, suggesting the SST model accurately predicts averaged turbulent field of the jet. Comparison of normalized mean velocity distribution for cold and hot jets with experimental data, h/D = 5.0.
Transient results – impingement tone radiation
Resolved turbulent flow structures in the jet shear layer and wall jet flow along the impingement surfaces can be seen in Figure 10 showing an iso-surface of the Q-criterion (Q = ½(S
2
− Ω
2
), S being the strain rate and Ω is the vorticity), colored by normalized transverse velocity magnitude. Turbulent structures by Q-criterion.
Instantaneous pressure contours on a plane cut through the jet axis in Figure 11 illustrate generation and near-field propagation of impingement acoustic waves at TR = 1.4. The pressure field in Figure 11 is plotted on a log scale to better visualize larger pressure perturbations of the mean flow, and much smaller scale acoustic waves. The feedback mechanism by which impingement tones are generated35,36 can be observed in the numerically predicted flow field (Figure 11). The jet shear layer at the nozzle exit is thin and susceptible to external excitations. Once excited, its perturbations propagate downstream as instability waves which quickly grow in amplitude by extracting energy from the mean flow. Impingement of these large-scale structures on the wall generates coherent pressure fluctuations, which then propagate as strong acoustic waves. Upon reaching the nozzle exit, these sound waves excite the shear layer again, which generates another set of instability waves and closes the feedback loop. Instantaneous pressure field plotted on a log scale.
The acoustic field between the lift plate and the impingement plane is complex due to the repeating reflections of generated impingement waves by the plate and the ground.
Narrowband spectra of calculated sound pressure levels (SPL) measured on the lift plate at x/D = 2.0 for cold and hot jet cases are compared with the experimental data in Figure 12. The frequency and SPL of the fundamental impingement tone at 6 kHz is well-predicted in the iso-thermal case. A higher frequency harmonic at around 12 kHz is also picked up in the simulation. The hot jet case shows a slight deviation of the dominant frequency from its experimental value of 6.9 kHz and somewhat lower SPL level. A similar behavior of numerical results in a simulation of the same hot jet was observed by other researchers
37
who stipulated that the discrepancy is likely due to the mismatch of nozzle exit flow and thermal profiled between the numerical setup and the experiment
2
where an inline heater was used to heat the airflow to a desired stagnation temperature. A later work
38
showed that indeed nozzle exit flow profiles have a notable effect on impinging jet tone frequencies. Comparison of calculated and experimentally measured noise spectrum at x/D = 2.0 location on the lift plate: (a) cold jet at TR = 1.0, and (b) hot jet at TR = 1.4.
The physics of the temperature effect on impinging tone frequencies has been covered in detail by others 3 and will not be restated here. The argument is based on classic Powell’s formula. 39 Increase in tone amplitude of the hot jet as compared to the cold jet tones can be explained by the original Lighthill’s observation 40 which used the dynamic similarity approach to derive an empirical correlation for the overall acoustic power. Lighthill pointed out that temperature field inhomogeneities amplify the sound produced by turbulent structures, and at the same time the effects of velocity and temperature cannot be separated. Based on this argument, it is expected that a supersonic hot jet to be louder as it exhibits greater temperature variations and exhausts at a velocity which is greater than that of the cold jet at the same fully expanded Mach number. Test data reported in experimental works1–3 confirms this observation.
Figure 13 shows a reasonable comparison of numerically predicted noise spectra at the near-field microphone location at x/D = 15 with experimental data for the hot jet case, and deviations can be explained by the aforementioned argument of the nozzle exit profile mismatch with the experiment. Narrowband spectrum at the microphone location for the cold jet at the specific h/D = 4 was not clearly provided in referenced works, and thus a comparison is not shown here. Comparison of calculated and experimentally measured noise spectrum at x/D = 15 microphone location for the hot jet at TR = 1.4.
Power spectral density (PSD) of temperature at the stagnation point x/d = 0 on the ground plane shows favorable comparison with test data
2
for the dominant peak and its higher harmonics for the hot jet case, as shown in Figure 14. The frequencies of temperature peaks are consistent with those of the acoustic tones measured at near-field microphone locations. Unlike the SPL frequency prediction, the dominant thermal peak frequency is in excellent agreement with the experiment, suggesting that the nozzle exit profile effects which might have affected the near field noise are far less pronounced in the jet column along its centerline than in jet shear layers. Comparison of calculated and experimentally measured temperature spectrum measured at the stagnation point x/D = 0 on the ground plane for the hot jet at TR = 1.4.
In all the presented results, higher frequency tones above about 15 kHz are not predicted in the simulation due to the lack of required spatial resolution.
Summary
The results reported in this work show that the pressure-based coupled solver (PBCS) proves to be a robust and effective method to resolve the physics and capture essential features of impinging jet flow and acoustic field, including impingement tone radiation. PBCS is less memory and CPU intensive than a traditional density-based approach, which makes it an attractive alternative to the density-based algorithm for obtaining transient solutions to supersonic jet noise problems.
The Wall Modeled LES model WMLES reduces the stringent and Reynolds number-dependent grid resolution requirements of classical wall-resolved LES in the near-wall region, thus making it an attractive choice for problems involving both free shear and wall bounded flows like the one of the impinging jet.
Steady-state predictions of hover lift loss and mean jet velocity distributions are in favorable agreement with experimental data. The predicted frequency and SPL of the dominant impingement tone matches the experimental data. Better predictions are observed in the cold jet case whilst hot jet tone prediction deviations are likely attributed to nozzle exit condition differences between experimental and numerical set ups.
The study of the effects of varying h/D on impinging jet tones will be a subject of a future investigation which may include a fully dynamic set up of the nozzle-underbody geometry moving away from the impingement surface driven by the thrust-gravitational force balance accounted for in the six-degrees-of-freedom framework.
Footnotes
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
The author(s) received no financial support for the research, authorship, and/or publication of this article.
