Several Cases for the Validation of Turbulence Models Implementation

This paper presents several test cases that were used to validate the implementation of two turbulence models in the UNS3D code, an in-house code. The two turbulence models used were the Shear Stress Transport model and the Spalart–Allmaras model. These turbulence models were explored using the numerical results generated by three computational fluid dynamics codes: NASA’s FUN3D and CFL3D, and UNS3D. Four cases were considered: a flat plate case, an airfoil near-wake, a backward-facing step, and a turbine cascade known as the Eleventh Standard Configuration. The numerical results were compared among themselves and against experimental data.


Introduction
Computational fluid dynamics (CFD) has become an integral part of the design systems in aerospace and other engineering fields. The increase in computational power in the last two decades allows us to tackle complex, engineering relevant problems. Using today's supercomputers, it has become customary to solve transport phenomena problems with tens or hundreds of millions of grid points and generate results overnight.
As the hardware improves, the models used for predicting transport phenomena should also improve. One standing challenge for flow simulation is the turbulence modeling. As direct numerical simulation is still too computationally expensive for engineering relevant flows, the governing mass, momentum, and energy equations must be closed using a turbulence model.
The widely used κ − turbulence model developed by Jones and Launder [1] cannot properly capture the turbulent boundary layer up to flow separation. The first turbulence model that accurately predicted separated airfoil flows was the Johnson-King model [2]. Being an algebraic model, the Johnson-King model was not easily extensible to threedimensional flow solvers. The κ − ω turbulence model proposed by Wilcox [3] improved the prediction for adverse pressure-gradient flows. In addition, the model has a simple formulation in the viscous sublayer. The κ − ω model, however, depends strongly on the freestream values of the specific turbulent dissipation rate, ω f . This issues was addressed by combining the κ − ω and κ − models, so that the former models the inner region of the boundary layer and the latter models outer part. The resulting model, called the Shear Stress Transport turbulence model [4], was developed to respond to the need for accurate prediction of flows with strong adverse pressure gradients and separation. A turbulence model that is less computationally expensive than the Shear Stress Transport turbulence model is the Spalart-Allmaras model [5]. This is a one-equation model that models an eddy viscosity-like variable,ν.
The purpose of this paper is to provide several cases that can be used to validate the implementation of the Shear Stress Transport and the Spalart-Allmaras turbulence models. The results generated using an in-house CFD code and two NASA codes, using a couple of turbulence models, are compared with experimental data. The next section presents the physical model and provides details on the governing equations. This is followed by a brief summary of the numerical method used to solve the governing equations. The results section presents four validation cases, ranging from a flat plate to a cascade of turbine airfoils.

Physical Model
The governing equations used by the CFD solvers include of the mass, momentum, and energy conservation. These equations are supplemented by the equations used for modeling the turbulence effects. This section briefly presents the integral form of the mass, momentum, and energy conservation equations, as well as the turbulence models used in this work.

Mass, Momentum and Energy Conservation Equations
For an arbitrary finite volume, Ω, bounded by a closed surface, ∂Ω, a small surface, dS, with a unit normal vector,n, pointing outward from the control volume, the mass conservation equation is where ρ is the density and v is the velocity vector. The momentum equation is where the scalar p is the pressure, µ is the dynamic viscosity coefficient, I is the identity tensor, and τ is the viscous stress tensor. The energy conservation equation for a calorically perfect gas is where E is the total energy, H is the total enthalpy, T is temperature, and k is thermal conductivity.

Turbulence Models
Turbulence models are essential for producing accurate results when using Reynolds averaging. They provide the closure needed for the Reynolds-averaged Navier-Stokes (RANS) equations. The following sections describe the two turbulence models used in this work: the Shear Stress Transport (SST) turbulence model and the Spalart-Allmaras (S-A) turbulence model. Both turbulence models employ the Boussinesq eddy-viscosity hypothesis, where the Favre-averaged turbulent shear stresses, τ F , is defined by where µ T is the turbulent eddy viscosity, S ij is the strain rate tensor, κ is the kinetic energy of the turbulent fluctuations, and δ ij is the Kronecker delta.

