Assessment of Turbulence Models in a Hypersonic Cold-Wall Turbulent Boundary Layer

In this study, the ability of standard oneor two-equation turbulence models to predict mean and turbulence profiles, the Reynolds stress, and the turbulent heat flux in hypersonic cold-wall boundary-layer applications is investigated. The turbulence models under investigation include the one-equation model of Spalart–Allmaras, the baseline k-ω model by Menter, as well as the shear-stress transport k-ω model by Menter. Reynolds-Averaged Navier-Stokes (RANS) simulations with the different turbulence models are conducted for a flat-plate, zero-pressure-gradient turbulent boundary layer with a nominal free-stream Mach number of 8 and wall-to-recovery temperature ratio of 0.48, and the RANS results are compared with those of direct numerical simulations (DNS) under similar conditions. The study shows that the selected eddy-viscosity turbulence models, in combination with a constant Prandtl number model for turbulent heat flux, give good predictions of the skin friction, wall heat flux, and boundary-layer mean profiles. The Boussinesq assumption leads to essentially correct predictions of the Reynolds shear stress, but gives wrong predictions of the Reynolds normal stresses. The constant Prandtl number model gives an adequate prediction of the normal turbulent heat flux, while it fails to predict transverse turbulent heat fluxes. The discrepancy in model predictions among the three eddy-viscosity models under investigation is small.


