Numerical Experiments of Subsonic Jet Flow Simulations Using RANS with OpenFOAM

This study investigates the accuracy of different numerical schemes of OpenFOAM software to simulate compressible turbulent jets. Both pressure-based schemes utilizing the implicit PIMPLE algorithm and density-based schemes relying on AUSM scheme and explicit Runge-Kutta time integration are considered. The results of the numerical tests are compared and validated against data from NASA ARN nozzle geometry. The choice of parameter setting of the schemes is discussed in depth and possible optimization strategies are proposed to increase accuracy of RANS simulations of turbulent jets.


Introduction
The numerical simulation of compressible jet flows is relevant to many areas of aerospace engineering, such as for example, the identification of noise sources in jet aeroacoustics [1] or the prediction of exhaust dispersion in the atmosphere [2]. Modeling high speed jet flow requires a compressible flow solver of Navier-Stokes equations. If Direct Numerical Simulations (DNS) of jet flows in real airframe and operating conditions is still out of reach for current computational resources-and it will be so in a foreseeable future-large-eddy simulations (LES) of jets in single or double stream flows including nozzles appeared recently in the literature [3]. LES offers a promising approach especially when the focus is on the unsteadiness of the jet. However, their computational cost can be extremely high due to the resolution required (both in space and time) to capture the turbulent features of the flow, and the large computational domains needed to compute the entire 3D flow-field even though the flow is statistically axisymmetric (at least for turbulent round jets). Indeed, recent reviews of experimental and numerical studies of turbulent round jets showed that DNS and LES satisfy the desired level of accuracy but the computational cost can be extremely high [3] [4].
RANS offers a valid alternative to LES when the interest is on the analysis of mean flow characteristics in complex geometries and airframes that are representative of actual engine nozzles, with or without mounted pylon and wings. In RANS models, turbulence is modeled via transport equations for Reynolds stresses in time-averaged Navier-Stokes equations. Spalart-Allmaras [5] one-equation model and k-ε [6], and k-ωSST [7] two-equation model are among the most turbulence models and have been used extensively for incompressible and compressible jets, in subsonic or supersonic conditions [8].
Even though RANS modeling significantly reduces the simulation computational cost, jet flow prediction accuracy using RANS models is limited [9] [10].
Miltner et al. (2015) performed a comprehensive analysis of turbulent, free jets modeled with several RANS models and compare the results with experimental data measured by Laser-Doppler-Anemometry for three-dimensional flow [11].
Their findings suggest that, while some models can offer higher accuracy for a specific application such as k-ωSST for free jet flow, they fail to provide the same level of accuracy for more complex problems. RANS modeling is more likely to exhibit limitations in the cases of relatively high temperature jets, when the flow cannot be treated as incompressible or when the axisymmetric conditions do not hold [12].
For applications such as aeroacoustics, accurate turbulence properties are required to assess jet noise emissions. , modified three different two-equation stress transport RANS models, tailored specifically for jet flow simulations. Their model modifications succeeded in predicting mean axial velocity with higher accuracy, but no improvement for the prediction of kinetic energy was detected [13]. In addition to the accuracy, the uncertainty of the turbulent models affects the simulation reliability. Mirsha and Iaccarino (2017), investigated discrepency between the RANS simulated data and high precision data to estimate models' uncertainty [14]. They investigated the uncertainty for subsonic and supersonic, and hot and cold jet flow simulations for four different nozzle geometries.
An advantage of simulating a jet flow in axisymmetric round jet with RANS models is that the computational mesh can be reduced by exploiting the axi-symmetry of the problem. Georgiadis et al (2006) used a quadrant of nozzle to simulate lobed nozzle flow fields due to symmetry in both x-y and x-z planes [15].
Depending on the symmetry level, the grid can be reduced further to a two dimensional plane.  [20]. They also tried to improve the numerical stability of central scheme by adopting an alternative form of the convective terms that conserves kinetic energy. Li and Paoli [21] analyzed the scalability of rhoCentrlaFoam and developed an explicit version of the scheme called rhoCentrlaRK4Foam that is based on Runge-Kutta temporal integration.
The main objective of this work is to investigate the accuracy of two schemes of OpenFOAM for RANS simulations of turbulent compressible jets. Specifically, this study utilizes rhoPimpleFoam, a pressure based solver from the Open-FOAM baseline distribution, and rhoEnergyFoam, a density based solver which has been developed by Modesti et al. [22]. The benchmark case for the validation is the NASA Acoustic Reference Nozzle (ARF) geometry [23]. The paper is organized as follows: Section 2 gives an overview of the numerical schemes and the computational setup while Section 3 presents the results of the numerical experiments. Finally, Sections 4 and 5 summarize the discussions and draw some conclusion of the study.

Numerical Methods
In the finite volume approach, the Navier-Stokes equations for a compressible ideal gas are integrated over an arbitrary control volume V: (1) where , i u f and v i f represent the conservative variables vector and the associated Eulerian and viscous fluxes, respectively. 0 , , (2) In OpenFOAM, the control volume is the volume of the grid cell, and the solvers differ in the way the numerical fluxes are discretized and the grid-average equations are evolved in time.

Compressible Solvers: rhoPimpleFoam
RhoPimpleFoam is a transient solver for turbulent compressible flow. It uses the 1) The inner-corrector sets the number of times the pressure is corrected in each iteration.
2) The outer-corrector, specifies the number of iterations the system of equa-