Shear Stress Transport Turbulence Model
The Shear Stress Transport turbulence [4] model is a two-equation turbulence model that is a combination of a high Reynolds number κ − turbulence model and the κ − ω turbulence model by Wilcox [6]. κ denotes the kinetic energy of the turbulent fluctuations, is the turbulence dissipation per unit mass, and ω is the root mean square fluctuating vorticity (or the rate of dissipation of energy in unit volume and time). The goal of combining the two models is to retain the best attributes of both models while fixing their weaknesses. For example, as the κ − ω model is sensitive to the freestream values of ω, the Shear Stress Transport model switches to the κ − model to resolve this issue. The Shear Stress Transport model is defined by the transport equations of the kinetic energy of the turbulent fluctuations, κ, and of the dissipation ω, f 1 is the blending function between the κ − and the κ − ω turbulence models. This blending allows the Shear Stress Transport turbulence model to switch between the two turbulence models depending on (i) the distance from the wall and (ii) the region of the turbulent boundary layer is being analyzed. Variables β * and σ ω2 are constant coefficients. Variables σ ω , σ κ , C w , and β are non-constant coefficients that utilize the blending function, f 1 [4].

Spalart-Allmaras Turbulence Model
The Spalart-Allmaras turbulence model [5] is a one-equation turbulence model that solves the transport equation of an eddy viscosity-like variable,ν, where P and D are the production and destruction source terms, Constant coefficients include σ, c b1 , c b2 , c w1 , and κ SA . The value d is the distance to the wall. Function f w is known as the wall function, as it is dependent on distance to wall, d. Viscous function f v2 depends heavily on the ratio between the molecular kinematic viscosity, ν L , and the eddy-viscosity-like variableν. The modified vorticity term,S, is a function of the magnitude of vorticity, wall distance, and viscosities [5]. The Spalart-Allmaras turbulence model is a popular turbulence model for turbomachinery flows as it was designed for wall-bounded flows. As the Spalart-Allmaras turbulence model only solves for one equation, it is also inherently less computationally expensive than Shear Stress Transport.
The original turbulence model [5] has an additional trip term, f t1 , used to model transition to turbulence. Allmaras et al. [7] found that the two trip terms-f t1 and f t2 -are unnecessary to simulate fully turbulent flows. The implementation used in this work does not use f t1 [8]. In addition, the implementation is based on the compressible version that includes an update of the modified vorticity term [7].

Numerical Method
The turbulence models tested herein were implemented in three RANS solvers: CFL3D, FUN3D, and UNS3D. The first two flow solvers were developed by NASA. The third flow solver was developed at Texas A&M University. As detailed information is available for the NASA codes, this section focuses on presenting the numerical method used in the UNS3D code. UNS3D stands for Unstructured, Unsteady Three-Dimensional flow solver.

Spatial Discretization
The UNS3D flow solver uses a finite volume method for its spatial discretization. The UNS3D and the NASA's FUN3D codes use a dual-mesh cell-vertex approach, while the NASA's CFL3D code uses a cell-centered approach. The cell-vertex approach has better adaptability with mixed element grids. In addition, as for large grids the number of nodes is smaller than the number of elements, the cell-vertex method is computationally more efficient than the cell-centered method.
3.1.1. Convective Flux UNS3D has several options for calculating the convective fluxes [9]. Herein, we selected an upwind scheme developed by Roe [10] with the additional entropy correction developed by Harten [11], known as the Roe's flux-difference splitting scheme with the Harten entropy fix. The upwind scheme was chosen for its stability and ability to capture shock waves and boundary layers accurately.

Diffusive Flux
The diffusive flux describes flux of φ entering the control volume Ω due to the diffusion of φ. The diffusive flux is proportional to the gradient of φ and is described as ∂Ω kρ(∇φ ·n)dS (10) where k is a diffusivity coefficient, ρ is the density, and ∇φ is the gradient of φ. As the gradients are stored at the nodes, a modified central scheme is used to calculate the edgebased gradients. An ordinary central scheme provides uncoupling between the local terms and the edge-based gradients. The following modification has been added to determine the edge-based gradients: To calculate the diffusive flux, and for second-order spatial discretization, the gradients at the nodes must be calculated. As UNS3D is an unstructured and cell-vertex dual-mesh flow solver, there are specific numerical methods that can be used to calculate the gradients. The following gradient methods are available in UNS3D: Green-Gauss, least-squares, least squares with QR decomposition, and the Weighted Essentially Non-Oscillatory method. The results presented in this paper were generated using the Least Squares with QR decomposition as the primary gradient calculation method.