Introduction
Accurate modeling of cold-wall hypersonic turbulent boundary layers (TBLs) is critically important to the prediction of the surface heat flux, and hence, to the design of thermal protection systems for hypersonic vehicles. Yet, the existing literature on hypersonic TBLs is rather limited, whether in regard to measurements, modeling, or numerical simulations [1][2][3][4][5][6][7][8][9][10]. Reynolds-Averaged Navier-Stokes (RANS) turbulence models have been widely used for simulating hypersonic turbulent boundary layers, without any specific considerations regarding the high-Mach-number effects. Standard one-or two-equation turbulence models, such as the Spalart-Allmaras (SA) Model and the Menter's Shear Stress Transport (SST) Model, have been developed largely based on studies of subsonic or moderately supersonic flows with adiabatic walls. Thus, careful assessment of model performance is necessary before such models can be extended to applications in the hypersonic, cold-wall regime. Despite many previous efforts to assess existing models for hypersonic applications (see, for example, the review by Roy and Blottner [2] and references therein), continued research to develop better physics-based compressible turbulence modeling is clearly needed for flows in the hypersonic regime, starting with attached boundary layers.
One of the primary objectives of this paper is to assess some commonly used eddy-viscosity turbulence models (Spalart-Allmaras, Menter's k-ω Baseline, and SST models) under hypersonic, cold-wall conditions, using a newly developed direct numerical simulation (DNS) database of supersonic and hypersonic turbulent boundary layers [10]. The flow configuration we focus on is a spatially-developing, zero-pressure-gradient, flat-plate turbulent boundary layer at, nominally, Mach 8, with a wall-to-recovery temperature ratio of T w /T r = 0.48. The current study contains detailed comparisons of mean and turbulence profiles against the DNS that are complementary to the few model evaluation efforts in the literature under high-Mach-number, cold-wall conditions that were limited to comparing wall quantities (skin friction, Stanton number, and wall pressure) with experimental data and empirical correlations [4,11].
Another objective of the current study is to use high-fidelity DNS to assess compressibility corrections for hypersonic, cold-wall turbulent boundary layers. Rumsey [11] provided a recent review of the most common classes of compressibility correction for the k-ω form of two-equation models. He also assessed the influence of compressibility corrections on turbulent skin friction in hypersonic boundary layers. The current study will extend such an assessment to mean velocity and temperature fields.

DNS Simulation of Hypersonic Turbulent Boundary Layers
To provide benchmark data for testing RANS models, the DNS of hypersonic turbulent boundary layers was conducted over a flat plate. The DNS simulation with flat plate is referred to as Case "DNS-FlatPlate". Targeted flow conditions in the free stream for Case DNS-FlatPlate are summarized in Table 1, including the freestream Mach number M ∞ , velocity U ∞ , density ρ ∞ , and temperature T ∞ . The wall is assumed to be isothermal, with a wall temperature of T w = 298.0 K. The corresponding wall-to-recovery temperature ratio is T w /T r = 0.48, with the recovery temperature estimated as based on a recovery factor of r = 0.89, and γ is the specific heat ratio. Case DNS-FlatPlate corresponds to the DNS reported (Case M8Tw048) in a previous paper by our group [10], in which the underlying methodology and validation are presented in detail. The DNS case falls within the perfect gas regime. The working fluid is nitrogen, and its viscosity was calculated by using the Keyes law [12]. A constant molecular Prandtl number of 0.71 was used for the DNS case. To perform DNS of turbulent boundary layers, the full three-dimensional compressible Navier-Stokes equations in conservation form are solved numerically in curvilinear Cartesian coordinates. The boundary layer is simulated in a rectangular box over a flat plate with spanwise periodic boundary conditions and a modified rescaling/recycling method for inflow turbulence generation [13]. The inviscid fluxes of the governing equations are computed using a seventh-order weighted essentially non-oscillatory (WENO) scheme. Compared with the original finite-difference WENO introduced by Jiang and Shu [14], the present scheme is optimized by means of limiters [15,16] to reduce the numerical dissipation. The viscous fluxes are discretized using a fourth-order central difference scheme, and time integration is performed using a third-order low-storage Runge-Kutta scheme [17]. The simulation involves a single domain with a long streamwise box, as illustrated in Figure 1. A detailed description of the problem formulation, the numerical scheme, and the initial and boundary conditions can be found in References [13,[18][19][20][21][22]. The validity of the numerical methods and procedures have been established in multiple previous publications [10,18,21,22].

RANS Simulation of Hypersonic Turbulent Boundary Layers
RANS simulations of high-speed turbulent boundary layers were conducted, with the flow conditions and thermodynamic equation of state matching that of the DNS case listed in Table 1. In RANS, ANSYS Fluent (Version 18.0, ANSYS, Inc., Canonsburg, PA, USA) [23] is used to solve the compressible Favre-averaged Navier-Stokes equations [24], which can be written as: where In the one-or two-equation turbulence models under investigation, the Reynolds stress term τ ij = −ρu i u j was modeled by the Boussinesq approximation: where S ij is the mean velocity strain, and µ t is the eddy viscosity obtained by a turbulence model to be discussed later, and k is the turbulent kinetic energy. The turbulent heat flux term was modeled based on the Reynolds analogy as: where the turbulent Prandtl number Pr t is assumed to be a constant of 0.85. The terms associated with molecular diffusion σ ij u i and turbulent transport − 1 2 ρu i u i u j in the Farve-averaged energy equation are neglected in the current study [25].
In the current study, two-dimensional planar RANS were conducted with a few representative turbulence models, selected among those in widespread use in the aerospace industry. The models chosen are the one-equation model of Spalart-Allmaras [26], the baseline k-ω model by Menter [27], and the shear-stress transport k-ω model by Menter [27]. Unless otherwise mentioned, the "standard" form" (i.e., the original published form) of the model equations was used, with the detailed information on the model formula and coefficients listed in Reference [28][29][30]. Similar to the DNS, ideal gas relations were used in the RANS simulations with nitrogen as the working fluid, and the Keyes law [12] was used for dynamic viscosity.
As far as the compressibility correction is concerned, a complete investigation of all forms of compressibility correction in the literature was not among the goals of this work. We chose to focus on the Wilcox form of compressibility correction [31], which was specifically developed for hypersonic applications. In Wilcox compressibility correction, correction is achieved through the modification of the coefficient in the k-ω destruction term [25,31]. For instance, the modified coefficient of the destruction term in the ω equation is: where ζ * = 1.5 and β is the original coefficient. The compressibility function F(M t ) is obtained as: where M t = √ 2k/a is the turbulent Mach number, M t0 = 0.25, and a = √ γRT is the speed of sound. Figure 2 shows a schematic of the computational domain for Case RANS-FlatPlate along with the boundary conditions. The streamwise (x) and wall-normal (y) domain sizes are L x /δ r × L y /δ r ≈ 93.5 × 19.8, where the reference length δ r = 35.3 mm is the boundary-layer thickness at the center of the domain. Here, the boundary-layer thickness δ is defined as the wall-normal height at which the streamwise velocity reaches 99% of the freestream velocity. Throughout the paper, the velocity components in the streamwise, wall-normal, and spanwise (z) directions are u, v, and w, respectively. A mesh of 551 × 293 grid is used in streamwise and wall-normal directions, respectively. A uniform grid distribution is used in the streamwise direction with a resolution of ∆x/δ r ≈ 0.17. In the wall-normal direction, a geometric distribution with a stretch ratio of less than 1.05 is used to cluster meshes near the wall. The wall-normal grid resolution is ∆y + ≈ 0.2 at the wall and ∆y + ≈ 12 near the boundary layer edge. The grid resolutions are normalized by the viscous length z τ at x/δ r = 66.0, where the turbulence statistics are reported, and the value of z τ is listed in Table 2. To monitor grid convergence, a grid study was conducted using the k-ω SST model on the baseline grid (551 × 293), along with two successively coarser grids for which every other grid point was removed in each coordinate direction (276 × 74 and 551 × 147) and one refined grid (1101 × 293) with two times higher grid resolution in the streamwise direction. Results are shown in Figure 3, where C f is the wall skin friction coefficient with τ w the wall shear stress, and C h is the wall heat flux coefficient with q w the wall heat flux and C p the heat capacity at constant pressure. Over most of the plate, the four grids yielded very close results. The difference in C f between the baseline grid and the finest grid is less than 0.4% at Re θ = 9000.

Results
In this section, the RANS of hypersonic turbulent boundary layers were conducted with different turbulence models. The performance of each model is assessed by comparing RANS against DNS. Figure 4 shows comparisons of the wall skin friction coefficient C f and wall heat transfer coefficient C h among DNS, RANS, and the empirical correlations of van Driest [32] and Spalding-Chi [33]. In general, the RANS cases give good predictions of surface skin friction and heat flux, especially when the k-ω SST model is used. For all RANS cases, the predictions of C f and C h lie within 5% of those of the DNS. Consistent with previous findings [11,34], the van Driest correlation gives significantly better predictions of skin friction and wall heat flux coefficients than the Spalding-Chi correlation.
Additionally, boundary-layer integral and wall parameters and profiles predicted by DNS and RANS were compared with a common friction Reynolds number of Re τ ≈ 500. Table 2 summarizes the boundary-layer parameters at the selected location for both cases, including the momentum thickness θ, shape factor H = δ * /θ (where δ * is the displacement thickness), boundary layer thickness δ, friction velocity u τ = τ w /ρ w , viscous length z τ = µ w /ρ w u τ , wall skin friction coefficient C f , Reynolds analogy factor RA = 2C h /C f (where C h is the wall heat transfer coefficient), dimensionless wall heat transfer rate B q = q w /(ρ w C p u τ T w ), and different definitions of the Reynolds number, namely Re θ = ρ ∞ U ∞ θ/µ ∞ , Re τ = ρ w u τ δ/µ w , and Re δ2 = ρ ∞ U ∞ θ/µ w . Figure 5 further shows a comparison in the mean boundary-layer profiles between RANS and DNS. The selected eddy-viscosity turbulence models, in combination with a constant Prandtl number model for turbulent heat flux, give good predictions of boundary-layer parameters and mean profiles. The discrepancy in model predictions among the different models is small. In terms of turbulence quantities, Figure 6 compares the Reynolds stresses between RANS and DNS cases. In RANS, the Reynolds stresses were derived using turbulent viscosity and the mean velocity gradient, according to the Boussinesq assumption (Equation (4)). Good comparison was achieved for the Reynolds shear stress u v . However, the Reynolds streamwise stress u u was significantly underpredicted by RANS. Figure 7 shows that the turbulent normal heat flux ρv T , calculated according to the Reynolds analogy assumption (Equation (5)), compares well with the DNS. The turbulent transverse heat flux (ρu T ), however, is not properly predicted by RANS.
Lastly, the Wilcox form of compressibility correction [31] is investigated. Figure 8 plots the predictions of boundary-layer mean profiles by Menter's k-ω SST model with and without compressibility correction. The comparison of RANS with the DNS shows that the Wilcox form of compressibility correction causes only a small change in the mean boundary-layer profiles, and no improvement is shown over the standard k-ω SST model for the current flow. A similar trend is found when the compressibility correction is used with the k-ω BSL model.

Conclusions
In this article, the performance of three commonly used eddy-viscosity models (Spalart-Allmaras, k-ω SST, k-ω BSL) was assessed under hypersonic cold-wall conditions. The model evaluation effort focuses on comparing mean and turbulence profiles with DNS for a spatially-developing, zero-pressure-gradient, flat-plate turbulent boundary layer at nominally Mach 8 with T w /T r = 0.48.
The study shows that the selected eddy-viscosity models provide good predictions for turbulent skin friction, wall heat flux, and the mean mass flux and temperature profiles. The use of Boussinesq and Reynolds analogy assumptions also provides acceptable predictions to the Reynolds shear stress, u v , and the turbulent normal heat flux, ρv T . Such assumptions, however, lead to wrong predictions of Reynolds normal stresses ( u i u i ) and the turbulent transverse heat flux (ρu T ). The Wilcox form of compressibility correction was found to cause only a small change in the mean boundary-layer profiles, and no improvement was shown over the standard k-ω model without correction.