Improved Supersonic Turbulent Flow Characteristics Using Non-Linear Eddy Viscosity Relation in RANS and HPC-Enabled LES

Abstract

A majority of the eddy viscosity models for supersonic turbulent flow are based on linear relationship between Reynolds stresses and mean strain rate. The validity of these models can be improved by introducing non-linearity in relation as RANS models offer advantages in terms of reduced turnaround times typical of industry applications. With these benefits, the present work utilizes quadratic constitutive relation (QCR) with Menter’s k omega SST model to characterize the flowfield of rectangular jets. The sensitivity of this model with QCR, weighted towards diffusion, dissipation, and a combination of both, is addressed. Viscous large eddy simulations (LES) with WALE subgrid scale models are employed for qualitative comparisons using a commercial solver. Massively parallel LES are enabled by the new in-house 1088-core computing cluster at the University of Cincinnati and are also used for benchmarking. The nearfield results are validated with available experimental data and show good agreement in both fidelities. Flow characteristics, including the shear layer profiles, Reynolds stresses, and turbulence kinetic energy (TKE) and its production are compared. LES reveal higher TKE production in the regions with highest Reynolds stresses. It is comparatively lower in QCR RANS. As a special case of TKE analysis in jets, a preliminary investigation of retropropulsion is outlined for rectangular nozzles for the first time. Improved flow behavior by implementation of a non-linear relationship between Reynolds stresses and mean strain rate is demonstrated.

1. Introduction

Turbulence is the most fascinating yet complex phenomenon that we encounter in our daily interaction with different types of fluids. Simulation of turbulent flows has been a topic of research for many years. While LES and DNS are the most computationally expensive techniques since they involve ‘less’ modeling and ‘more’ resolved flow, RANS models have become common in many industries due to their reduced turnaround time. One- and two-equation eddy viscosity RANS models are based on Boussinesq approximation, which assumes that the Reynolds stress tensor is linearly proportional to the mean strain rate and is correlated by turbulent eddy viscosity [1]. The validity of this hypothesis has been examined by Schmitt [1] to demonstrate its limitations. Therefore, the present work attempts to consider non-linear eddy viscosity relation through multi-fidelity simulations for rectangular jets. High-speed flows are dominant in aerospace applications where turbulence becomes complicated due to compressibility effects. One such example is jets emanating from exhaust nozzles. Their flow structure is of particular importance due to shear-layer mixing, where the mean flow contributes to the turbulent fluctuations. Traditionally, axisymmetric jets have been studied exhaustively due to their wide variety of applications in many industries. With integrated aircraft-body concepts, non-axisymmetric jets have gained attention as they offer frame-integration benefits, along with reduced noise characteristics. From the early 1960s until now, many researchers have contributed to the fundamental understanding of turbulent jet mixing. Sforza et al. [2] experimentally investigated turbulent 3D jets with different orifice shapes, including rectangular, triangular, elliptical, and axisymmetric geometries, and concluded that they decay to axisymmetry far downstream of the jet exit. Some notable experimental work on rectangular jets include the studies reported in references [3,4,5,6,7,8]. Zaman [7] highlighted the significance of vorticity dynamics in the development of rectangular jets. While these experimental works aided the understanding of flow physics, some researchers also approached the problem computationally. One such early work is by Bobba and Ghia [9]. They calculated the flow field of a rectangular jet using a two-equation turbulence model and noted that the conditions at the nozzle exit play a dominant role in flow development downstream. In late 1990s, Moin and Kim [10] emphasized the use of supercomputers for tackling turbulence in various applications ranging from automotive to aerospace. With significant advancements in compute power over the last ~20 years, many researchers have conducted computationally intensive large eddy simulations to investigate the aerothermodynamics and acoustics of compressible jets [11,12,13,14]. These studies utilized various strategies for LES, such as the compressible Euler solver reported in [13]. Meanwhile, Bellan [14] used the dynamic Smagorinsky model to achieve closure. Others [15,16] reported the experimental results on flow and noise characteristics of rectangular jets. Bridges and Wernet [17] experimentally investigated the turbulence statistics associated with rectangular jets at low to high aspect ratios, which provided some insights on crossflow turbulence.

