Turbulence Model Assessment in Compressible Flows around Complex Geometries with Unstructured Grids

Abstract

One of the key factors in simulating realistic wall-bounded flows at high Reynolds numbers is the selection of an appropriate turbulence model for the steady Reynolds Averaged Navier–Stokes equations (RANS) equations. In this investigation, the performance of several turbulence models was explored for the simulation of steady, compressible, turbulent flow on complex geometries (concave and convex surface curvatures) and unstructured grids. The turbulence models considered were the Spalart–Allmaras model, the Wilcox k-$\omega$ model and the Menter shear stress transport (SST) model. The FLITE3D flow solver was employed, which utilizes a stabilized finite volume method with discontinuity capturing. A numerical benchmarking of the different models was performed for classical Computational Fluid Dynamic (CFD) cases, such as supersonic flow over an isothermal flat plate, transonic flow over the RAE2822 airfoil, the ONERA M6 wing and a generic F15 aircraft configuration. Validation was performed by means of available experimental data from the literature as well as high spatial/temporal resolution Direct Numerical Simulation (DNS). For attached or mildly separated flows, the performance of all turbulence models was consistent. However, the contrary was observed in separated flows with recirculation zones. Particularly, the Menter SST model showed the best compromise between accurately describing the physics of the flow and numerical stability.

1. Introduction

Wall-bounded turbulent flows at high Reynolds numbers are mainly characterized by a wide range of time and length scales (Pope [1]). In Direct Numerical Simulation (DNS), the full Navier–Stokes (NS) equations are generally solved to capture all time scales and eddies in the flow. This is an accurate approach that supplies extensive information, but it demands the use of a very fine mesh, and significant computational resources, to properly capture features in the order of the Kolmogorov scales [2]. Filtered Navier–Stokes equations, with an additional sub-grid scale stress term, are employed in Large Eddy Simulation (LES), in which large scales or eddies are directly resolved and the effects of the small scales are modeled [3]. The LES approach still requires high mesh resolution in boundary layer flows, particularly, in the near-wall region. Furthermore, the use of hybrid approaches (RANS-LES) to overcome the fine resolution needed in wall-bounded flows is an emerging and promising area [4]; nevertheless, there is still much ground to cover for industrial applications. From the perspective of performing turbulent industrial flow simulations, the workhorse is still the use of the Reynolds Averaged Navier–Stokes equations (RANS), which are obtained by time averaging the full NS equations. In this approach, almost the full power spectra of velocity fluctuations and turbulent kinetic energy are modeled, while capturing the very large turbulent scales of low frequencies. These equations require a model or closure to compute the Reynolds stresses, which arise from the convective terms of the NS equations after applying the time averaging process. The selection of an appropriate turbulence model is crucial to accurately represent important physical aspects of the flow, such as boundary layer separation and shock–boundary layer interaction [5]. Catalano and Amato [5] tested five different turbulence models in 2D and 3D classical aerodynamic applications: Spalart–Allmaras, Myong and Kasagi k-$\epsilon$, Wilcox k-$\omega$, Kok turbulent/non-turbulent (TNT), and Menter SST, where k is the turbulent kinetic energy, $\epsilon$ is the turbulent dissipation and $\omega$ is the specific rate of dissipation. They concluded that the Menter shear stress transport (SST) model exhibited the best balance between the physical capabilities and the numerical robustness for the transonic and high-lift flows explored. Knight et al. [6] performed a review of recent CFD simulations of shockwave turbulent boundary layer interactions. In addition, Knight et al. [6] tested five different configurations: 2D compression corner, 2D shock impingement, 2D expansion-compression corner, 3D single fin and 3D double fin in structured and unstructured grids by considering several two-equation turbulence models. In general, they reported a fair agreement between RANS results to experiments, except for heat transfer predictions of strong separated shock wave turbulent boundary layer interactions where turbulence models still need to make some improvements. In addition, Manzari et al. [7] successfully tested the k-$\omega$ turbulence model in 3D compressible flows on unstructured tetrahedral grids. Furthermore, the use of multi-grid techniques for solution acceleration in 3D turbulent flows on unstructured meshes constitutes a challenging problem [8].

In the present study, the one-equation Spalart–Allmaras model and the two-equation Wilcox k-$\omega$ and Menter SST turbulence models were implemented within an innovative unstructured mesh finite volume flow solver. We computed and validated not only global quantities (e.g., lift, drag and moment coefficient) but also local and boundary layer parameters such as displacement/momentum thickness, pressure coefficient, skin friction coefficient and streamwise velocity profiles. Consequently, this study represents an extensive and thorough analysis of popular turbulence models in large scale systems at experimental Reynolds numbers in compressible flows.

2. Governing Equations

The unsteady compressible Navier–Stokes equations are expressed, over a 3D domain $\Omega \subset {\Re }^3$ with closed surface $\Gamma$, in the integral form

$${\int }_{\Omega }\frac{\partial U_i}{\partial t}d\mathbf{x}+{\int }_{\partial \Omega }{F}_{ij}{n}_{j}d\mathbf{x}={\int }_{\partial \Omega }{G}_{ij}{n}_{j}d\mathbf{x},$$

(1)

where $\overrightarrow{n}=(n_1,n_2,n_3)$ is the unit normal vector to $\Gamma$. In addition, the unknown vector of conservative variables is expressed as

$$U_i=\left[\begin{array}{c}\rho \\ \rho u_i\\ \rho ϵ\end{array}\right]$$

(2)

where $\rho$ is the fluid density, $u_i$ denotes the ith component of the velocity vector and $ϵ$ is the specific total energy and the inviscid and viscous flux vectors are expressed as

$$F_{ij}=\left[\begin{array}{c}\rho u_j\\ \rho u_i{u}_j+p{\delta }_{ij}\\ u_i(\rho ϵ+p)\end{array}\right]G_{ij}=\left[\begin{array}{c}0\\ {\tau }_{ij}\\ u_k{\tau }_{kj}-q_j\end{array}\right]$$

(3)

respectively. The quantity

$${\tau }_{ij}=-\frac{2}{3}\mu \frac{\partial u_k}{\partial x_k}{\delta }_{ij}+\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right),$$

(4)

is the deviatoric stress tensor, where $\mu$ is the dynamic viscosity and ${\delta }_{ij}$ is the Kronecker delta. The quantity $q_j=-k\partial T/\partial x_j$ is the heat flux, where k is the thermal conductivity and T is the absolute temperature. The viscosity varies with temperature according to Sutherlands’s law and the Prandtl number is assumed to be constant and equal to 0.72. In addition, the medium is assumed to be calorically perfect.

3. Solution Procedure

3.1. Hybrid Unstructured Mesh Generation