Second-Order Spatial Discretization
The state variables and their gradients, which are located at the nodes, are required to define a piecewise, linearly-varying approximation for a continuous flow field [12]. To achieve a linear variation, the nodal values across edge "ij" are reconstructed and defined as "left" and "right" state variables.
where Ψ is a limiter function. The limiter functions Ψ range from 0 to 1, and prevent nonphysical oscillations and spurious solutions in the vicinity of large gradients. The solution limiters accomplish this by enforcing a monotonicity preserving scheme. The solution limiters used in this work were those proposed by Venkatakrishnan [13], Carpenter [14], and Dervieux [15].

Temporal Discretization
The terms of the Navier-Stokes equations can be grouped into spatial and temporal derivatives. The spatial derivatives consist of the source terms, and the convective and diffusive flux terms. These spatial derivative terms are grouped together to form the residual R, thus the Navier-Stokes equations can be written as where Q i is the state vector, Ω i is the control volume, and R i is the residual at node i. The UNS3D solver can integrate (13) using either an explicit or implicit method. The results presented in this paper were obtained using only explicit time integration. Explicit time-integration utilizes either a first-order, forward finite difference approximation for the time derivative or a second-order scheme. The second-order time integration uses a four-stage Runge-Kutta method. Implicit residual smoothing is also available in the UNS3D solver. For unsteady simulations, UNS3D utilizes either single-or dual-time stepping.

Boundary Conditions
The implementation of boundary conditions can greatly affect the accuracy and convergence of the solution. The boundary conditions in the UNS3D solver have weak implementations. Weak implementation of boundary conditions refers to the scenario when the boundary conditions are not being directly applied to the nodal values, but instead to the boundary nodes created by the dual mesh. These boundary nodes are used in determining the fluxes from the boundary, and it is in these fluxes that the boundary conditions are satisfied.

Flow Boundary Conditions
For the subsonic inlet, the four conditions imposed from upstream infinity are the two components of the angle of attack, the total pressure, and the total temperature. The condition imposed from the interior of the domain is the outgoing Riemann invariant, u + 2c/(γ − 1). For subsonic outlets, only one condition is imposed from downstream infinity, which is a user-specified static pressure. The conditions enforced from the interior of the domain are the entropy; tangential velocity; and the incoming and outgoing Riemann invariants, u + 2c/(γ − 1) and u − 2c/(γ − 1), respectively.
The wall boundaries are either inviscid and viscous wall boundaries. All wall boundaries enforce the no-penetration condition [16]. This is the only condition for inviscid walls. For viscous flows, the boundary conditions consist of the no-penetration and no-slip conditions [17], p. 80.
Symmetry boundary conditions are used to truncate the computational domain whenever there is a plane of symmetry. Symmetry boundary conditions are also used to simulate a two-dimensional (2D) flow field with a three-dimensional (3D) flow solver, creating what is known as a quasi-3D simulation. There are two conditions that must be satisfied with symmetry boundaries: no-penetration and the zero-gradient of state variables on the planes of symmetry. For this paper, a strong implementation of the zero-gradient condition was enforced, that is, both the boundary state vectors and state vectors on the symmetry boundary had the zero-gradient condition applied [18].

Shear Stress Transport Boundary Conditions
The Shear Stress Transport turbulence model primarily utilizes Neumann boundary conditions. The boundary conditions used are the recommended values from the original reference for the Shear Stress Transport turbulence model [4]. These boundary conditions are listed in Table 1.