The authors’ previous work on this topic addressed the nozzle aspect ratio effects and the role of multiphysics modeling in design cycles [18,19,20,21] that used the k-omega SST turbulence model with Boussinesq approximation. The model assumes the tensor of turbulent stresses is proportional to the mean strain-rate tensor [22]. This assumption has its own deficiencies, as discussed by Wilcox in his famous book on turbulence modeling [22]. To circumvent these deficiencies, non-linear constitutive relations have been suggested in the literature [22]. Among such relations, quadratic constitutive relation is employed in the present work which was first suggested by Spalart [23]. The motivation is to examine the turbulence statistics associated with rectangular supersonic jets using LES and RANS on the common grounds of non-linear eddy viscosity models. Some researchers have attempted to characterize the turbulence in a tip-leakage flow of a single blade using non-linear RANS and zonal LES, as shown by Monier et al. [24]. While Chen et al. [25,26] demonstrated the vorticity dynamics-based flow diagnosis for a 1.5 stage high-pressure compressor using RANS simulations to improve compressor aerodynamic performance, Wu and Porté-Agel [27] presented the atmospheric turbulence effects on wind turbine wakes using LES. While these studies employed LES, Siddappaji [28] demonstrated TKE is equal to viscous dissipation using momentum, vorticity, and entropy transport through RANS simulations in adiabatic process with an in-house analysis framework. Although these references provided some background on turbulence characterization for various applications, the same has not yet been shown for rectangular jets using such models. Therefore, the goal is to use verified and validated high-fidelity simulations to compare with RANS. Note that the primary focus is to capture the non-linearity in the relationship between Reynolds stress and mean strain rate through multi-fidelity simulations. Since the approach is general, preliminary analysis of retro propulsive flow physics is also discussed briefly in Section 3.4.

The nozzle geometry is shown in Figure 1. It has a rectangular cross-section with aspect ratio (AR) = 2, equivalent diameter = 20.65 mm, and design Mach number = 1.5. The sharp throat causes a shock wave just downstream of the throat. The remaining of paper is outlined as follows. The methodology section details the non-linear RANS and LES simulations, the sensitivity analysis of the k-omega SST model towards diffusion and dissipation, the verification and validation and kinetic energy spectra for the LES case, and the HPC runtime statistics. The subsequent sections describe the results, primarily focusing on turbulence characterization and TKE. The final subsection briefly describes the retropropulsion flow physics of rectangular nozzles as a preliminary analysis.

2. Methodology

The present work utilizes RANS with quadratic constitutive relations (QCR), which is called QCR RANS from here onwards. A full-fledged LES with WALE (wall-adapting local eddy viscosity) subgrid scale model is used for flow diagnosis and comparing with QCR RANS solution. We focus on the generality of this approach.

2.1. RANS Details

3D steady RANS equations are solved for compressible flow of air, modeled as an ideal gas that uses Menter’s k-omega SST turbulence model with QCR. All simulations are conducted using the commercial solver Simcenter Star-CCM+ 15.04.008-R8 [29]. A coupled implicit solver with second-order spatial discretization is used. All simulations correspond to design operating conditions i.e., nozzle pressure ratio (NPR) = 3.67, with Mach number = 1.5 and unheated jet. NPR is defined as the ratio of total pressure at nozzle inlet to freestream pressure. Nozzle inlet is modeled as stagnation inlet, nozzle walls are modeled as no-slip adiabatic walls, and domain boundaries are modeled as freestream, representing atmospheric pressure and temperature. The computational domain extends 100xDe downstream of the nozzle exit, 15xDe radially, and 10xDe upstream. To capture the shock-wave formation and the plume, volumetric refinements are used for polyhedral mesh. A total of 15 prism layers with expansion factor of 1.5 are used for the boundary layer. Wall, y+, is maintained ~ 1 along the nozzle walls.