The computational domain was represented by an unstructured hybrid mesh for viscous problems. The process of mesh generation began with the discretization of the domain boundary into a set of triangular or quadrilateral meshes that satisfy a mesh control function specified by the user [9]. Additionally, the distribution of the mesh parameters [10], such as spacing, stretching and direction of stretching, were described by a background mesh as well as point, line and planar sources. For viscous flows, the boundary layers were generated using the advancing layer method [11] and the rest of the computational domain was filled with tetrahedral elements using a Delaunay incremental Bowyer Watson point insertion [12]. Hybrid unstructured meshes were constructed by merging certain elements of this tetrahedral mesh in the boundary layers.

3.2. Spatial Discretization

The discretization of the governing equations can be performed in a number of different ways. A computationally efficient approach was implemented, consisting on a cell vertex finite volume method. This involved the identification of a dual mesh, with the medial dual being constructed as an assembly of triangular facet, ${\Gamma }_I^K$, where each facet was formed by connecting edge midpoints, element and face centroids in the basic mesh, in such a way that only one node was contained within each dual mesh cell. Hence, the dual mesh cells formed the control volumes for the finite volume process. When hybrid meshes are employed, the method for constructing the median dual has to be modified to ensure that no node lies outside its corresponding control volume. To perform the numerical integration of the fluxes, a set of coefficients was calculated for each edge using the dual mesh segment associated with the edge. The values of the internal edge coefficients, $C_j^{IJ}$, and the boundary edge coefficients, $D_j^{IJ}$, are defined as follows,

$$C_j^{IJ}=\sum _{K\in {\Gamma }_{IJ}}{A}_{{\Gamma }_I^K}{n}_j^{{\Gamma }_I^K},D_j^{IJ}=\sum _{K\in {\Gamma }_{IJ}^B}{A}_{{\Gamma }_I^K}{n}_j^{{\Gamma }_I^K},$$

(5)

where $A_{{\Gamma }_I^K}$ is the area of facet ${\Gamma }_I^K$, $n_j^{{\Gamma }_I^K}$ is the outward unit normal vector of the facet, ${\Gamma }_{IJ}^B$ is the set of dual mesh faces on the computational boundary touching the edge between nodes I and J and $n_j^{{\Gamma }_I^K}$ denotes the normal of the facet in the outward direction of the computational domain. The numerical integration of the fluxes over the dual mesh segment associated with an edge was performed by assuming the flux to be constant, and equal to its approximated value at the midpoint of the edge, i.e., a form of mid point quadrature. The calculation of a surface integral for the inviscid flux over the control volume surface for node I is defined as follows,

$${\int }_{\partial {\Omega }_I}{\mathbf{F}}_j{n}_{j}d\mathbf{x}\approx \sum _{J\in {\Lambda }_I}\frac{C_j^{IJ}}{2}\left({\mathbf{F}}_j^I+{\mathbf{F}}_j^J\right)+\sum _{J\in {\Lambda }_I^B}{D}_j^{IJ}{\mathbf{F}}_j^I,$$

(6)

where ${\Lambda }_I$ denotes the set of nodes connected to node I by an edge and ${\Lambda }_I^B$ denotes the set of nodes connected to node I by an edge on the computational boundary. Thus, the last term is non-zero only in a boundary node. Similar formula can be implemented for the viscous fluxes. The use of edge based data structure has become widely used due to its efficiency in terms of memory and CPU requirements, compared to the traditional element based data structure, especially in three dimensions.

The resulting discretization is basically central difference in character. Therefore, the addition of a stabilizing dissipation is required for practical flow simulations. This was achieved by replacing the physical flux function by a consistent numerical flux function, such as the Jameson–Schmidt–Turkel (JST) flux function [13] or the Harten–Lax–van Leer-Contact (HLLC) solver [14]. Discontinuity capturing may be accomplished by the use of an additional harmonic term in regions of high pressure gradients, identified by using a pressure switch.

3.3. Relaxation Scheme

The FLITE3D flow solver [15] can simulate both unsteady and steady problems. However, in this investigation, a steady flow analysis was performed. Therefore, the relaxation scheme employed consisted of a three-stage Runge–Kutta approach with local time stepping. The viscous and artificial dissipation terms were only computed one time at every time step to reduce the computational requirements of the scheme. The stability limit of this scheme or maximum CFL is 1.8, where CFL is the Courant–Friedrich–Levy number. The local time step values for the momentum and energy equations were computed based on inviscid and viscous term contributions [16]. Similarly, local time steps were implemented for the solution of the corresponding transport equations of the turbulence model approaches (either one- or two-equation models in the present study) with contributions of the inviscid and viscous parts as well as a source term correction. More details can be found in [16].

4. Solver Parallelization

A parallel implementation of the flow solver was achieved by using the single program multiple data concept in conjunction with standard Message Passing Interface (MPI) routines for message passing. In the parallel implementation, at the start of each time step, the interface nodes of the sub-domain with the lower number obtained contributions from the interface edges. These partially updated nodal contributions were then broadcasted to the corresponding interface nodes in the neighboring higher-number domains. A loop over the interior edges was followed by the receipt of the interface nodal contributions and the subsequent updating of all nodal values. The procedure was completed by sending the updated interface nodal values from the higher domain back to the corresponding interface nodes of the lower domain. In addition, the computation and communication processes took place concurrently. The employed approach was a global procedure in which the agglomeration was performed across the parallel domain boundaries. Global agglomeration ensures equivalence with a sequential approach for any number of domain, but increases the communication load in the inner–grid mapping [8,16].

5. Turbulence Modeling in RANS

To obtain the compressible RANS equations, unsteady equations (Equation (1)) were averaged in time to smooth the instantaneous turbulent fluctuations in the flow field, while still allowing the capture of the time-dependency in the large time scales of interest. In many engineering problems, this assumption is valid, but this averaging procedure breaks down if the time scale of the physical phenomena of relevance is similar to that of the turbulence itself. For compressible flows, the density weighted Favre averaging procedure is mostly employed [16]. The Favre averaging procedure applied to Equations (1) generates the extra convective term:

$${\tau }_{ij}^R=-\overline{\rho u_i^{″}{u}_j^{″}},$$

(7)

which is the Favre averaged Reynolds stress tensor. The most straightforward approach is to associate the unknown Reynolds stresses with the computed mean flow quantities by means of a turbulence model or closure. If the Boussinesq hypothesis is applied, this results in a linear relationship to the mean flow strain tensor through the eddy viscosity ${\mu }_t$ [17],

$${\tau }_{ij}^R={\mu }_t\left(\frac{\partial U_i}{\partial x_j}+\frac{\partial U_j}{\partial x_i}-\frac{2}{3}\frac{\partial U_k}{\partial x_k}{\delta }_{ij}\right)-\frac{2}{3}\rho k{\delta }_{ij},$$