Boundary κ ω
Wall The symmetry boundary condition for the Shear Stress Transport turbulence model follows the same structure as the governing equations, where the zero-gradient normal is applied for κ and ω. For the freestream and inlet boundaries, Dirichlet boundaries are enforced for a specific range. The freestream values of ω are obtained from turbulent intensities [4,19]. For the outlet, the values of κ and ω are simply extrapolated. Periodic boundaries are applied to κ and ω in the same manner as to the state vectors in the full-order model.

Spalart-Allmaras Boundary Conditions
The Spalart-Allmaras turbulence model utilizes both the Dirichlet and Neumann boundary conditions for its working variableν. Dirichlet boundary conditions are applied to the inlet/freestream and wall boundaries, while Neumann boundary conditions are applied for symmetry boundaries. Table 2 lists the boundary conditions for the Spalart-Allmaras turbulence model.

Boundary
Dirichlet Neumann The symmetry boundary condition for the Spalart-Allmaras turbulence model follows the same structure as that of the governing equations, where the zero-gradient in the normal direction is applied forν. For the freestream and inlet boundaries, aν/ν L value between 3 and 5 was used to simulate fully turbulent flows. For tripping the laminar-toturbulent transition,ν/ν L must be less than 1. For the outlet, the value ofν is extrapolated. Periodic boundaries are applied toν in the same manner as the state vectors in the mass, momentum, and energy conservation equations.

Results
The cases used for assessing the Spalart-Allmaras turbulence model follow the NASA's turbulence modeling resource website [20]. There are two versions of the UNS3D code: a sequential version, UNS3D-SEQ, and a parallel version, UNS3D-PAR. Both the parallel and sequential versions of UNS3D implemented the Spalart-Allmaras turbulence model and underwent the same tests. The Shear Stress Transport turbulence model was also assessed using the same test cases. The results generated using UNS3D were compared against experimental data and numerical results obtained using FUN3D and CFL3D. The results of the sequential version of UNS3D used a coarser grid than that used by the parallel UNS3D and the NASA codes FUN3D and CFL3D.

Turbulent Flat Plate
The first case investigated was the turbulent flow over a flat plate. The turbulent flat plate case allowed one to validate the turbulence model against Coles' Law of the Wake [21]. The velocity profile and distance from the wall were nondimensionalized to highlight the different regions of the turbulent boundary layer such as the log-region sublayer. Both the sequential and parallel versions of UNS3D were compared against the results from NASA's FUN3D and CFL3D [20]. These results include dimensionless turbulent viscosity µ T /µ ∞ , κ, and ω contours.
The computational domain and boundary conditions are illustrated in Figure 1. Two levels of mesh refinement were used for this case: where the parallel code used the finest grid and the sequential code used the third finest grid [20]. Figure 2 shows the coarse mesh. Table 3 specifies the freestream boundary conditions for the various turbulence models and codes. The freestream turbulence boundary conditions used in NASA's FUN3D SST and CFL3D SST were above the recommended range [4], which is The values used by NASA FUN3D correspond to Therefore, FUN3D used turbulent freestream conditions that were approximately two orders of magnitude higher than the recommended range. UNS3D-PAR SST used values of κ and ω that were in the recommended range: We also generated UNS3D-PAR SST results using the turbulence boundary conditions used by the NASA codes. In this case, UNS3D-PAR SST produced κ and ω contours that were approximately 10 percent smaller in magnitude than the NASA results. The results included in this paper for UNS3D-PAR SST used the turbulent boundary conditions in the recommended range, as stated in Table 3.

Code
Meshν       Figure 6 shows that the main differences between the ω contours generated by UNS3D-PAR and FUN3D using the Shear Stress Transport model were located in the region near the inlet of the domain. FUN3D ω contours started at high values and dissipated as the flow went out of the domain. This difference was to be expected as the turbulent boundary conditions for FUN3D were at least one order of magnitude higher than those of the UNS3D-PAR. Overall, however, the UNS3D-PAR results were in good agreement with those of FUN3D in the region near the wall.  Figure 7 shows the variation of dimensionless velocity u + = u/u * , where u * is the friction velocity, as a function of y + . A good agreement is observed among the results generated by UNS3D-SEQ, UNS3D-PAR, and CFL3D codes, for both the Spalart-Allmaras and Shear Stress Transport models.  Figure 8 shows the skin friction coefficients predicted using the Spalart-Allmaras and Shear Stress Transport models. For the Spalart-Allmaras model, the skin friction coefficient was calculated using UNS3D-SEQ, UNS3D-PAR, and FUN3D codes. The agreement among the results generated by these three codes was excellent. For the Shear Stress Transport model, the skin friction coefficient was calculated using UNS3D-SEQ, UNS3D-PAR, and FUN3D codes. The skin friction coefficient predicted by the UNS3D-SEQ code was slightly smaller than the results predicted by the other two codes. Most likely, this difference is due to the fact that UNS3D-SEQ used a coarser grid than the other two codes, as shown in Table 3.