In the k-omega SST model, the linear relation between the Reynolds stresses and the mean strain rate tends to strongly underpredict the anisotropy of turbulence [29]. To overcome this issue, non-linear constitutive relations have been suggested, and a brief account can be found in reference [22]. In this work, the quadratic constitutive relation suggested by Spalart [23] is used since it accounts for the non-linear turbulence production without introducing additional partial differential equations.

2.2. LES Details

LES with WALE subgrid scale model is conducted as it accounts for both strain and rotation tensor [30]. A compressible viscous solver that solves the energy equation is used. More details on LES in Star-CCM+ can be found in reference [30]. Nozzle inlet is specified as stagnation inlet corresponding to NPR = 3.67 and unheated jet. The filtered governing equations of mass, momentum, and energy are given by

$$\frac{\partial \rho }{\partial t}+\nabla .\left(\rho \tilde{v}\right)=0$$

(1)

$$\frac{\partial \left(\rho \tilde{v}\right)}{\partial t}+\nabla .\left(\rho \tilde{v}\otimes \tilde{v}\right)=-\nabla .\tilde{p}I+\nabla .(\tilde{T}+T_{SGS})+f_b$$

(2)

$$\frac{\partial (\rho \tilde{E})}{\partial t}+\nabla .(\rho \tilde{E}\tilde{v})=-\nabla .\tilde{p}\tilde{v}+\nabla .(\tilde{T}+T_{SGS})\tilde{v}-\nabla .\tilde{q}+f_b\tilde{v}$$

(3)

In above equations, $\rho$ is the density, $\tilde{v}$ is the filtered velocity, $\tilde{p}$ is the filtered pressure, I is the identity tensor, $\tilde{T}$ is filtered stress tensor, $T_{SGS}$ is subgrid scale stress tensor, $f_b$ is resultant of body forces, $\tilde{E}$ is filtered energy per unit mass, and $\tilde{q}$ is the filtered heat flux.

Obtaining precursor converged RANS solution is crucial before launching LES case. Therefore, much attention was given to the setup of the RANS case, which captured the nozzle internal wall-bounded flow and jet region. Appropriate grid refinements are used in these areas to keep wall y+ ~1 along nozzle walls. This ensures a converged RANS with adequate quality of agreement of the results. The LES are run for a total time = 6 × (100 D_e/u_j) seconds, i.e., six flow-through times with fixed time step size = 10⁻⁶ s and constant CFL. Three volumetric grid refinement zones are used for nozzle, near-field, velocity decay region, and details are shown in

Table 1. Grid refinement details.

Mesh Size	Number of Cells (Million)	Zone I	Zone II	Zone III
Coarse	28	De/34	De/41	De/34
Medium	37	De/41	De/69	De/41
Fine	73	De/51	De/82	De/51

. Zone I corresponds to the refinement in the nozzle, and Zone II corresponds to the region spanning from nozzle exit to five diameters downstream, which covers the shear layers on both minor and major axis. Zone III starts immediately after the end of Zone II and lasts 25 diameters downstream. These zones are in cylindrical frustrum shape, in line with the direction of jet spread. Figure 2 shows the computational grid on minor and major-axis symmetry plane, illustrated in black and blue, respectively. The bottom image shows a close-up view of the grid on minor axis symmetry plane.

LES cases were run on the new in-house HPC cluster (described in Appendix A) at University of Cincinnati’s Advanced Research Computing (ARC) Center and are also used as a suite of load tests. The cases were run on a 17-node cluster with a total of 1088 cores. Note that this is the first time such heavy computations have been conducted with massive parallelization using an in-house cluster at University of Cincinnati. HPC statistics for coarse, medium, and fine-mesh LES cases are shown in Figure 3a,b. Although medium mesh is heavier in terms of number of cells, it took less time to complete than the coarse mesh case. We anticipate that this could be due to the cluster not taking 100% case load, as well as the recommended scaling for parallelization, which is ~50,000 cells per core for Star-CCM+ [29]. The total CPU time is calculated as—Number of cores × total wall-clock time.