(8)

where k is the turbulent kinetic energy. The eddy viscosity depends on the velocity and the length scales of the turbulent eddies, i.e., ${\mu }_t\sim k^{1/2}\ell$, where ℓ is the turbulence length scale. In this study, one- and two-transport-equation models were considered, in which additional partial differential equations were solved to describe the transport of the eddy viscosity. Consequently, nonlocal and history effects on ${\mu }_t$ were taken into account. In the one-equation turbulence model of Spalart–Allmaras [18], a combination of the turbulent scales, through the viscosity-like variable, $\overline{\nu }$, is obtained by solving an empirical transport equation. Two-equation turbulence models are complete, because transport equations are solved for both turbulent scales, i.e., the velocity and the length scale. The original k-$\omega$ model [17] exhibits a freestream dependency of $\omega$, which is generally not present in the k-$ϵ$ model. Menter [19] combined the advantages of both models by means of blending functions, which permits the switching from k-$\omega$, close to a wall, to k-$ϵ$, when approaching the edge of a boundary layer. A further improvement by Menter [19] is a modification to the eddy viscosity, based on the idea of the Johnson–King model, which establishes that the transport of the main turbulent shear stresses is crucial in the simulations of strong Adverse Pressure Gradient (APG) flows. This new approach is called the Menter shear–stress transport model (SST). In particular, the Menter SST turbulence model is well-known for its good performance in boundary layer flows subjected to APG, with eventual separation.

5.1. Turbulent Transport Equations

In this section, a brief description of the governing equations for the three turbulence models is presented. The differential equations are presented in the normalized form. A convenient scaling guarantees a unit order of all variables, which decreases round-off errors of calculations. The reader is referred to Appendix B of Sørensen’s thesis [16] for the scaling and the normalized version of momentum and energy equations employed in the present study; however, the normalization symbol * is dropped here for simplicity. Furthermore, the dimensionless Spalart–Allmaras transport equation [18] is expressed as;

$$\begin{array}{c}\underset{\mathrm{transient}\mathrm{term}}{\underbrace{\frac{1}{S{t}_{\infty }}\frac{\partial \tilde{\nu }}{\partial t}}}+\underset{\mathrm{conv}.\mathrm{term}}{\underbrace{U_j\frac{\partial \tilde{\nu }}{\partial x_j}}}=\underset{\mathrm{production}\mathrm{term}}{\underbrace{c_{b1}[1-f_{v1}]S\tilde{\nu }}}+\underset{\mathrm{diffusion}\mathrm{term}}{\underbrace{\frac{1}{R{e}_{\infty }\sigma }\left\{\frac{\partial }{\partial x_j}\left[\left(\nu +\tilde{\nu }\right)\frac{\partial \tilde{\nu }}{\partial x_j}\right]+c_{b2}\frac{\partial \tilde{\nu }}{\partial x_j}\frac{\partial \tilde{\nu }}{\partial x_j}\right\}}}\hfill \\ -\underset{\mathrm{dissipation}\mathrm{term}}{\underbrace{\frac{1}{R{e}_{\infty }}\left(c_{w1}{f}_w-\frac{c_{b1}}{{\kappa }^2}{f}_{t2}\right){\left(\frac{\tilde{\nu }}{d}\right)}^2}}+\underset{\mathrm{transition}\mathrm{term}}{\underbrace{R{e}_{\infty }{f}_{t1}{(\Delta U)}^2}},\hfill \end{array}$$

(9)

In Equation (9), the normalized quantity $\tilde{\nu }$ is related to the eddy viscosity by

$${\mu }_t=\rho \tilde{\nu }\frac{{\chi }^3}{{\chi }^3+c_{v1}^3},$$

(10)

where $\nu$ is the molecular kinematic viscosity of the fluid, $\chi =\tilde{\nu }/\nu$, and $c_{v1}$ is equal to 7.1.

Additionally, $S{t}_{\infty }=U_{\infty }t/L$ and $R{e}_{\infty }=\rho U_{\infty }L/{\mu }_{\infty }$ are the Strouhal and Reynolds numbers, respectively; $U_j$ represents the time Favre-averaged velocity; and ${\mu }_{\infty }$ is the freestream dynamic viscosity. The rest of the parameters and constants are defined as follows:

$$S=\left|\omega \right|+\frac{\tilde{\nu }}{R{e}_{\infty }{\kappa }^2{d}^2}{f}_{v2},$$

(11)

$$f_{v2}=1-\frac{\chi }{1+\chi f_{v1}},$$

(12)

$$f_{\omega }=g{\left(\frac{1+c_{w3}^6}{g^6+c_{w3}^6}\right)}^{\frac{1}{6}},$$

(13)

$$g=r+c_{w2}r\left(r^5-1\right),$$

(14)

$$r=\min \left(\frac{\tilde{\nu }}{R{e}_{\infty }S{\kappa }^2{d}^2},10\right),$$

(15)

$$f_{t2}=c_{t3}exp\left(-c_{t4}{\chi }^2\right),$$

(16)

$$f_{t1}=c_{t1}{g}_{t}exp\left[-c_{t2}\frac{{\omega }_t^2}{{(\Delta U)}^2}\left(d^2+g_t^2{d}_t^2\right)\right],$$

(17)

$$g_t=\min \left(0.1,\frac{\Delta U}{{\omega }_t\Delta x}\right),$$

(18)

where $\left|\omega \right|$ is the vorticity magnitude, d is the distance from a given point to the nearest wall, $c_{b1}$ = 0.1355, $\sigma$ = 2/3, $c_{b2}$ = 0.622, $\kappa$ = 0.41, $c_{\omega 1}=c_{b1}/{\kappa }^2+(1+c_{b2})/\sigma$, $c_{\omega 2}$ = 0.3, $c_{\omega 3}$ = 2, $c_{t1}$ = 1, $c_{t2}$ = 2, $c_{t3}$ = 1.1, and $c_{t4}$ = 2. Moreover, the variable $\Delta U$ is the difference in velocity between the point and the associated trigger point, $d_t$ indicates the closest distance to such a trigger point or curve, ${\omega }_t$ is the vorticity magnitude at the associated trigger point, and $\Delta x$ is the surface grid spacing at this point. The convective term is conveniently rewritten in conservative form,

$$U_j\frac{\partial \tilde{\nu }}{\partial x_j}=\frac{\partial \left(U_j\tilde{\nu }\right)}{\partial x_j}-\tilde{\nu }\frac{\partial U_j}{\partial x_j},$$

(19)

which introduces an extra source term. Nevertheless, the contribution of this source term is negligible [20]; thus, it is ignored in the present formulation.

