1. Introduction
Despite decades of research and progress, accurate and efficient computational fluid dynamics (CFD) simulation of wall-bounded compressible turbulent flows still represents a challenging task, e.g., [
1]. Within this framework, impinging supersonic jet flows continue to receive great attention, due to their technological importance as well as high complexity, in both theoretical and applied engineering research. In particular, the impingement of a turbulent supersonic jet on an inclined flat plate appears in various aerospace engineering applications, including rocket launchers, vertical/short-takeoff-and-landing aircraft operation, and jet-engine thrust control [
2,
3]. Moreover, this particular jet flow configuration has been also employed in gas-liquid separators, heat exchanges, and mixing devices. In all these cases, the jet flow characteristics play a crucial role in enhancing mixing, as well as heat and mass transfer between the fluid phases, e.g., [
4].
Generally, using applied CFD techniques makes it less expensive to examine a wide range of flow parameters and configurations, compared to experimental works. In this particular case, notwithstanding the geometric simplicity, the interaction between the jet flow and the inclined solid surface results in highly intricate fluid dynamic phenomena, such as flow separation, reattachment, vortex dynamics, and complex interactions between shocks and boundary layers. The significant computational complexity has directed research efforts primarily towards the examination of perpendicular jet impingement, e.g., [
5]. However, recent advancements in computational power accessible to industrial research, coupled with the evolution of sophisticated numerical methods for applied CFD, have led to an increased focus on supersonic jet flow dynamics, including more complex configurations, e.g., [
6]. Specifically, engineering research has taken into account the variations in nozzle-to-wall distance, nozzle exit diameter, and impingement angle [
7]. In fact, these geometric parameters, in conjunction with the nozzle pressure ratio, play a crucial role in determining the characteristics of impinging jet flow. For instance, an increase in the angle of inclination of the plate results in heightened surface pressure oscillations, due to more pronounced flow separation and reattachment phenomena, e.g., [
8].
The present study specifically examines the scenario of near-wall impingement, characterized by a ratio of
L/
, where
L denotes the distance from the nozzle to the plate (measured along the centerline of the nozzle) and
refers to the diameter of the nozzle exit section. In this context, the pressure exerted on the plate is significantly augmented, owing to the intricate flow structures that emerge from the intensified interactions between shock waves and expansion waves. Indeed, experimental research indicates that the pressure load resulting from the supersonic jet striking an inclined solid surface can exceed that observed at normal incidence, primarily due to the presence of multiple shock waves between the nozzle exit and the plate [
9]. Furthermore, various studies have emphasized that the compressible flow structures formed at the zone of impingement are also affected by the fluid turbulence, e.g., [
10].
In the initial experimental studies focused on near-wall impingement, pressure taps were employed to assess the pressure distribution on the plate. Typically, these measurements yield pressure fields characterized by low spatial resolution, which restricts further analysis, particularly in small-scale experiments. Additionally, in certain instances, the use of physical probes for measurement is not feasible, as they may significantly modify the flow field. In response to these limitations, the pressure-sensitive paint (PSP) technique was introduced to ascertain the pressure distribution with relatively high spatial resolution, thus serving as an effective tool for the examination of jet impingement [
11]. For example, Nakai and collaborators [
12] explored the behavior of supersonic jets impacting inclined flat plates at various nozzle-to-plate distances (
), impingement angles (
), and nozzle pressure ratios (
). The latter parameter is defined as
, where
and
represent nozzle exit and ambient pressures, respectively. By integrating the high-resolution pressure maps obtained from PSP technology with Schlieren imaging, a methodology was developed to anticipate the overall impinging jet configuration. These experimental findings revealed the existence of various flow patterns, which are dependent on the characteristics of the impinging jet [
13]. Additionally, Ito et al. [
14] conducted experiments on jet impingement at elevated plate angles, systematically varying the position and inclination of the plate within the ranges of
and 60° ≤
θ ≤ 90°, respectively.
It is crucial to note that CFD analysis effectively replicates the surface pressure distribution over the plate with the desired spatial resolution. The precise assessment of the pressure field greatly enhances flow diagnostics and, thus, industrial design, particularly in predicting and mitigating pressure loading and improving system reliability, which is especially vital for aerospace and defense applications. For example, considering the particular application to the takeoff of an aircraft (or a missile), the local increase in velocity results in a decrease in local pressure near the impingement area, which can lead to diminished lift. This reduction in lift may result in either a failed takeoff or increased fuel consumption. In recent years, several experimental studies have been conducted to explore the lift loss induced by this phenomenon, utilizing a simplified test configuration that includes a baffle plate surrounding a single air jet. For instance, a series of experiments examining various nozzle-to-ground heights were presented in [
15], indicating that a closer proximity to the ground surface correlates with a greater loss of lift.
In this demanding context, the primary objective of this work is to conduct a computational investigation of the impinging jet flow field utilizing a modern turbulence model. Unlike earlier numerical investigations, such as those presented in [
16], which relied on conventional unsteady Reynolds-averaged Navier–Stokes (RANS) models, this research also employs a hybrid turbulence modeling procedure known as detached-eddy simulation (DES). Indeed, traditional RANS simulations are limited in their ability to capture spectral content in the flow field, even when sufficient spatio-temporal resolution is achieved, as they compute the (statistical) mean flow, thereby eliminating all turbulence information from the resolved flow field through the averaging process. Although multi-scale RANS models that can accurately predict flow separation details and anisotropic turbulence are still in development, e.g., [
17], the growing accessibility of high-performance computing (HPC) resources for industrial research enables more precise evaluation of turbulent supersonic jet impingement. Specifically, the current DES formulation, which is based on the hybridization of the two-equation eddy-viscosity
shear-stress transport (SST) model [
18], akin to large-eddy simulation (LES), facilitates the emergence of a resolved turbulent spectrum, allowing for partial resolution of turbulent fluctuations. To examine the capability of the two different techniques to reproduce the complex flow physics arising from the jet impingement, various computations are conducted for varying inclination angles.
Furthermore, contrasting with similar investigations where commercial proprietary software is usually employed, the use of an open-source CFD solver allows us to overcome a couple of significant limitations. In fact, these solvers do not provide access to the underlying code, so users cannot modify the algorithms or scrutinize the underlying physical-mathematical model. Moreover, even a basic-level solver license can be prohibitively expensive, restricting the possibility of performing high-fidelity simulations. The ability of the selected CFD solver to correctly simulate wall-confined supersonic flows is tested by comparing the present results against both numerical data obtained with a commercial code, and experimental findings. Particular attention is devoted to the very complex flow structures created by the jet impact on the inclined plate, and the consequent pressure distribution on the solid surface, while examining the effect of varying the plate inclination.
The remainder of the paper is organized as follows.
Section 2 introduces the overall computational model, including the meshing procedure and the main numerical settings, while the two different turbulence models employed are given in
Section 3. The results of the various calculations are presented and discussed in
Section 4. Finally,
Section 5 provides some conclusions and perspectives of this study.
3. Turbulence Modeling
Two different turbulence modeling procedures were employed to assess the predictive capabilities of the present CFD model to simulate the jet flow impingement on the inclined flat plate, namely, unsteady RANS and DES. The compressible flow governing equations, based on the
k-
formulation, can be formally written the same way for both approaches [
27]. Specifically, the balance equations for mass, momentum, and energy read:
where
,
, and
represent the balanced variables per unit volume. The energy Equation (
3) ensures conservation of the total energy that is
, including internal energy
e, kinetic-energy, and turbulent kinetic-energy. Note, that in the present formulation, the energy transport equation is based on the enthalpy variable
and, once this equation is solved, the temperature field can be directly estimated.
In the above equations,
and
are the viscous and turbulent stresses,
the turbulent eddy-viscosity,
and
the molecular and turbulent Prandtl numbers. Here, according to the eddy-viscosity assumption, the unknown turbulent stresses are approximated as:
where
is the resolved strain-rate tensor, and
stands for the Kronecker delta, with
being determined in terms of the resolved turbulence variables.
3.1. Unsteady RANS Simulation
In this work, the high Reynolds number version of the
k-
turbulence model proposed by Wilcox [
27] was utilized. The turbulence modeling procedure involves the solution of two additional evolution equations for the turbulence variables
k and
, specifically:
Compared to the original version [
18], such a model variant uses the additional cross-diffusion term at the right-hand side of the
equation, as well as a stress limiter. These extra terms are applied using switches so that the model formulation maintains its simplicity, without the need for blending functions. In the above equations, the production of turbulent energy is represented by
, whereas the corresponding dissipation rate is given by
, while
. Moreover,
,
,
,
,
, and
are constant coefficients for which standard values were selected. For a comprehensive presentation of this two-equation eddy-viscosity model and relevant discussion, the interested reader is referred to [
28].
3.2. Detached-Eddy Simulation
The DES method represents a more sophisticated hybrid turbulence modeling formulation that can be adopted in conjunction with various RANS models. In this work, the
SST-based DES model is actually employed. The
SST model was among the different pure RANS models tested in the previous CFD study used as reference [
16]. Therein, a correct surface pressure distribution for the impinging jet flow at low incidence (
θ = 30°) was attained, due to the proven ability of this model to provide accurate mean flow results, especially for wall-bounded separated flows [
18].
In the current approach, the evolution equations for the two turbulence variables are slightly modified as follows:
with
given by:
where the switching functions
and
are the same as in the standard
SST model [
18]. In the above definition,
represents the resolved turbulence length scale,
is the largest (linear) dimension of the FV cell, with
standing for a calibration coefficient.
In practice, RANS and LES solutions are coupled by introducing a hybrid turbulent length scale, which is modified with respect to the pure RANS approach, where it holds
. This way, a single turbulence closure acts as either a sub-grid scale model, where the local grid resolution is fine enough for LES, or a RANS model, where that is not true. The two different regions are actually separated by a transition zone where the eddy-viscosity level is sufficiently low to allow the rapid regeneration of turbulent eddies. Basically, DES aims at the prediction of separated flows in fluids engineering research at a manageable cost. More details about the DES formulation can be found in the original work by Strelets [
29].
4. Results
As a preliminary step, the grid convergence analysis presented in
Section 2.2 is demonstrated in
Figure 3 by comparing the pressure distributions on the plate mid-section for the different mesh resolutions with each other, as well as against the experimental data in [
12]. The three different FV grids are made of 0.66 (coarse), 2.7 (medium), and 6.8 (fine) million cells for the whole computational domain. In this plot, the spatial coordinate along the plate is non-dimensionalized by the nozzle exit radius (
), whereas the pressure variable is non-dimensionalized by the nozzle chamber (stagnation) pressure that is
kPa. By inspection of these data, the numerical resolution corresponding to the finest grid can be regarded, a posteriori, as adequate for capturing most of the complex flow features. Therefore, the various calculations for the two different turbulence modeling approaches were conducted using this particular grid.
The pressure contours over the plate surface predicted by the two different methods are presented in
Figure 4. These results have been obtained once the statistically steady state is reached, for different inclination angles that are
= 30°, 40°, and 50°, respectively. Note, that practically since the near-wall region is modeled using the RANS approach for both turbulence models the shape of the pressure pattern of impingement is pretty similar for the two solutions, particularly at low plate incidences. As the inclination angle increases, slight differences start appearing. In fact, the higher the inclination angle, the stronger the jet deviation imposed by the plate (where
θ = 90° would correspond to perpendicular jet impingement), and the more complex is the separated compressible flow field. This way, the LES capability of the hybrid model leads to a better description of the resolved pressure field with increasing plate incidence. Apparently, at higher inclinations, the DES solution shows a more elongated impingement zone over the plate, with the stagnation flow region being slightly larger, compared to RANS.
The present qualitative results show a good agreement with reference numerical data in [
16], which were obtained using a commercial code. Expectedly, since the plate center is perfectly aligned to the nozzle axis, the maximum pressure level is achieved at the mid-section of the plate, as illustrated in
Figure 4. This way, when performing a quantitative analysis and comparison of pressure load using different methods, the plate section at
is selected.
Figure 5 shows the pressure distribution on the mid-section of the plate for the two different solutions, and varying plate incidence. Making a comparison with reference experimental data [
12], both models result in quite correctly reproducing the pressure peak and its location on the plate. At the lowest inclination angle, the numerical solutions are able to predict the pressure oscillations away from the high-pressure region, even though the values and locations of local extrema are not well captured.
The two different CFD solutions are further analyzed by looking at the resolved velocity and temperature fields, for varying plate inclination.
Figure 6 shows the normalized speed contours in the mid-plane
for the space region between the nozzle exit and the plate surface, while the corresponding thermal contours are depicted in
Figure 7. In any case, a stagnation zone is visible in the impingement region, together with the consequent boundary layers development in both directions along the solid surface. It is worth noting that both solutions are able to capture the jet boundary and the jet shock away from the wall with acceptably high fidelity. However, a slower and more consistent velocity field along the plate is obtained with DES, as the RANS computation seems to overpredict the fluid velocity in the wall jet region. Moreover, the DES approach is able to more clearly reveal the presence of the expansion waves that are produced at the exit of the under-expanded nozzle and their reflection from the jet boundary. Such waves turn into compression waves that coalesce to form a barrel-shaped shock. Increasing the plate angle, in accordance with the experimental findings, the flow structure changes and the flow complexity grows. Apparently, as the jet spreads along the plate, the fluid velocity tends to decrease, while the boundary layer thickness increases. Moreover, the stagnation zone tends to shift slightly upwards. This flow scenario is in agreement with previous numerical investigations, where more sophisticated, and more costly, LES techniques were adopted [
30]. However, the drift of the stagnation bubble should be further investigated, due to its crucial impact on the evaluation of heat transfer and wall shear stress over the plate. By inspection of the temperature contours in
Figure 7, it seems that unsteady RANS tends to over-predict the temperature in the jet core region, as well as in the wall proximity, while the thermal field exhibits a more uniform behavior for DES. Moreover, the temperature in the far-field region appears slightly higher for the RANS solution.
In accordance with previous investigations, the pure RANS solution is shown to provide too slow decay of the mean velocity profile along the jet axis, including the core region, which can be attributed to the underprediction of the turbulence mixing in the shear layer. The nearby temperature field is overestimated, leading to the misrepresentation of the heat transfer rate in the impingement zone. On the contrary, the hybrid RANS/LES approach produces faster mixing, due to the capability to take into account the vortex breakup process in the shear layer of the jet and in the near-wall region, which leads to a more correct estimation of the thermal field in the proximity of the plate [
31].
A more detailed description of the complex flow structures arising from the jet impingement is presented in the following, using the numerical Schlieren technique. In
Figure 8, the various flow scenarios associated with different plate inclinations are visualized by reporting the contour maps of the density gradient modulus, making a direct comparison against experimental Schlieren images provided by Nakai et al. [
12]. For the numerical data, the contour level varies between zero (white color) and 0.053 (yellow color). Indeed, based on the experimental findings, it was found that the complex flow patterns can be classified into different types, namely: Type I (with/without bubble), Type II (with/without bubble), and Type III, as discussed in [
13,
14]. The actual flow type can be predicted in advance, once the plate angle, the nozzle-to-plate distance, and the shock cell length are given and known. Apparently, at
θ = 30°, both solutions correctly resemble the experimental flow field, with no bubble being present. In this case, the flow can be classified as Type III. A good agreement with experiments is also found at
θ = 40°, where the plate shock forming at the upper side of the plate is predicted by both solutions, confirming what was found in [
16]. According to these results, the flow can be classified as Type II (with bubble). Note, however, that the plate shock size seems thicker for DES than for unsteady RANS. Lastly, at
θ = 50°, the flow type is still of Type II (with bubble). In fact, the intermediate tail shock and the stagnation bubble are clearly recognized in the present Schlieren pictures, consistent with the pressure peak reported in
Figure 5. The presence of the bubble on the upper part of the impingement zone is more evident here than for previous numerical investigations [
16]. Note, that the current results are fully consistent with the jet flow classifications reported in [
14].
Overall, both approaches seem to correctly capture the essential flow structures, compared to experiments. However, DES is shown to predict some fluid flow structures more accurately, both near and away from the plate. For example, jet boundary and intermediate tail shock are captured with high detail by the hybrid model, whereas the traditional RANS model shows a lower level of resolution. Importantly, enhanced turbulence resolution is attained with DES by utilizing the same computational grid and, thus, with a limited increase in computational complexity.
The DES approach, which partly resolves the turbulent fluctuations, is able to provide an estimation for the turbulent Reynolds stresses. A qualitative analysis of the normal stresses
and
is presented in
Figure 9, where the stress level is non-dimensionalized by the inlet velocity squared
. The axial component is higher in magnitude, in agreement with previous numerical and experimental investigations [
32,
33]. This fact can be attributed to the compressible effects at the exit of the under-expanded supersonic nozzle. Moreover, the amplitude of the velocity fluctuations in the axial direction tends to increase along the plate, as the inclination gets higher. Practically, a rise in the plate incidence leads to an increase in the level of disturbance for the impinging jet flow, resulting in a higher magnitude for the transversal component
along the plate.
5. Conclusions and Perspectives
The open-source software OpenFOAM was utilized for the computational evaluation of a turbulent supersonic jet impinging on an inclined flat plate at various incidences. The proposed applied CFD model was based on the unsteady RANS approach using wall functions, thus avoiding excessive resolution requirements in the wall region. In addition, the modern DES technique, which integrates RANS and LES methods, was tested for this particular application in fluids engineering, with the present results being validated against reference numerical data and experimental findings. According to the analyses carried out, both qualitatively and quantitatively, the current models are shown to correctly reproduce the surface pressure distribution on the inclined plate, depending on the plate inclination. Although more sophisticated turbulence models can be used, this study demonstrates that simpler RANS-based models provide a way to create computationally less expensive simulations, where the mean flow can be predicted with engineering accuracy and moderate computational complexity. However, the DES method was able to better reproduce the very complex structure of the impinging jet flow, combining RANS-like turbulence resolution in the wall proximity with partial capture of the turbulent fluctuations away from the plate.
As far as perspectives are concerned, the deeper analysis of the temperature field will be the next step of this continuing research, in order to thoroughly investigate the heat transfer at the plate. The importance of such an analysis is actually twofold. On one hand, in comparison to other conventional systems, jet impingement cooling has proven to attain higher heat transfer rates. In this context, the computational study contributes to determining the jet flow conditions (including the geometry) ensuring the desired thermal energy drain. On the other hand, in a number of industrial applications, the increasing temperature of the plate has to be controlled, to avoid the undesired deterioration of the thermophysical properties of the solid surface coating. Further investigations will be performed by testing innovative unsteady RANS modeling approaches based on adaptive numerical methods [
34], along with enhanced DES formulations, like delayed detached-eddy simulation (DDES) [
1,
35]. Future research may also focus on applied aeroacoustics, specifically examining the acoustic waves generated by fluid-structure interaction, which can be estimated using the pressure fluctuations predicted by the LES component of the hybrid solution.