A direct numerical simulation (DNS) of spanwise-rotating turbulent channel flow was conducted for four rotation numbers: , 0.2, 0.5, and 0.9 at a Reynolds number of 8000 based on laminar centerline mean velocity and Prandtl number 0.71. The results obtained from these DNS simulations were utilized to evaluate several turbulence closure models for momentum and heat transfer transport in rotating turbulent channel flow. Four nonlinear eddy viscosity turbulence models were tested and among these, explicit algebraic Reynolds stress models (EARSM) obtained the Reynolds stress distributions in best agreement with DNS data for rotational flows. The modeled pressure–strain functions of EARSM were shown to have strong influence on the Reynolds stress distributions near the wall. Turbulent heat flux distributions obtained from two explicit algebraic heat flux models (EAHFM) consistently displayed increasing disagreement with DNS data with increasing rotation rate.
Turbulent channel flow subject to rotation in the spanwise direction is characterized by reduced turbulence levels near one wall and elevated turbulence levels near the opposite wall; these regions are known as the suction and pressure sides, respectively . Subsequently, the symmetric profiles of mean velocity and Reynolds stress distributions in the nonrotating channel become asymmetric with respect to the channel centerline. The effects of spanwise rotation on momentum transport in turbulent channel flow have been well documented in previous DNS studies by Grundestam et al.  and Wu and Kasagi . Rotation-induced body forces (Coriolis and centrifugal) generate a secondary cross flow, and consequently, a particular type of complex turbulent flow regime is developed with more than one mean flow gradient. The analyses of such complexities on the structure and parameterization of turbulence have relevance to engineering applications, such as gas turbine blade and rotating turbomachinery design, especially with respect to surface heat transfer and skin friction within the internal cooling passages [3,4].
Due to significantly reduced computational costs compared to DNS, Reynolds-averaged Navier–Stokes (RANS) and turbulent heat flux models are highly desirable if they can accurately parameterize turbulent transport. The strong influence of rotational forces on the Reynolds stress and turbulent heat flux transport equations imposes significant challenges for model development especially at high rotation numbers. In this work, we employ DNS to integrate the Navier–Stokes and energy equations and utilize this database to assess four RANS models proposed by (a) Reif et al. , (b) Speziale and Gatski , (c) Girimaji , and (d) Grundestam et al. . We also evaluate two algebraic heat flux models proposed by (e) Younis et al.  and (f) Abe and Suga . In addition, the pressure–strain functions proposed in Refs.  and  are investigated for their influence on the modeled Reynolds stress distributions.
Numerical Method and Case Descriptions
where , ν is the kinematic viscosity, κ is the thermal diffusivity, and the vector is composed of three velocity components in the x (streamwise), y (wall-normal), and z (spanwise) directions, respectively. The rotation number (or Rossby number) is defined as , where Ω is the spanwise angular rotation vector. p is the nondimensional effective pressure () which combines the static pressure (p0) and centrifugal force, and r represents the nondimensional distance away from the axis of rotation. Also, t is the nondimensional time and θ is the nondimensional temperature () with and representing dimensional temperatures on the upper and lower walls, respectively. Equation (2) is uncoupled from Eq. (3), and buoyancy effects are neglected such that present DNS thermal statistics can be compared to the results of similar heat transfer studies involving rotation where negligible buoyancy effects were also assumed. The Prandtl number (Pr) was kept constant at 0.71.
The flow geometry of the present DNS is shown in Fig. 1. Both the velocity and temperature flow fields are assumed to be statistically homogeneous in the streamwise (x) and spanwise (z) directions allowing periodic boundary conditions in those directions. No-slip conditions are imposed on the rigid channel walls (y = 0 and 2). An isothermal boundary condition on the dimensionless temperature θ ( and ) was imposed at the upper and lower walls. Temperature coupling was initiated once the turbulence field had reached a quasi-steady-state solution, and a linear distribution for the temperature field was assumed as an initial condition in the wall-normal direction.
The fractional-step time integration scheme used is a modified semi-implicit Adams–Bashforth/Crank–Nicolson method, which is second-order accurate in time and fourth-order accurate in space. Spatial derivatives were computed using fourth-order central finite differences in order to facilitate parallelization. The wall-normal direction incorporated an exponentially stretched and staggered mesh that clusters points near the two wall boundaries. All the solutions for the elliptic pressure equation and the spatially elliptic momentum equations are obtained using iterative methods available on the Portable, Extensible Toolkit for Scientific Computation Library at the Texas Advanced Computing Center (TACC). To validate the accuracy of the code, numerical solutions were compared with the hydrodynamic linear stability theory for plane channel (Poiseuille) flow. By measuring the growth rate and phase speed of a primary two-dimensional disturbance from the computation and comparing these values to the eigenvalues of the Orr–Sommerfeld equation, the accuracy of the code was assessed. Excellent agreement between the theoretical predictions and DNS results was observed (less than deviation in amplitude), which gives confidence to the viability of the code for solving the unsteady incompressible Navier–Stokes equations. More information on the code scheme and validation may be found in Refs. [11–14].
The full listing of simulations and their corresponding grid resolutions are found in the case descriptions (Table 1). Four simulation cases A–D were conducted for Rossby numbers , 0.2, 0.5, and 0.9, such that models could be tested for a wide range of rotation rates. The Reynolds number based on the laminar centerline velocity was kept constant which resulted in a friction Reynolds number for the no-rotation case (). The low Reynolds number was chosen such that reasonable comparisons could be made between the present DNS and other DNS studies of rotating turbulent channel flow with similar Reynolds numbers [1,15]. The asymmetric velocity distributions due to the rotational effects decreased the value of , which was calculated as an average between the two walls. Simulations were performed at constant mass flux which resulted in a bulk Reynolds number , where is the mean bulk velocity defined as
where denotes a plane-averaged quantity. The bulk Reynolds number of the present DNS simulations was 3800. The domain lengths were selected as , and for all the simulations such that two-point spatial autocorrelations in the streamwise and spanwise directions converged to zero at the largest separations.
Mesh independence was established by designating two high-resolution cases with grid numbers for case A () and for case D (). For both cases, the distributions of mean velocity, Reynolds stresses, and turbulent kinetic energy budgets compared very favorably to those of the original simulations, demonstrating that the selected meshes of the present DNS cases were mesh invariant. The grid spacing is also comparable to other DNS studies of spanwise-rotating turbulent channel flow [1,15].
Superscript “+” refers to nondimensionalization by the friction velocity, or friction temperature, . The global friction velocity for the rotational cases is denoted as , where and are the local friction velocities at the suction and pressure walls, respectively, an equivalent calculation of the global friction temperature was used for the thermal statistics. Unless otherwise specified, all the coordinate directions are nondimensionalized by δ. For all the cases, the governing equations were integrated until both the friction Reynolds number and friction temperature converged availing a sufficiently long time window () to calculate statistics.
In fully resolved turbulent flow, energy spectra should demonstrate an energy density roll-off of several decades with increasing wavenumber and negligible aliasing at high wavenumbers . Sample streamwise power spectra (E) for case A () are shown in Fig. 2. The energy spectra demonstrate sufficient energy drop-off and no aliasing at high wavenumbers in addition to proper characterization of both the inertial and dissipative energy scales.
In Fig. 3, turbulence quantities for simulation case A () are compared with the experimental hot-film data in Ref.  and particle-tracking velocimetry data (PTV) in Ref. ; a PTV distribution was not available for . The present DNS results are shown to compare favourably, especially with respect to the PTV data for and the hot-film data for . In addition, good agreement was also shown for rotating simulation cases B–D with DNS results in Refs.  and . Thermal statistics obtained for simulation cases A–D also compared favorably with the DNS results in Ref. .
Rotational effects on the mean velocity and temperature profiles are shown for simulation cases A–D in Figs. 4 and 5, respectively. In Fig. 4, the mean velocity profile is symmetric about the channel centerline (y = 1) for case A (). With system rotation, the mean velocity distributions become asymmetric as the flow regime is separated into the pressure and suction regions. For rotational cases B–D, a laminarlike (parabolic) profile is observed near the suction wall (y = 0) which is characteristic of the suppressed turbulence in the suction region. As the flow progressively relaminarizes with increasing rotation number, the suction region expands. Near the pressure wall (y = 2), a constant gradient of is shown in the mean velocity profiles at all the rotation rates, which is consistent with previous DNS results in Refs.  and . In Fig. 5, the thickness of the thermal diffusive sublayer, characterized by large near-wall temperature gradients, is broader near the suction wall than the pressure wall for rotating simulation cases B–D. As the rotation number increases, the size of the diffusive layer in the suction region increases and the mean temperature profile shifts toward the pressure wall .
In Fig. 6, maps of fluctuating streamwise vorticity are shown for a y–z crossflow section at for simulation cases A–C. In Fig. 6(a), there are spanwise arrays of alternating high- and low-speed streaks near both channel walls for case A (), which also represent the local intensity of turbulence. In Figs. 6(b) and 6(c), the coherence of the high- and low-speed structures is disrupted near the suction wall for cases B () and C (). With increasing rotation number in the suction region, a significant reduction in the number of turbulence structures is observed along with the expansion of a near-wall region of completely suppressed turbulent activity. In contrast to the suction region, the pressure region continues to demonstrate high levels of turbulent activity, and the coherence of turbulence structures is preserved for all the rotation rates.
In Fig. 7, time-averaged spanwise and wall-normal velocity vectors are shown for a y–z crossflow section at for simulation cases A–C. In Fig. 7(a), no spanwise flow is observed in case A () as expected. In Figs. 7(b) and 7(c), large-scale vortical structures (roll cells) appear from the secondary mean flow advected by the Coriolis force for rotating cases B () and C (). These counter-rotating pairs of roll cells, also known as Taylor–Gortler vortices , are observed to circulate flow throughout the pressure region. With increasing rotation number, the roll cells become more numerous and the size of the circulation region is reduced from increased relaminarization.
To illustrate the effects of rotation on wall shear stress and heat transfer, the dimensionless friction Reynolds () and Nusselt (Nu) numbers for both channel walls are provided for the present simulation cases A–D in Table 2. The introduction of rotation is shown to initially decrease on the suction wall while increasing on the pressure wall. At higher rotation numbers, these trends are shown to significantly weaken or even reverse in case D () for on the pressure wall. These results correspond well with Ref.  which demonstrated that on both walls trended toward convergence at high rotation numbers until the eventual full relaminarization of the flow regime. For the Nusselt number, a dimensionless number generally used to represent surface heat transfer, increasing system rotation is shown to continually decrease Nu on both channel walls. This trend was similarly observed in Ref. .
RANS Turbulence Models
RANS Models: Overview.
The engineering approach for the calculation of turbulent flows generally consists of the use of Reynolds-averaged (time-averaged) Navier–Stokes (RANS) equations which require heuristic closure approximations to model the Reynolds stresses that appear in these equations as a result of the time-averaging process. Accurate parameterization of these quantities in terms of the fundamental variables of the problem, often with the application of an eddy (turbulent) viscosity coefficient, is essential for an accurate prediction of turbulent transport. RANS models are often incorporated into engineering applications but have difficulties especially for the modeling of complex flows with extra mean rates-of-strain. In this work, the current state of turbulence models with respect to spanwise-rotating turbulent channel flow is evaluated through the examination of four commonly used turbulence models in their ability to accurately characterize turbulent transport with rotational effects.
In this work for model testing, nonlinear eddy viscosity models are used due to their advantages over linear eddy viscosity models, such as the k–ε model  for complex turbulent flows. Linear eddy viscosity models invoke a linear constitutive relationship between the Reynolds stresses and mean flow straining field (Boussinesq approximation) and although these models have produced satisfactory predictions in two-dimensional thin shear flows, they perform poorly in complex turbulent flows with more than one mean flow velocity gradient such as the spanwise-rotating turbulent channel flow considered here . In Fig. 8, the Reynolds stress distributions of the k–ε model were compared to the present DNS data and a sample nonlinear eddy viscosity model  for case D (). In comparison to the nonlinear eddy viscosity model, the k–ε model demonstrated significantly poorer accuracy especially with respect to proper characterization of the pressure region (). The k–ε model does produce a very similar distribution with the nonlinear eddy viscosity model for but, in general, a clear disadvantage for the linear eddy viscosity model is observed at high rotation number compared to the nonlinear model.
represents a turbulent timescale.
Most of these closure formulations have been shown to produce significant errors in predicting the mean velocity and temperature distributions for rotating turbulent channel flow (e.g., see, Ref.  for the GWJ model). Hence, in this work, we examine the details of these models to assess the likely sources of such errors in the modeled mean and fluctuating quantities.
RANS Models: Results.
The normal Reynolds stress and primary shear Reynolds stress distributions calculated using the four aforementioned RANS-based models are compared directly to the present DNS results for simulation cases A–D. For the no-rotation case A, the comparisons are shown in Fig. 9. The approximations of the GWJ and GI models are observed to be similar for all the Reynolds stress components. The SG and PRDO models demonstrate the best agreement with DNS results for the normal and shear Reynolds stress distributions, respectively. Also noted in Ref.  is that the SG model is observed to not satisfy realizability as the modeled distribution of contains nonphysical negative amplitudes near both channel walls. Realizability is a set of mathematical and physical principles that need to be satisfied to prevent the turbulence model from generating nonphysical results. The three realizability conditions proposed by Schumann  are
where det is the tensor determinant.
For case B (), the modeled Reynolds stresses and DNS results are shown in Fig. 10 and all the four turbulence models demonstrate proper characterization of the suppressed and enhanced turbulence levels in the suction and pressure regions, respectively. However, the three EARSM-type models (GWJ, SG, and GI) display better agreement with DNS results than the PRDO model, which displayed identical distributions of . The poor performance of the PRDO model is likely the result of the linear relationship between the Reynolds stress and mean strain-rate tensors (Eq. (5)). With respect to the DNS, the SG model showed the best agreement for the normal Reynolds stresses but all the four turbulence models failed to predict the Reynolds shear stress profile in the pressure region. For all the cases, the PRDO and GWJ models produce identical distributions for spanwise Reynolds stress , which was unexpected due to the significant differences in the modeled PRDO and GWJ expressions shown in Eqs. (5) and (10), respectively. The SG model again fails to satisfy realizability as negative amplitudes of are observed near the pressure wall (y = 2) in Fig. 10(b).
The DNS and modeled Reynolds stresses for cases C () and D () are shown in Fig. 11 and 12, respectively. The invariance of the normal Reynolds stress distributions from the PRDO model continues to result in poor agreement with the DNS in relation to EARSM. The GWJ model demonstrates the best correspondence with DNS data for cases C and D. For case C, the modeled SG Reynolds stress distributions develop disparities with the DNS near the pressure wall as the peak modeled amplitude of is significantly overestimated in Fig. 11(d) and an unexpected decrease of amplitudes is observed in Fig. 11(c). For case D, these issues worsen and the modeled SG distributions demonstrate significant deviations from the DNS for all the Reynolds stress components except . For cases C and D, the SG model is shown to satisfy the previously violated realizability condition as the modeled distributions of in Figs. 11(b) and 12(b) remain non-negative throughout the entire channel.
For no-rotation, all the four tested turbulence models produced Reynolds stress distributions in close agreement with the present DNS results. For rotating simulation cases B–D, the EARSM (GWJ, SG, and GI) had significantly better agreement with the DNS than the PRDO model, a different type of nonlinear eddy viscosity model. This observation supports the current trend of turbulence model design for rotational favoring EARSM-type closures . One notable deficiency of EARSM is generally manifested in the near-wall regions where significant deviations of from the DNS distributions are observed. It is therefore instructive to examine the parameterization used for the pressure–strain term as this term is responsible for intercomponent energy transfer.
EARSM: Pressure–Strain Modeling.
The terms on the right-hand side of Eq. (15) represent, respectively, the production (), pressure–strain (), dissipation (), Coriolis (), and diffusion terms ().
In Figs. 13 and 14, the distributions of the pressure–strain closure formula for the SG and GI models are compared to the pressure–strain distributions directly extracted from the present DNS data to evaluate possible sources of error specifically for cases A () and B ().
Figure 13 displays results for case A (), the modeled pressure–strain function values for both models are similar to the DNS results in the interior of the channel. However, near both walls where most of the intercomponent energy transfer takes place, the modeled pressure–strain amplitudes are drastically overestimated for both the SG and GI models. This discrepancy is not unexpected as the SSG function is a single-point closure which captures only the local effects of the pressure–strain correlations, neglecting nonlocal effects such as wavevector information . In Fig. 14, similar comparisons are presented for case B (). It is observed that the modeled pressure–strain distributions have better agreement with the DNS data with increasing rotation number near the suction wall (y = 0) and the only major disparity occurs near the pressure wall with suppressed amplitudes from relaminarization in the suction region.
In Figs. 13 and 14, it is observed that is a loss term for and and are gain terms for their corresponding Reynolds normal stresses. This is expected because the main contribution of the pressure–strain correlation term is to isotropize (return-to-isotropy) the normal Reynolds stress components. In nonrotating turbulent channel flow, the production term contributes only to , so the pressure–strain tensor isotropizes the normal Reynolds stresses through a negative contribution to and positive contributions to and . It is also observed that these pressure–strain models extract too much energy from and transfer excess energy to and , which is directly reflected in the corresponding modeled Reynolds stress distributions in Figs. 9 and 10.
Heat Transfer Models
Heat Transfer Models: Overview.
The terms , and represent heat flux diffusion, production of heat flux through shear forces, production of heat flux through rotational forces, pressure–temperature gradient correlations, and viscous dissipation, respectively.
Heat Transfer Models: Results.
Similar to previous RANS model testing, the YWL and SA modeled turbulent heat flux distributions are directly compared to DNS results for simulation cases A–D. In Figs. 15(a) and 15(b), the DNS and model distributions of the streamwise and wall-normal turbulent heat fluxes, and , respectively, are shown for the no-rotation case A (). Good agreement is displayed with the DNS results for both and . In Figs. 15(c) and 15(d), inconsistencies develop in the YWL and SA modeled heat flux distributions for case B (). Both models show good agreement with the DNS for in the suction region but overestimate the peak amplitudes near the pressure wall (y = 2). In addition, neither model approximates the linear distribution of at the channel center shown in the DNS results.
In Fig. 16, the modeled heat flux distributions from the YWL and SA models are shown and compared with the DNS results for cases C () and D (). In Fig. 16(a), a nonphysical drop-off in distributions from both the SA and YWL models is observed for case C near y = 0.6. In Fig. 16(c), this discrepancy worsens at a higher rotation number in case D (). Increasing departure from DNS data was also observed for increasing rotation number in the modeled distributions of . With an increase in Rossby number, the DNS amplitudes of are shown to become more suppressed in the suction region but the modeled distributions instead display significant amplitude increases of throughout the channel. Through comparison of turbulent heat flux distributions with DNS results for system rotation, both heat transfer models displayed significant deviations that increased with increasing rotation number. As pressure–strain correlations were shown to influence EARSM accuracy, it is instructive to examine the parameterization used for the pressure–temperature gradient correlations in EAHFM.
EAHFM: Pressure–Temperature Gradient Correlation Modeling
n and f denote a unit vector and wall-damping function, respectively.
In Fig. 17, the modeled pressure–temperature gradient distributions from the YWL and SA models are compared to DNS results for selected cases A () and B (). For case A, and are gain terms for the corresponding turbulent heat fluxes, in agreement with Ref. . In Figs. 17(a) and 17(b), both models demonstrate significant overestimations of in the near-wall regions, similar to the tested pressure–strain models. There is significantly better agreement between the modeled and DNS distributions for , and the SA model shows exceptional accuracy.
For case B, continues to be a source term for but changes to become a loss term for in the suction region. This feature is not captured by the YWL and SA models, resulting in significant model disparity from DNS results in the suction region. This major error in the pressure–temperature gradient correlation term is shown to contribute to inaccuracies in the corresponding modeled heat flux distributions in Fig. 15(d). is incorrectly characterized by both models as a gain term for case B, resulting in significant overestimations of modeled in the suction region. The examination of error contributions from and functions to their corresponding models demonstrated the profound effects of pressure correlations on turbulent transport in nonrotating and rotating turbulent channel flow.
For spanwise-rotating turbulent channel flow, the Reynolds stress distributions produced from a linear and nonlinear eddy viscosity model were compared to demonstrate improved accuracy for the nonlinear model over the linear model. In a comparison of four nonlinear eddy viscosity models with DNS data, EARSM were the most compatible with DNS results in modeling the Reynolds stresses for turbulent channel flow subject to spanwise rotation. The Speziale–Gatski (SG) model was shown to be the most compatible for zero and low rotation numbers but displayed significant deviations near the pressure wall at high rotation numbers. The GWJ model showed the best agreement with the DNS data at high rotation numbers. Heat flux models were in good agreement with DNS data for the no-rotation case but with system rotation, the models deviated from the DNS, increasing at higher rotation rates.
The pressure–strain models of two EARSM (Girimaji and SG) were shown to have significant disagreements with the DNS data in the near-wall regions and the pressure–temperature gradient models of two EAHFM (YWL and SA) demonstrated inaccurate characterization of the suction region with system rotation. The errors in the modeled contributions from these terms resulted in degeneration of the predictive capabilities of their respective closure models. These errors contributed to inaccurate Reynolds stress amplitudes in the near-wall regions and an inaccurate modeled distribution shape for wall-normal turbulent heat flux in EARSM and EAHFM, respectively. The present results indicate that correct characterization of pressure fluctuations is a crucial factor in both EARSM and EAHFM design for spanwise-rotating turbulent channel flow.
This work was supported in part by NSF Grant IDR No. 1131802 and by the Air Force Office of Scientific Research under Grant No. FA9550-15-0495. The authors acknowledge the NSF Extreme Science and Engineering Discovery Environment (XSEDE) and Texas Advanced Computing Center (TACC) at the University of Texas at Austin for providing the computational resources that have contributed to the results reported within this paper under Grant No. TG-CTS120002. The authors thank Cesar Galan funded through the University of Colorado Undergraduate Research Opportunity Program for his contributions to this work.
- bij =
- Cij =
diffusion (heat flux transport)
diffusion budget term (heat flux transport)
production from rotational forces (heat flux transport)
- k =
turbulent kinetic energy,
- Nu =
- p =
nondimensional effective pressure,
- Pij =
production from shear forces (heat flux transport)
- p =
quantities on the pressure wall
- p0 =
- Pr =
- qw =
mean heat flux at channel walls,
- r =
nondimensional distance from axis of rotation
- Rec =
centerline Reynolds number,
friction Reynolds number,
- Rob =
bulk Rossby number,
- Roc =
centerline Rossby number,
- s =
quantities on the suction wall
- Sij =
mean strain-rate tensor,
- t =
- ub =
mean bulk velocity
- uc =
mean laminar centerline velocity
- Ui =
mean velocity in xi direction [U, V, W]
fluctuating velocity in xi direction 
turbulent heat flux
- ν =
- Wij =
mean rotation tensor,
- xi =
Cartesian coordinate [x, y, z]
temperature difference between upper/lower walls
- δ =
- δij =
- ε =
turbulent dissipation rate
cyclic permutation tensor
- εij =
viscous dissipation (heat flux transport)
- θ =
- κ =
- τ =
turbulent time scale,
- Πij =
pressure–temperature gradient correlations (heat flux transport)
- Ωi =
angular rotation vector