Furthermore, the corresponding normalized transport equations in the Wilcox k-$\omega$ model [17] for the turbulent kinetic energy, k, and the specific dissipation rate, $\omega$, in compressible flows read as follows:

$$\underset{\mathrm{transient}\mathrm{term}}{\underbrace{\frac{1}{S{t}_{\infty }}\frac{\partial \left(\rho k\right)}{\partial t}}}+\underset{\mathrm{conv}.\mathrm{term}}{\underbrace{\frac{\partial \left(\rho U_{j}k\right)}{\partial x_j}}}=\underset{\mathrm{production}\mathrm{term}}{\underbrace{{\tau }_{ij}\frac{\partial U_j}{\partial x_j}}}-\underset{\mathrm{dissipation}\mathrm{term}}{\underbrace{{\beta }_k\rho \omega k}}+\underset{\mathrm{diffusion}\mathrm{term}}{\underbrace{\frac{1}{R{e}_{\infty }}\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]}}$$

(20)

$$\underset{\mathrm{transient}\mathrm{term}}{\underbrace{\frac{1}{S{t}_{\infty }}\frac{\partial \left(\rho \omega \right)}{\partial t}}}+\underset{\mathrm{conv}.\mathrm{term}}{\underbrace{\frac{\partial \left(\rho U_j\omega \right)}{\partial x_j}}}=\underset{\mathrm{production}\mathrm{term}}{\underbrace{\alpha \frac{\omega }{k}{\tau }_{ij}\frac{\partial U_j}{\partial x_j}}}-\underset{\mathrm{dissipation}\mathrm{term}}{\underbrace{{\beta }_{\omega }\rho {\omega }^2}}+\underset{\mathrm{diffusion}\mathrm{term}}{\underbrace{\frac{1}{R{e}_{\infty }}\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\omega }}\right)\frac{\partial \omega }{\partial x_j}\right]}},$$

(21)

where ${\beta }_k$ = 0.09, ${\beta }_{\omega }$ = 3/40, ${\sigma }_k$ = 2, ${\sigma }_{\omega }$ = 2 and $\alpha$ = 5/9. The normalized eddy viscosity is defined as:

$${\mu }_t=R{e}_{\infty }\rho k/\omega .$$

(22)

In the Menter SST model for the $\omega$ equation, an extra term is considered in Equation (21), which is called the cross-diffusion term,

$$\underset{\mathrm{cross}-\mathrm{diffusion}\mathrm{term}}{\underbrace{2(1-F_1)\frac{\rho {\sigma }_{\omega 2}}{\omega }\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j}}},$$

(23)

where ${\sigma }_{\omega 2}$ = 0.856. $F_1$ is a blending function defined as,

$$F_1=\tanh \left\{\left[\min \left(\max \left(\frac{\sqrt{k}}{{\beta }_k\omega y},\frac{500\nu }{y^2\omega R{e}_{\infty }}\right),\frac{4\rho {\sigma }_{\omega 2}k}{C{D}_{k\omega }{y}^2}\right)\right]\right\},$$

(24)

and

$$C{D}_{k\omega }=\max \left(\frac{2\rho {\sigma }_{\omega 2}}{\omega }\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j},{10}^{-10}\right).$$

(25)

Thus, the blending function $F_1$ generates values close to one far from the wall (k-$ϵ$ model) and almost zero values near the edge of the boundary layer (k-$\omega$ model). More details can be found in [19]. The main differences of Menter SST equations with respect to the Wilcox k-$\omega$ equations [17] can be summarized as follows: (i) the consideration of a cross-diffusion term in the $\omega$ equation (k-$ϵ$ model), which makes the model insensitive to the freestream boundary condition of $\omega$; (ii) the implementation of a stress limiter for the maximum value of $\omega$ as well as a production limiter to impede the build-up of turbulence in stagnation zones; and (iii) the application of a blending function to compute the corresponding constants of the k-$ϵ$ and k-$\omega$ models.

5.2. Discretization of the Turbulent Transport Equations

The turbulence model equations involve the solution of partial differential equations, which are discretized in a similar manner as the governing equations of the flow. Nevertheless, to avoid instabilities from the convective term, a first order upwind discretization was used plus the addition of a term to introduce adequate dissipation for stabilization purposes:

$$\begin{array}{c}{\int }_{\partial {\Omega }_I}{U}_{j}tu{n}_{j}d\mathbf{x}\approx \sum _{J\in {\Lambda }_I}\frac{C_j^{IJ}}{2}\left(U_j^{I}t{u}_I+U_j^{J}t{u}_J\right)-\left(C_j^{IJ}\frac{U_j^{J}t{u}^J-U_j^{I}t{u}^I}{t{u}_J-t{u}_I}\right)(t{u}_J-t{u}_I)+\hfill \\ \sum _{J\in {\Lambda }_I^B}{D}_j^{IJ}{\mathbf{F}}_j^I,\hfill \end{array}$$

(26)

where tu stands for turbulence unknowns ($\tilde{\nu }$, k and $\omega$). Furthermore, the volume integrals were calculated using the midpoint rule and the gradients appearing in the model were calculated as follows:

$${\int }_{{\Omega }_I}\frac{\partial U_i}{\partial x_j}d\mathbf{x}={\int }_{\partial {\Omega }_I}{U}_i{n}_{j}d\mathbf{x},$$

(27)

Equation (27) in discrete form reads,

$$\frac{\partial U_i}{\partial x_j}{|}^I\approx {\partial }_j^h{U}_i^I\equiv \frac{1}{V_I}\left[\sum _{J\in {\Lambda }_I}\frac{C_j^{IJ}}{2}\left(U_i^I+U_j^J\right)+\sum _{J\in {\Lambda }_I^B}{D}_j^{IJ}{U}_i^I\right],$$

(28)

where $V_I$ is the volume of the control volume. Equations (27) and (28) are expressed only for the velocity but can be applied to any flow parameter. The second-order diffusion term was calculated using the compact stencil form according to Equation (3.38) in [16], where the gradient along the edges are evaluated by means of the compact finite difference scheme.

6. Initial and Boundary Conditions

Generally speaking, all run tests were started from freestream conditions. Due to the similarity of the two-equation models, the steady solution by the Wilcox k-$\omega$ model was used as starting condition for the Menter SST model to save computational resources. The boundary conditions for the flow parameters (i.e., velocity, density and total energy) are discussed in [15]. In this section, focus is given to the corresponding boundary conditions of the turbulence variables. Furthermore, the boundary conditions employed in this investigation were classified as: solid wall, inflow, outflow, symmetry and engine boundaries.

6.1. Solid Wall Boundary