Two-Dimensional Airfoil Near-Wake
This test case is based off Model A airfoil, a convential airfoil that was experimentally investigated by Nakayama [22]. Model A airfoil is a 10%-thick conventional airfoil with a 61 cm chord. The angle of attack was 0 deg, the wind tunnel speed was 30 m/s, and the freestream turbulence level was 0.02%. The boundaries of the computational domain, shown in Figure 10, were 20 chords away from the airfoil. A detail of the mesh near the airfoil is shown in Figure 11. The input flow parameters used in the numerical simulation are listed in Table 4.   As with the turbulent flat plate case, the freestream values of κ and ω for UNS3D-PAR SST were different from the values used for the FUN3D simulation. Figure 12 shows the contour plots of the dimensionless turbulent viscosity µ T /µ ∞ at the airfoil's trailing edge predicted using the Spalart-Allmaras model. There is good agreement between the results generated by the UNS3D-PAR and FUN3D codes. The turbulent viscosity predicted by the UNS3D-SEQ code, which used the coarse grid, was similar to the two codes; however, the lack of grid refinement at the trailing edge leads to a somewhat thicker wake.  Figures 13 and 14 show the contour plots of the turbulent kinetic energy, κ, and the specific dissipation rate, ω, predicted using the Shear Stress Transport model. The agreement among the results generated by the three codes is good. While there are few differences on the turbulent kinetic energy contour plots, the specific dissipation rates are almost identical.   Figure 16 compares the measured turbulent shear stresses −u v against the predicted ones at the same seven locations in the wake. Both the Spalart-Allmaras and Shear Stress Transport models were used in the UNS3D-PAR code. The two turbulence models predicted the general shape of the shear stresses, although there were some differences between the predicted results and the experimental data.

Two-Dimensional Backward-Facing Step
The backward-facing step is a case proposed by Driver and Seegmiller [23] to test how well the turbulence model simulates separated flows. The backward-facing step case is a challenging problem as turbulence models generally have difficulty modeling separated flows.
The test configuration had a 100 cm long × 15.1 cm wide × 10.1 cm high rectangular inlet duct followed by a backward-facing step of H = 1.27 cm. The geometry of the computational domain of the backward-facing step, shown in Figure 17, is defined as a function of the step height, H. The experiment took place at atmospheric total pressure and temperature. The freestream velocity was 44.2 m/s. The Reynolds number based on step height was Re H = 36,000. The boundary conditions are also shown in Figure 17. A detail of the mesh near the step is shown in Figure 18. Two grids were used for this test case: (1) the second coarsest and (2) the finest grid provided by the NASA Turbulence Modeling [20]. The results generated with the UNS3D-SEQ and UNS3D-PAR codes were compared against those generated with the NASA FUN3D and CFL3D codes which used the finest grid. The sets of results presented herein include velocity profiles, turbulent shear stresses, skin friction coefficients, and reattachment points. Figure 19 compares the measured and predicted velocity profiles at four locations downstream from the step: x/H = 1, 4, 6, and 10. The velocity profiles were calculated using the UNS3D-SEQ, UNS3D-PAR, and CFL3D codes with the Spalart-Allmaras and Shear Stress Transport models. There was a good agreement among the results of the three codes, for both turbulence models. The Shear Stress Transport model matched the experimental data better than the Spalart-Allmaras model near the step, at x/H = 1. Away from the step, at x/H = 4, 6, and 10, the Spalart-Allmaras model captured the velocity profiles better than the Shear Stress Transport model. Figure 20 compares the measured and predicted friction coefficients. There is a relatively good agreement among the results produced by the UNS3D-SEQ, UNS3D-PAR, and CFL3D codes, using both Spalart-Allmaras and Shear Stress Transport models.
The experimental investigation found the reattachment point to be at 6.26 ± 0.01 [23]. The predicted location of the reattachment point varied depending on the turbulence model. The Spalart-Allmaras model predicted an earlier reattachment, while the Shear Stress Transport model predicted a delayed reattachment. The predicted locations of reattachment point are summarized in Table 5.