It does not include the time taken for precursor RANS simulations, LES setup, cluster idle time, job restart time due to wall-clock time limit, and post-processing time. It includes the runtime for completion of six flow-through times only. The quantities of interest were automatically extracted after a certain number of time steps after the completion of ~2–3 flowthrough times to ensure the accuracy without effects from precursor simulation.

2.3. Sensitivity of K-Omega SST with QCR towards Diffusion and Dissipation

Menter’s k-omega SST model uses blending function, which is a hyperbolic tangent function that appears in the diffusion and dissipation terms. In free shear layers, jet mixing is governed by turbulence mechanisms where diffusion and dissipation play a key role. Therefore, their sensitivity is analyzed in this section to examine the effect on centerline velocity prediction. Equations (4) and (5) represent the transport equations of TKE ($k$) and specific dissipation rate ($\omega$). In below equations, $\overline{v}$ is the mean velocity, $\mu$ is the dynamic viscosity, ${\sigma }_k$ and ${\sigma }_{\omega }$ are the model coefficients, $P_k$ and $P_{\omega }$ are the production terms, $f_{{\beta }^{*}}$ is the free shear modification factor, $f_{\beta }$ is vortex stretching modification factor, and $S_k$ and $S_{\omega }$ are the user defined source terms. The first term on the right-hand side of Equations (4) and (5) represents diffusion, and the second-last term represents dissipation. Three variations of blending function are chosen, which are categorized as highly diffusive, low diffusive and a combination of diffusion and dissipation. Figure 4 shows the jet centerline velocity prediction for all. It can be seen that the slope of centerline velocity decay is sensitive to these terms, although the amplitude of velocity variations is unchanged until x/De ~ 7.

$$\frac{\partial }{\partial t}\left(\rho k\right)+\nabla .\left(\rho k\overline{v}\right)=\nabla .\left[\left(\mu +{\sigma }_k{\mu }_t\right)\nabla k\right]+P_k-\rho {\beta }^{*}{f}_{{\beta }^{*}}\left(\omega k-{\omega }_0{k}_0\right)+S_k$$

(4)

$$\frac{\partial }{\partial t}\left(\rho \omega \right)+\nabla .\left(\rho \omega \overline{v}\right)=\nabla .\left[\left(\mu +{\sigma }_{\omega }{\mu }_t\right)\nabla \omega \right]+P_{\omega }-\rho \beta f_{\beta }\left({\omega }^2-{{\omega }_0}^2\right)+S_{\omega }$$

(5)

2.4. Validation

Experimental data for jet centerline velocity and TKE along lip lines are reported in references [15,16], which are based on PIV measurements conducted at University of Cincinnati’s Gas Dynamics Propulsion Lab, and are used for validation. PIV is a non-intrusive measurement technique with minimal distortion to the flow, and therefore, is an obvious choice for the validation of results. Figure 5 shows that the coarse and medium mesh matches well at the initial shock cells, but further downstream, the effect of mesh resolution plays a significant role, and fine mesh LES captures the velocity decay quite well. The improvements in centerline velocity prediction can be seen by comparing baseline RANS case [19] with Boussinesq approximation, denoted as ‘baseline BSQ’ vs. QCR RANS. The baseline RANS case overpredicts the length of potential core, and therefore, the centerline velocity decay starts at a later location compared to experimental data. On the other hand, QCR RANS captures the trends exhibited by experimental data quite well, with diffusion and dissipation playing an important role, as discussed in the previous section. Figure 6a,b show the comparison of TKE along the minor and major-axis lip lines. On minor axis, LES underpredicts TKE magnitudes because the turbulence in the shear layer is underpredicted in LES as compared to experimental data. This is primarily because the boundary layer inside the nozzle develops to become turbulent until the nozzle throat, and after the throat, it undergoes shock-induced separation and reattachment. However, as one moves further downstream, the turbulence levels are in good agreement with experimental data. The TKE on major axis matches quite well with experimental data. Length of potential core in LES is also predicted well compared with experimental data from reference [16], which is, x/De = 7. This indicates that LES has captured the necessary trends exhibited in the experiments. QCR RANS overpredicts the TKE along minor-axis lip line, while it is adequate on major-axis lip line.