In the Spalart–Allmaras model, the eddy viscosity was set to zero since no velocity fluctuations were present at such locations. Similarly, the turbulent kinetic energy, k, was assigned a zero value at the wall. The specific dissipation rate, $\omega$, did not possess a natural boundary condition at the wall. However, based on the asymptotic solution given by Wilcox [17], the following value for $\omega$ was prescribed according to the implemented normalization:

$${\omega }_o=\frac{6{\nu }_w}{\beta y_w^{2}R{e}_L},$$

(29)

where ${\nu }_w$ is the laminar kinematic viscosity at the wall, $\beta$ is a constant (= 3/40), $y_w$ is the local first off-wall point and $R{e}_L$ is the Reynolds number.

6.2. Inflow and Outflow Boundaries

In the Spalart–Allmaras turbulence model, the independent variable, $\tilde{\nu }$, was set to ten percent of the laminar viscosity at the inflow and outflow. In a similar way and based on recommendations by Spalart and Rumsey [21], the following values for $k_{\infty }/U_{\infty }^2$ and ${\omega }_{\infty }/(U_{\infty }/L)$ were adopted in the two-equation turbulence models at these boundaries: $1\times {10}^{-6}$ and 5, respectively.

6.3. Symmetry Boundary

A symmetry condition is defined as a boundary where the normal velocity component is zero, being the density and total energy normal derivative zero, as well. Furthermore, a symmetry condition can also be compared to an inviscid wall (streamline). Hence, the turbulent eddy viscosity was set to zero at this location in turbulence models. If the configuration of the problem possesses a symmetry plane, this represents a very convenient way to significantly reduce the computational resources required.

6.4. Engine Boundary

Some test cases in the present study considered jet engines, such as the F15 aircraft. The internal parts of the engines were assumed to be outside the computational domain. Therefore, it was necessary to prescribe the corresponding boundary conditions at the engine inlets and outlets. For engine inlets, the mass flow, $m_e$, is defined as follows:

$$m_e=-{\int }_{{\Gamma }_e}\rho u_j{n}_j^{e}d\mathbf{x},$$

(30)

where ${\Gamma }_e$ is the inlet engine surface and $n_j^e$ is the corresponding unit normal. Furthermore, the independent variables of the turbulence model equations were prescribed fixed values at engine inlets and outlets. In general, freestream values for the turbulence variables were set at the engine inflow and outflow.

7. Results and Discussion

The Spalart–Allmaras, Wilcox k-$\omega$ and Menter SST turbulence models were implemented in the FLITE3D flow solver [15]. The results of classical aerodynamic test cases are presented in this section.

7.1. Supersonic Flat Plate

Although the flat plate shows a very simple geometry without a streamwise pressure gradient, it is very appropriate to evaluate the performance of any turbulence model due to the extensive experimental data, theoretical/empirical correlations and high-resolution numerical data available from the literature. In this study, a supersonic flat plate at a freestream Mach number, $M_{\infty }$, equal to 2 was computed. The Reynolds number, based on the streamwise x-coordinate was $1\times {10}^7$ per unit length. The computational domain dimensions shown in Figure 1a were 10.4 m × 0.78 m × 0.78 m along the x-streamwise, y-vertical and z-spanwise directions, respectively. The boundary layer thickness by the end of the flat plate was ${\delta }_{ref}\approx$ 0.06 m, therefore, the computational box in terms of the reference boundary layer thickness was 173${\delta }_{ref}\times 13{\delta }_{ref}\times 13{\delta }_{ref}$. Hence, the computational box was wide enough to eliminate any influence from the lateral faces on the flow statistics. Furthermore, numerical results shown here were taken from the vicinity of the central longitudinal plane, far from lateral surfaces. Additionally, the height of the computational domain was sufficient (∼ 13${\delta }_{ref}$) to allow a natural streamwise developing of the turbulent boundary layer without interference of the top boundary condition. The hybrid mesh consisted of 1.1 million tetrahedral and 0.5 million prismatic elements. The boundary layer consisted of 24 layers with the closest point to the wall located at a distance of approximately $y^{+}\approx 0.7$ in wall units. A close-up of the boundary layer elements can be observed in Figure 1b. A region with a bottom symmetry condition was imposed upstream of the flat plate leading edge. Downstream of the flat plate (where the no-slip condition was assumed for the velocity field), a freestream condition was prescribed. Similarly, at inflow, outflow and top surfaces, freestream values were imposed. An isothermal wall condition Was assumed for the temperature field with a wall to freestream temperature ratio ($T_w/T_{\infty })$ equal to 1.75. This was a 3D case with a symmetry condition assumed in lateral faces. Iso-contours of the turbulent kinetic energy, $k^{*}=k/U_{\infty }^2$, are shown in Figure 1c) from the Menter SST model. The starting and ending points of the flat plate were represented by two cross-sectional $YZ$ cutting planes. The full $X-$length of the flat plate was about 87${\delta }_{ref}$. Figure 1c presents the natural evolution of $k^{*}$ in the direction of the flow, which was mainly responsible for triggering turbulence.

Figure 2 shows the streamwise variation of the skin friction coefficient, $C_f$, as a function of $R{e}_x$. Generally speaking, the corresponding values obtained by Spalart–Allmaras, Wilcox k-$\omega$ and Menter SST turbulence models depicted an excellent agreement with theoretical correlations from White [22] and Schlichting [23], particularly by the end of the flat plate. However, the Menter SST model exhibited a shorter transition and the $C_f$ profile quickly tended to realistic fully-turbulent values downstream from the leading edge. Additionally, Direct Numerical Simulation (DNS) data from Araya and Jansen [24] are also included for a supersonic isothermal flat plate ($T_w/T_{\infty }=$ 2.25) at a freestream Mach number, $M_{\infty }$, of 2.5 and low Reynolds numbers. Again, the Menter SST turbulence model exhibited a good agreement of $C_f$ with high spatial/temporal numerical results. Generally speaking, the mean streamwise velocity along the boundary layer followed the DNS profile of Araya and Jansen [24] and the 1/7 power-law distribution in the outer region (i.e., $y/\delta >$ 0.1), as shown in Figure 3a for $R{e}_x=3.4\times {10}^7$. Nevertheless, all velocity profiles computed from turbulence models depicted some deviations in the inner region (where the 1/7 power law did not work properly) from the DNS profile. The thermal boundary layer plays a crucial role in the transport phenomena of compressible wall-bounded flows; consequently, its accurate understanding and modeling are very important. In other words, the coupled behavior of the velocity and thermal field should be properly represented. In Figure 3b, the computed temperature distributions, $T/T_{\infty }$, as a function of the mean streamwise velocity, $U/U_{\infty }$, are plotted at a streamwise station where $R{e}_x=3.4\times {10}^7$. For comparison, the theoretical Crocco–Busemann relation (see page 502 in White [22]) is also included. The Crocco–Busemann relation assumes a linear variation of the total enthalpy across the boundary layer in zero pressure gradient with a unitary turbulent Prandtl number. Furthermore, present RANS thermal profiles showed good agreement with the quadratic Crocco–Busemann relationship. It is worth mentioning that the Crocco–Busemann relation is a function of the local wall-temperature; therefore, a slightly different theoretical profile was obtained for each turbulence model.