Compressible Solvers: rhoEnergyFoam
RhoEnergyFoam is a density-based solver of transient compressible Navier-Stokes equations. It is designed obtain low diffusive solutions in complex geometries [22]. With local use of the AUSM flux splitting method, rhoEnergyFoam preserves total kinetic energy and maintains properties conservation on unstructured mesh. In addition, the AUSM diffusive flux provides the solver with shock-capturing capabilities [24]. The Eulerian fluxes are split into convective and pressure terms.
is the numerical flux calculated at the interface of two adjacent cells, denoted as Owner and Neighbor, respectively (see Figure 1). Each term of the Eulerian flux consists of a central part and a diffusive part.
The central part of the convective flux is interpolated with the midpoint scheme, whereas for the pressure flux, the central part is interpolated with the standard central scheme. The midpoint scheme allows the solver to preserve the flow's kinetic energy. The convective and pressure central fluxes are calculated as follows:

Computational Grid
The axisymmetric nature of the RANS solution for round jets allows the geometry to be reduced to a two-dimensional grid. In order to replicate a two dimensional     Table 2 shows the grid details for each mesh.

Simulation Setup
The subsonic axisymmetric jet flow exits into quiescent ambient air with 0.

Results
Section 3.1 discusses a series of tests conducted using rhoPimpleFoam with different turbulence models. Then the grid dependency test is conducted using this pressure-based solver. In Section 3.2, the medium size mesh and the k-ωSST

rhoPimpleFoam Solver
The Spalart-Allmaras, k-ωSST and k-ε turbulence models are used for the 2D  Figure 4(a) and Figure 4(b) show the streamwise contours of velocity and kinetic energy for the medium grid resolution using the k-ε turbulence model.

Figures 5(a)-(c) report the axial velocity along the centerline normalized by
jet exit velocity for different grid resolutions. The Spalart-Allmaras, k-ωSST, and k-ε turbulence models are used for this test. It can be seen that for all three turbulence models, the medium and fine grids show no significant difference, which indicates that the medium grid is sufficiently resolved for the RANS simulations.  The decay error for each turbulence model is calculated using the root mean square error, RMSE, for k-ωSST, k-ε, and Spalart-Allmaras are 3%, 3.6%, and 9.1% respectively. Figure 6 shows the normalized radial velocity, axial velocity, and kinetic energy profiles in the radial direction at 5D j , 10D j , 15D j , and 20D j axial distances from the jet exit. Figure 6(a) shows that the peak radial velocity occurs farer from the centerline compared to the experimental data, and the radial difference between peaks becomes wider when moving away from the nozzle exit. In addition, Figure 6(a) shows that, as the axial distance increases, the simulations fail to accurately predict the location where the radial velocity changes sign (from outward to inward). with experimental data. Along the streamwise direction, the k-ε and k-ωSST results tend to underpredict the velocity at the centerline more, whereas they reach a better agreement with the experimental data between 1 to 1.5 diameter distances from the centerline. Moreover, all three models represent the experimental data accurately at a distance of 5D j from the jet exit. The accuracy deteriorates when moving further away from the exit for k-ωSST and Spalart-Allmaras models, which is in line with the trend observed in Figure 5.

rhoEnergyFoam Solver
Based on the grid dependency test results in Section 3.1, the medium resolution grid is sufficient to insure grid convergence and is then used for the simulations  Table 3 presents the potential core length and velocity decay errors with respect to the experimental data for all tests conducted in this section. Figure 7 compares the results obtained using rhoPimpleFoam and rhoEner-gyFoam with the NASA subsonic cold jet experimental data. The figure shows that the mode-C setting of rhoEnergyFoam is the most dissipative scheme and produces the highest decay rate in the streamwise direction compared to the experimental data. The first step to reduce the numerical dissipation is by switching to mode-B configuration. However, as anticipated in Sec. 2.2, this leads to higher oscillations in the core region as shown in Figure 7. The overprediction of the potential core length halves by simply disabling the convective diffusive flux (see Table 3).
The next approach to reduce the generated numerical diffusion of rhoEner-gyFoam solution is to alter the threshold value *     Figure 9 compares the centerline velocity in the streamwise direction using different models and smoothness factors. Increasing the smoothness tends to reduce the overprediction of the potential core length. The decay rate is also improved compared to the baseline mode-C (see Table 3).
With the aim of minimizing both numerical dissipation and oscillations, we finally considered different combination of the threshold and smoothness factor.
The results of these tests showed that the values * 0.1 θ = and 0.9 smooth = represent an optimal combination in that they improve the decay rate substantially. Table 3 shows that with this combination of parameters, the error in the decay rate was reduced to 2.9%, the lowest among all tests performed in the study. However, this combination cannot completely eliminate the oscillations in the potential core region. The same test was then repeated by reducing the time-step by one order of magnitude, which leads to a substantial damping of numerical oscillations as shown in Figure 10.     increase from the centerline. This behavior is expected to appear when using the K-ωSST turbulence model which was also observed in Figure 6. The location where radial velocity peaks was well captured at all four axial locations. It also worth noting that while all simulated cases predicted the change of velocity vector direction at 5 j x D = accurately, they gradually overpredict it further away from the jet exit.

Discussion
In Section 3.1, the solution's dependency on the grid resolution was investigated.
The grid dependency test showed that the medium size grid has the necessary and sufficient resolution for the solution to be independent of the grid size.
Based on the centerline velocity comparison in Figure 5(d), the choice of turbulence model depends on whether the study prioritizes a more accurate prediction of the potential core length or velocity decay rate. A closer look at the normalized centerline velocity in Figure 10 shows that the density-based solver tends to generate oscillation even when running in mode C.
It appears that the slightest change in the runtime configuration can magnify any existing oscillation. To overcome this stability issue, the time-step was lo-

Conclusion
This