The TKE in experimental data [16] is defined as $TKE=\frac{1}{2}\left(\overline{u^{\prime }{u}^{\prime }}+2\overline{v^{\prime }{v}^{\prime }}\right)$, and therefore, this definition is used to plot the TKE in LES. The results are agreeable in quality and sufficient to move forward with the objective of the present study.

2.5. Kinetic Energy Spectrum in LES

Figure 7 shows the fast Fourier transform (FFT) of time history of kinetic energy at x/De = 20 along the jet centerline. It verifies that the present LES has captured the necessary length scales according to Kolmogorov energy spectrum. The time history is recorded after every 20 time steps. A Python-based tool was written to automatically extract the data and perform FFT. Please note that the spectrum was verified for all three cases, and the results for fine mesh are shown.

3. Results

The results are mainly focused on turbulent flow characterization, and the final section briefly describes the retropropulsion flow physics.

3.1. Turbulent Kinetic Energy

Turbulent fluctuations in the flow arise due to the inherent differences between minor and major-axis boundary layer to shear layer transformation. Therefore, the primary source of asymmetry is through the boundary-layer growth. TKE magnitudes are dependent on the definition used, meaning whether the three components of velocity fluctuations are considered or whether the definition assumes $\overline{v^{\prime }{v}^{\prime }}=\overline{w^{\prime }{w}^{\prime }}$. Computationally, it is possible to access all three components of TKE, while it may be a challenge in an experimental setup. The experimental data [16] were based on the assumption of $\overline{v^{\prime }{v}^{\prime }}=\overline{w^{\prime }{w}^{\prime }}$. The validity of this assumption is examined by highlighting the differences in both definitions in Figure 8a,b. They reveal the TKE magnitudes on the minor and major axis planes, as well as on the crossflow planes. The difference is dominant on the major axis plane, as can be seen in Figure 8b. This is because radial and spanwise components of TKE are not equal at all locations.

This is further confirmed by the Reynolds normal stresses shown in Figure 9. Minor-axis lip line stresses from LES are shown in Figure 9b. This indicates the dominance of the axial component of stresses, while radial and spanwise components are comparable. The major-axis lip line in Figure 9d shows the asymmetry in radial and spanwise components, which is also seen in the corresponding QCR RANS solution from Figure 9c. This explains the differences in TKE on the major-axis lip line seen in Figure 8b. In all plots, QCR RANS predicts the asymmetry of Reynolds stresses. Past 30 De, the Reynolds stresses damp out monotonically.

3.2. Shear Layer Profiles

Figure 10 compares the radial velocity profiles in QCR RANS and LES on the minor-axis plane at various streamwise locations. They are in good agreement with one another until x/De = 4.9. However, the potential core in LES ends at x/De = 7, and therefore, the profiles differ past that location. This difference is clearly visible from Figure 10c,d, indicating that QCR RANS predicts the radial profiles in close agreement with LES up to a few diameters downstream of the nozzle exit.

3.3. Turbulent Kinetic Energy Production and Reynolds Stresses

Production of turbulent kinetic energy is one of the main mechanisms of turbulence transport. As the jet exits the nozzle, the velocity fluctuations reach peak values. Recalling TKE transport Equation (4) in k-omega SST model, the first term on the right-hand side is the diffusion term, combining turbulence diffusion and molecular diffusion. $P_k$ is the production term and is the sum of turbulent production, $G_k$, buoyancy production, $G_b$, and non-linear production, $G_{nl}$. The third term on the right-hand side represents dissipation, and the fourth term is a user-defined source term. At RANS level, the Reynolds stress tensor is given as the sum of linear and non-linear terms, ${\mathit{T}}_{RANS}={\mathit{T}}_{RANS,L}+{\mathit{T}}_{RANS,NL}$. QCR, which is included in the production term, activates the antisymmetric normalized rotation tensor in addition to the linear part and accounts for the anisotropy of turbulence by adding non-linear functions of strain and vorticity tensors [29]. When QCR is activated in the solver, Reynolds stresses are internally calculated at RANS level. $G_k$ and ${\mathit{T}}_{RANS,NL}$ are given by