7.2. Transonic RAE 2822 Airfoil

This evaluation test consisted on the simulation of steady turbulent transonic flow around the RAE2822 airfoil, which is the standard test Case 9 [25] for turbulence modeling. The freestream Mach number was $M_{\infty }=0.73$, the Reynolds number was $6.5\times {10}^6$ and the angle of attack was $\alpha =$ 2.79°. This problem was solved using a 3D hybrid and unstructured mesh. The mesh consisted of 4.2 million tetrahedral elements, 1.9 million prismatic elements and 2 thousand pyramid elements. The total number of elements was approximately 6.1 million. The boundary layer consisted of 40 layers with the closest point to the wall located at a distance of $y/c=0.5\times {10}^{-6}$ or approximately $y^{+}\approx$ 0.1–0.4 in wall units.

Table 1. Lift, drag and moment coefficients in RAE2822.

	SA	k-$\mathit{\omega }$	SST	Exp. [25]	Num. [26]	Num. [27]
$C_L$	$0.801$	$0.811$	$0.792$	$0.803$	$0.794-0.824$	$0.824$
$C_{Dtotal}$	$0.0175$	$0.0178$	$0.0182$	$0.0168$	$0.0164-0.0165$	$0.0165$
$C_{Dfriction}$	$0.0054$	$0.0054$	$0.0064$	$N/A$	$0.0054-0.0055$	$0.0055$
$C_{m25\%}$	$-0.1121$	$-0.1144$	$-0.1099$	$-0.099$	$-0.0874$	$N/A$

shows the computed values of lift, drag and moment coefficients for the Spalart–Allmaras, Wilcox k-$\omega$ and Menter SST turbulence models, together with experimental data from the AGARD AR–138 report [25] and other numerical simulations in 2D and structured meshes [26,27]. In Table 1, $N/A$ stands for Not Available. For this simulation case, all turbulence models produced values of $C_L$ within 11 lift counts (1 lift count = 0.001) with respect to experiments by Cook et al. [25]. Similarly, the predicted total drag coefficients, $C_{Dtotal}$, showed good agreement with experiments with a maximum discrepancy of 1.5 drag counts (1 drag count = 0.0001), at most, for the Menter SST model; nevertheless, it generated the most accurate moment coefficient at 25% of the chord, $C_{m25\%}$.

In general, the three models are in good agreement with the experimental results and other simulation values, as shown in Table 1. Furthermore, the pressure ($C_p$) and skin friction ($C_f$) coefficients are depicted in Figure 4, together with the corresponding experimental data from the AGARD AR–138 report [25]. The predicted values of $C_p$ in the upper surface looked very similar in all turbulence models. The weak shock at $x/c\approx 0.05$ was not fully captured by turbulence models with discrepancies in the order of 9% in peaks of $C_p$. Downstream, a quasi-zero pressure gradient zone with almost constant values of $C_p$ and about half-chord in length was observed. The second and stronger shock located at $x/c\approx 0.55$ was accurately captured by all models. Perhaps, the Menter SST slightly separated from the experimental values by the end of the shock (i.e., $x/c\approx 0.57$). The strong shock was characterized by a sharp increase of wall pressure. The presence of a very strong APG induced a small flow separation zone, which is discussed below. Beyond this zone of “pockets” with supersonic flow, the wall pressure kept increasing towards the trailing edge, but at a moderate rate. In addition, the numerical predictions of $C_p$ on the lower surface were almost identical to the experimental values given by all turbulence models. The skin friction coefficient is defined as $C_f=\frac{{\tau }_w}{1/2{\rho }_e{U}_e^2}$, where ${\tau }_w$ is the wall shear stress, and ${\rho }_e$ and $U_e$ are the density and velocity at the edge of the boundary layer, respectively. From the results in Figure 4b, it was inferred that the closest values to experiments were produced by the Menter SST model in the upper side, particularly at the strong shock location (i.e., $x/c\approx 0.55$). Furthermore, the Spalart–Allmaras and Wilcox k-$\omega$ models were observed to perform similarly; however, the Wilcox k-$\omega$ was the only turbulence model that predicted a shock-induced separation due to the presence of a very strong APG, with a small recirculating zone around 0.55 $<x/c<$ 0.60, as shown in Figure 4b. All models slightly underpredicted the experimental value of the skin friction on the lower surface. Figure 5 shows iso-contours of the Mach number for the Wilcox k-$\omega$ model. This transonic regime was characterized by the formation of regions or “pockets" with supersonic flow. Since the flow must return to its freestream conditions by the time of reaching the trailing edge, this caused the formation of a strong shock by $x/c\approx 0.55$ where subsonic flow was observed beyond that point.

During the mesh generation, a good resolution was ensured by clustering the first-off wall point well inside the linear viscous layer (i.e., $y^{+}<4$) to accurately predict the skin friction coefficient. Accordingly, the distance distribution of the first off-wall point in wall units, $\Delta z_w^{+}$, along the RAE2822 airfoil and based on the Menter SST model is shown in Figure 6 together with the airfoil coordinates in the right vertical axis. It can be seen that $\Delta z_w^{+}$ was around 0.4 in the vicinity of the leading edge for the lower surface; however, $\Delta z_w^{+}$ was lower than 0.1 in the upper surface, where the flow faced more complex pressure changes. Consequently, the high spatial resolution of the employed mesh was demonstrated. In Figure 7, profiles of the mean streamwise velocity downstream of the strong shock in the upper surface (i.e., at $x/c=$ 0.65 and 0.75) are depicted. Notice that the mean streamwise velocity was normalized by the local velocity at the edge of the boundary layer, $U_e$. In general, for velocity profiles shown in Figure 7a,b, the Menter SST model produced the most accurate velocity profiles of the three turbulence models tested in this study, when compared to experimental values from Cook et al. [25]. It is interesting to point out that both velocity profiles at $x/c=$ 0.65 and 0.75 were located in a zone of very strong APG or increasing pressure, as shown in Figure 4a. Furthermore, the Shear Stress Transport (SST) formulation by Menter [19] focused on improving the performance of this turbulence model in adverse pressure gradient flows. Thus, the results depicted in Figure 7a,b support the previous statement. In addition, the corresponding displacement thickness, ${\delta }^{*}$, momentum thickness, ${\theta }^{*}$, and, shape factor H ($={\delta }^{*}/{\theta }^{*}$) were computed at these stations. These integral boundary layer parameters (i.e., ${\delta }^{*}$ and ${\theta }^{*}$) are defined as follows for compressible flows,