Eleventh Standard Configuration
The Eleventh Standard Configuration contains two test cases from the experimental data obtained in the annular test cascade at EPF-Lausanne on a two-dimensional turbine cascade: (1) a subsonic, attached flow case, and (2) a transonic, separated flow case [24]. Herein, only the subsonic, attached flow case with stationary airfoils was simulated.
The airfoil chord is 77.8 mm, the pitch is 56.55 mm, and the stagger angle is −40.85 • . The inflow angle is 15.2 • . The input flow conditions are listed in Table 6.  Figure 22 shows the computational mesh for two adjacent airfoils. The flow, however, was calculated over a single passage by using periodic boundary conditions to simulate the rest of the cascade. A multiblock O4H-grid padded with inlet and outlet H-grids was used to discretize the computational domain. The number of grid nodes varied between 39 k on the coarse mesh to 182 k on the fine mesh. Figure 22 shows the coarse grid. In the x − y plane, the fine mesh used 401 × 101 nodes on the O-grid, 62 × 49 nodes on the block prior to the leading edge, 82 × 29 nodes on the block downstream from the trailing edge, 305 × 25 above and below the airfoil, 160 × 97 at inlet, and 200 × 77 at outlet. One cell was used the spanwise direction, as the flow was modeled as two-dimensional and the flow solver was three-dimensional.  Figure 23 shows a comparison of the measured and predicted Mach number over the airfoil. The UNS3D-PAR code was used with the Spalart-Allmaras and Shear Stress Transport turbulence models. The results generated by the coarse and fine meshes were almost identical. Furthermore, the results generated by the two turbulence models were almost identical. On the pressure side, the UNS3D-PAR code results matched the experimental data well. Over the last one-quarter of chord, the predicted flow, however, had a higher velocity than the measured one. On the suction side, the predicted flow had a higher velocity than the measured one over 60% of the airfoil. The UNS3D-PAR code results matched the experimental data better than the numerical results reported in [24].

Conclusions
Four test cases were used to validate the implementation of the Shear Stress Transport and Spalart-Allmaras turbulence models in the UNS3D code. Two versions of the UNS3D code were considered: a sequential code, UNS3D-SEQ, and its parallel version, UNS3D-PAR. The validation of the implementation was done using experimental data and results generated by two NASA codes: FUN3D and CFL3D. The four cases considered were (i) a flat plate case, (ii) an airfoil near-wake, (iii) a backward facing step, and (iv) a turbine cascade known as the eleventh standard configuration. In all cases, the residuals of the solutions were less than 10 −11 . The numerical results generated by the UNS3D codes with the two turbulence models compared well with the experimental data and the results of NASA codes FUN3D and CFL3D. The two turbulence models produced similar results. Compared to the experimental data, the Spalart-Allmaras turbulence model did not predict the turbulent fluctuations and skin friction coefficient as well as the Shear Stress Transport turbulence model. The Spalart-Allmaras model predicted the velocity profiles better than the Shear Stress transport model as the flow moved away from the backward-facing step. Overall, the results generated with the Shear Stress Transport model were closer to the experimental data, the only exception being the prediction of the reattachment point on the backward-facing step, where the Spalart-Allmaras model outperformed the Shear Stress Transport model. The computational time of the Shear Stress Transport model solution exceeded that of the Spalart-Allmaras model by 4% to 38%. In addition, the solution generated by the Shear Stress Transport model was more sensitive to the far-field boundary conditions than that of the Spalart-Allmaras model.