$$G_k={\mu }_t{f}_c{S}^2-\frac{2}{3}\rho k\nabla .\overline{v}-\frac{2}{3}{\mu }_t{\left(\nabla .\overline{v}\right)}^2$$

(6)

$${\mathit{T}}_{RANS,NL}=-2{\mu }_{t}0.04645\left(\mathit{O}.\mathit{S}-\mathit{S}.\mathit{O}\right)$$

(7)

where $\mathit{S}$ is the strain-rate tensor, ${\mu }_t$ is the turbulent eddy viscosity, and $\mathit{O}$ is normalized-rotation tensor [23,29]. Favre-averaged TKE transport equation [23] is given by

$$\overline{\rho }\frac{\partial k}{\partial t}+\overline{\rho }\widetilde{u_j}\frac{\partial k}{\partial x_j}=\overline{\rho }{\tau }_{ij}\frac{\partial \widetilde{u_i}}{\partial x_j}-\overline{t_{ji}\frac{\partial u_i^{″}}{\partial x_j}}+\frac{\partial }{\partial x_j}\left[\overline{t_{ji}{u}_i^{″}}-\overline{\rho u_j^{″}\frac{1}{2}{u}_i^{″}{u}_i^{″}}-\overline{p^{\prime }{u}_j^{″}}\right]-\overline{u_i^{″}}\frac{\partial P}{\partial x_i}+\overline{p^{\prime }\frac{\partial u_i^{″}}{\partial x_i}}$$

(8)

In the above equation, the first term on right-hand side is the production term, where ${\tau }_{ij}$is the Reynolds stress tensor. The second term is the viscous dissipation term. The third term represents the sum of molecular diffusion, turbulence diffusion, and pressure diffusion. The fourth term is the pressure work, and last term is pressure dilatation. LES simulations provide access to the turbulent quantities required to calculate the transport terms. Therefore, the production term is plotted using the definition from Equation (8). Note that minor and major-axis planes correspond to XY and XZ planes, and therefore, the Reynolds stress definition changes.

As seen in Figure 11, TKE production is the highest in LES in the region of highest Reynolds stresses on minor axis. Comparatively lower values of TKE production are seen in QCR. Additionally, it is steeper in LES, while gradual trends are seen in QCR in Figure 11a. Higher TKE production indicates more contribution from mean flow to turbulent fluctuations, which causes better mixing, and the potential core does not sustain longer in LES. Meanwhile, in QCR RANS, lower production values and a gradual rate cause the potential core to sustain longer. The production term reaches its peak on minor axis in LES where the Reynolds stresses are highest, as shown in Figure 12a,b.

TKE production obtained from LES is visualized in Figure 13, with isosurfaces of Q criteria colored by TKE production normalized by ${\rho }_j{u}_j^3$. Figure 14 shows the isosurfaces of Q criterion colored by streamwise vorticity, which are obtained from LES solution.

3.4. Retropropulsion Flow Physics—A Special Case Analysis

Retropropulsion is gaining momentum, although it has been a topic of research for last fifty years or so [31,32,33]. It is a phenomenon where engines are fired against an opposing freestream during descent in order to slow down. Soft landing, vertical takeoff, landing, and re-entry are some of the applications. SpaceX has recently demonstrated this technology in Earth’s atmosphere for their reusable Falcon 9 orbital rocket [32]. With the latest advancements in the space-exploration industry, retropropulsion has become crucial. This kind of technology is still relatively new and is useful not only on Earth but also on other planets for vertical takeoff, landing, and re-entry. Traditionally, the literature, as well as technology demonstrations, has employed circular nozzles for retropropulsion applications [33]. Since non-axisymmetric nozzles offer their own benefits, we feel the need to explore this area, and as a preliminary analysis, retropropulsion flow physics with rectangular nozzles is discussed in this section. Flowfield characterization is the main goal rather than a mission-specific application of retropropulsion. Note that the effect of angle of attack is not in the current scope. Three operating conditions are chosen to account for deceleration from the supersonic to subsonic regime:

• Supersonic nozzle-exit Mach and subsonic freestream;
• Subsonic nozzle-exit Mach and subsonic freestream:

  (a)

  high subsonic;

  (b)

  low subsonic.

The relationship between nozzle-exit thrust and freestream dynamic pressure can be quantified by thrust coefficient, which is defined as $C_T=\frac{F_T}{q_{\infty }A}$, where $F_T$ is the nozzle-exit thrust, $q_{\infty }$ is the dynamic pressure, and $A$ is the reference area.

Table 2. Operating conditions.

Case	Nozzle-Exit Mach Number	Freestream Mach Number	${\mathit{C}}_{\mathit{T}}$
a	1.6	0.95	1.24
b	0.7	0.55	0.74
c	0.7	0.3	1.65

describes the three cases with boundary conditions

Figure 15 shows the contour plots of density-gradient magnitude, defined as $\Vert \nabla \rho \Vert =\sqrt{{\left(\frac{\partial \rho }{\partial x}\right)}^2+{\left(\frac{\partial \rho }{\partial y}\right)}^2+{\left(\frac{\partial \rho }{\partial z}\right)}^2}$ and TKE normalized by $u_j^2$ for the three cases. For the first case, which is Figure 14a, $M_{exit}=1.6$ and $M_{\infty }=0.95$, a normal shock wave is present at location x/De = 1, and a recirculation zone occurs due to the opposing freestream. The TKE magnitudes increase in this zone, as seen in the corresponding TKE contour plots. For the second case with $M_{exit}=0.7$ and $M_{\infty }=0.55$, a similar flow structure can be seen. However, the recirculation zone and the normal shock wave are pushed further towards the nozzle exit. This is because the plume does not grow further due to an increase in opposing dynamic pressure. For the third case with $M_{exit}=0.7$ and $M_{\infty }=0.3$, the recirculation zone is larger. This is because the lower opposing dynamic pressure allows the plume to grow further downstream of the nozzle exit. The third scenario reveals an unstable flowfield as the plume shows flapping motion. In all three cases, the TKE is present in the recirculation zones, as opposed to a forward-propulsion case where the TKE is present along the shear layers/jet boundaries.

4. Conclusions

In this paper, non-linear eddy viscosity relation is employed in RANS and LES to characterize the turbulent flow in supersonic rectangular jets using a commercial solver. Menter’s k omega SST model with QCR is used. The sensitivity of k omega SST model with QCR weighted towards diffusion, dissipation, and a combination of both is discussed, which indicates that they play a key role in jet centerline velocity decay. Massively parallel compressible viscous LES are conducted on the new in-house 1088-core computing cluster at the University of Cincinnati. LES are also used for benchmarking the cluster. Kinetic energy spectrum is verified for LES that follow the Kolmogorov energy spectrum. Nearfield results are validated against the experimental data from available literature and are in good agreement for both fidelities. Flow characteristics, such as the Reynolds stresses and TKE and its production are analyzed. The asymmetry of turbulent kinetic energy is demonstrated through the TKE production term to show the differences in RANS and LES. The production is comparatively higher along minor axis plane in LES in the region where Reynolds stresses are highest. Relatively lower values of production and gradual rate cause the potential core to sustain longer in QCR RANS simulations. The assumption of $\overline{v^{\prime }{v}^{\prime }}=\overline{w^{\prime }{w}^{\prime }}$ is examined by using both definitions of TKE, i.e., the one with all three components and another with the assumption that radial and spanwise components are equal. This reveals that they are not equal at all locations. Reynolds normal stresses reveal asymmetry in both QCR RANS and LES. As a preliminary TKE analysis, retropropulsion flow physics is discussed for the first time concerning rectangular nozzles. Qualitative analysis of RANS and LES through non-linear eddy viscosity relation is demonstrated, which shows improved flow characteristics.