$${\delta }^{*}={\int }_0^{\delta }\left(1-\frac{\rho U}{{\rho }_e{U}_e}\right)d\mathbf{z},$$

(31)

$${\theta }^{*}={\int }_0^{\delta }\frac{\rho U}{{\rho }_e{U}_e}\left(1-\frac{U}{U_e}\right)d\mathbf{z},$$

(32)

where subscripts e stands for values at the edge of the boundary layer. In

Table 2. Boundary layer parameters at $x/c$ = 0.65 and 0.75 in the RAE2822 airfoil.

Location	Parameter	SA	k-$\mathit{\omega }$	SST	Exp. [25]
$x/c=0.65$	${\delta }^{*}$	$0.004782$	$0.005029$	$0.005304$	$0.004956$
	${\theta }^{*}$	$0.001816$	$0.001794$	$0.002010$	$0.002043$
	H	$2.633$	$2.803$	$2.638$	$2.426$
$x/c=0.75$	${\delta }^{*}$	$0.006571$	$0.006630$	$0.007276$	$0.006884$
	${\theta }^{*}$	$0.003032$	$0.003013$	$0.003204$	$0.002999$
	H	$2.167$	$2.201$	$2.271$	$2.296$

, the computed values of ${\delta }^{*}$, ${\theta }^{*}$ and H are exhibited for the three turbulence models together with experimental data from Cook et al. [25]. Furthermore, it was established that the three models analyzed in this investigation gave a similar level of accuracy on the calculation of the boundary layer parameters; nevertheless, the Menter SST was slightly superior to the other models with an average discrepancy of 5% with respect to experiments. In addition, the knowledge of the shape factor, H, in turbulent flows is important in assessing how strongly the APG is imposed. In fact, the higher is the value of H, the stronger is the APG. Similarly, it is well known that any flow subjected to strong deceleration is prone to separation, where most of the turbulence models find enormous difficulties. In this opportunity, we picked up two stations in the RAE2822 airfoil (i.e., $x/c$ = 0.65 and 0.75) with strong APG or high shape factors. The idea was to evaluate the performance of the turbulence models in extreme conditions. Generally speaking, the Menter SST model showed the best performance in describing the physics of the flow subjected to strong APG.

7.3. ONERA M6 Wing

The third test case considered the simulation of the flow over the ONERA M6 wing [28], where the freestream Mach number was $M_{\infty }=0.84$, the angle of attack was $\alpha ={3.06}^{\circ }$ and the Reynolds number based on the mean geometric chord was $1.2\times {10}^7$. The hybrid mesh employed consisted of 6.3 million tetrahedral elements, 2.5 million prisms and 11.6 thousand pyramids. The surface mesh consisted of 1.1 million triangular elements. The mesh had 35 viscous layers and the first off-wall point was located at $y^{+}\approx 0.4$ in wall units.

The values of lift and drag coefficient for the Spalart–Allmaras, Wilcox k-$\omega$ and Menter SST turbulence models are shown in

Table 3. Lift and drag coefficients in ONERA M6 wing at $\alpha ={3.06}^{\circ }$.

Aerodynamic Coeff.	SA	k-$\mathit{\omega }$	SST	Num. [29]	Num. [30]
$C_L$	$0.260$	$0.262$	$0.253$	$0.270$	$0.253$
$C_{Dtotal}$	$0.0175$	$0.0179$	$0.0189$	$0.0174$	$0.0168$
$C_{Dfriction}$	$0.0048$	$0.0051$	$0.0057$	$0.0050$	N/A

together with numerical results from Le Moigne and Qin [29] and Nielsen and Anderson [30]. A maximum of nine lift counts deviation in the lift coefficient $C_L$ between the three turbulence models was observed. In addition, a maximum of 14 drag counts difference in the drag coefficient $C_D$ was computed, of which nine drag counts were due to the friction contribution to the drag coefficient. Furthermore, the present numerical results of lift and drag coefficients are in good agreement with numerical values by Le Moigne and Qin [29] and Nielsen and Anderson [30], who employed the Baldwin–Lomax and Spalart–Allmaras turbulence model, respectively. The pressure coefficient profiles $C_p$ at different sections across the wing span are depicted in Figure 8. All three models captured reasonably well the corresponding location of both shocks at each wing span section, given by abrupt modifications on $C_p$. The most significant discrepancies were observed at $y/(b/2)=0.8$, where none of the three models were able to appropriately represent the shock within 0.25 <$x/c$< 0.35. However, a similar behavior of the Spalart–Allmaras and Wilcox k-$\omega$ turbulence models at $y/(b/2)=0.8$ was reported by Huang et al. [31], Jakirlić et al. [32] and Nielsen and Anderson [30]. The skin friction coefficient $C_f$ is plotted in Figure 9 at spanwise stations $y/(b/2)$ = 0.9 and 0.95. These figures show the good agreement among all turbulence models about the location of the flow separation point $x/c\sim$ 0.25 and 0.2 for $y/(b/2)$ = 0.9 and 0.95, respectively. The downstream recovery of $C_f$ was different in all turbulence models, with the Wilcox k-$\omega$ and Menter SST turbulence models predicting similar reattachment lengths and Spalart–Allmaras producing the largest reattachment length. The Menter SST model induced higher values for the skin friction downstream of the reattachment point, which was consistent with the largest values of $C_{Dfriction}$ shown in Table 3.

Figure 10 shows iso-surfaces of negative streamwise velocities normalized by the freestream velocity $U_{\infty }$. The recirculating flow bubble in dark blue (identified by a white circle) predicted by the Spalart–Allmaras model was slightly larger than that of predicted by the Menter SST model. Iso-contours of the streamwise velocity, normalized by the freestream velocity $U_{\infty }$, are depicted in Figure 11 in Menter SST at $y/(b/2)=0.9$, where a strong shock wave can be clearly observed around $x/c\sim$ 0.25 of the upper surface. Furthermore, the interaction of this shock with the turbulent boundary layer provoked a recirculating zone represented by the dark blue area.

7.4. Generic F15 Aircraft Configuration

The last example corresponds to a generic F15 aircraft configuration, which represents a very complex and challenging geometry for CFD simulations with concave and convex surface curvatures. The freestream Mach number was $M_{\infty }=0.9$, the Reynolds number, based on the airplane length ($L_x\sim$ 18 m), was $3.58\times {10}^8$ and the angle of attack was five degrees. The hybrid mesh consisted of 12.5 million nodes and 49 million elements (36 million tetrahedra, 12.7 million prisms and 48.5 thousand pyramids). Due to the spanwise symmetry nature of the problem, only half of the model was used for the simulation. In Figure 12a, the employed surface mesh is plotted by mirroring the half-aircraft grid across the symmetry plane. In addition, details of the viscous layers can be seen by means of the cut plane at $x/L_x\approx 0.67$ in Figure 12b. Additionally, the boundary layer was made of 44 layers and the first off-wall point was prescribed at $1\times {10}^{-5}$ m, which ensured that the location of this point was in the linear viscous layer for an accurate skin friction computation. A great deal of thought was given to the design of all computational meshes and dimensions in the present study to ensure grid-independent numerical results. The selected mesh spacing distribution, shown in this section, was the result of grid convergence study that was carried out to obtain grid independent results. The numerical simulations were run on 128 processors.

Table 4. Lift and drag coefficients in F15 aircraft at $\alpha ={5.0}^{\circ }$.

Aerodynamic Coeff.	Spalart–Allmaras	Wilcox k-$\mathit{\omega }$	Menter SST
$C_L$	$0.181$	$0.174$	$0.175$
$C_{Dtotal}$	$0.0397$	$0.0374$	$0.0398$

shows the lift coefficient, $C_L$, and the drag coefficient, $C_D$, for each turbulence model. The numerical results of $C_L$ obtained by the three turbulence models were consistent with a maximum difference of seven lift counts. Furthermore, the largest discrepancy observed in the $C_D$ coefficient was 24 drag counts between the Menter SST and Wilcox k-$\omega$ turbulence models.

In Figure 13a, iso-contours of the pressure coefficient $C_p$ over the entire F15 aircraft configuration are depicted. Strong shock waves were observed on the canopy over the cock pit, leading edge of wings and between the vertical tails. The cut plane in Figure 13b at $x/L_x\approx 0.67$ depicted high values of $C_p$ toward the wing tip. In addition, a top view of iso-contours of $C_p$ are shown in Figure 14. The value of the pressure coefficient at four points (labeled as I–IV in Figure 14) were compared with flight and wind tunnel data from Webb et al. [33], as shown in Figure 15. In general, the present numerical results are in good agreement with the experimental wind tunnel data. Furthermore, it can be noticed that the pressure coefficients given by the Menter SST model showed a better match with experimental data than those calculated by Spalart–Allmaras and Wilcox k-$\omega$ models. In addition, the most significant discrepancies of the numerical results of $C_p$ with experiments were found in the region $0.3<x/L_x<0.5$, just above the turbine intake. The discrepancy can be attributed to the difference between the estimated engine mass flow prescribed and the used engine mass flow during the flight test. This was emphasized by the significant scattering between the experimental data (in flight and wind tunnel), particularly, at Points I and II. In this zone, the different turbulence models predicted different separation regions, which are discussed below. Hence, this resulted in different local flow patterns that induced different pressure coefficients.

The skin friction coefficient $C_f$ at a spanwise location of $y=1m$ or $y/L_y=0.15$, where $L_y=6.52$ m is the semi-span length of the aircraft, is plotted in Figure 16. The study of the drag coefficient due to friction showed that this section had the largest discrepancy between the three turbulence models, as observed in Figure 16a. Figure 16b indicates that the Wilcox k-$\omega$ model predicted two separation bubbles at $x/L_x\approx 0.325$ and 0.36, respectively, while the other two models predicted a single separated region (i.e., at $x/L_x\approx 0.35$) with the Menter SST predicting a much larger separation region than that computed by the Spalart–Allmaras model. Furthermore, the iso-contours of streamwise velocity (normalized by the freestream velocity $U_{\infty }$) by the Menter SST model, as shown in Figure 17, indicated the presence of a lengthy recirculation zone with negative values of the velocity above the turbine intake, given by the shock wave-boundary layer interaction.

Figure 18 and Figure 19 show the pressure coefficients on the upper and lower surfaces of the wing at two different spanwise locations, $y/L_y=0.36$ and 0.59, respectively. In addition, the vertical coordinates, z, of the local profile are represented by the vertical axis on the right. It is shown in Figure 18a that, at $y/L_y=0.36$, the $C_p$ profiles predicted by the three turbulence models showed similar upper peak values in the vicinity of the leading edge, i.e., $C_p\approx -1$. Hence, it can be concluded that the pressure difference between the surface static pressure to the freestream static pressure was equal to the freestream dynamic pressure in absolute value at this location. On the other hand, the Wilcox k-$\omega$ model predicted a stronger second shock wave located at $x/L_x\approx 0.7$. In the lower side, a moderate APG zone was observed in 0.55 $<x/L_x<$ 0.6, as shown in Figure 18b, followed by a nearly ZPG region. Furthermore, Figure 19a shows similar peak values of $C_p$ in the upper surface and close to the leading edge at station $y/L_y=0.59$ by all three turbulence models. However, the Spalart–Allmaras model predicted a much stronger second shock located further downstream of the leading edge (at $x/L_x\approx 0.7$) than those predicted by the Menter SST and Wilcox k-$\omega$ model. All turbulence models exhibited similar $C_p$ distributions in the lower surface, as shown in Figure 19b.

8. Conclusions

An assessment of three popular turbulent models was performed in 3D aerodynamic cases. The FLITE3D flow solver, in conjugation with the Spalart–Allmaras, the Wilcox k-$\omega$ and the Menter SST turbulence models, was applied to the supersonic flat plate, RAE2822 airfoil, ONERA M6 wing, and F15 aircraft.

The Menter SST model exhibited a short developing section in the supesonic flat plate downstream of the edge. The transonic RAE2822 airfoil possessed a strong shock–boundary layer interaction with an induced separation. The ONERA M6 wing exhibited a very complex flow phenomena such as double shocks, strong APG and streamline curvature effects. In addition, the computed global lift and drag coefficients ($C_L$ and $C_D$) over the F15 aircraft by all turbulence models were very consistent; however, local pressure coefficients by the Menter SST model were in better agreement with experimental values in [33].

Generally speaking, numerical results have been very consistent for the three turbulence models. For attached flows or middle separated flows, there was not a clear superiority of the two-equation turbulence models over the one-equation Spalart–Allmaras model. However, the contrary was observed in separated flows with recirculation zones. Particularly, the Menter SST model showed the best compromise between accurately describing the physics of the flow and numerical stability. The Shear Stress Transport (SST) formulation might be the reason for this. Future investigation will focus on unsteady flow simulations. This could be more appropriate in realistically capturing the complex turbulent structures governed by shock-boundary layer interactions and highly-separated flows; in particular, at the transonic regime.