On the analytical behavior of the $k$--$ω$ turbulence model in buoyant-driven thermal convection

TL;DR

The paper addresses the lack of analytical guidance for incorporating buoyancy in two‑equation RANS closures by deriving an analytical solution for the one‑dimensional Rayleigh–Bénard configuration within the standard $k$-$ ext{ω}$ framework, revealing explicit $Nu$–$Ra$–$Pr$ scalings. It then formulates two algebraic buoyancy corrections that restore observed trends without disrupting the baseline closure, calibrating them against DNS and experimental data. The corrected model is validated across internally heated convection, unstable Couette flow, 2D RB in a square cavity, and side‑heated natural convection, showing substantial improvements in predictive accuracy while preserving numerical robustness. These corrections integrate seamlessly with the existing closure and provide a principled path to improved predictions for buoyancy‑driven turbulence in engineering applications. Limitations relate to gradient‑diffusion closures and dataset coverage; future work could extend the approach to rotation and higher‑order closures.

Abstract

The representation of buoyancy-driven turbulence in Reynolds-averaged Navier--Stokes models remains unresolved, with no widely accepted standard formulation. A key difficulty is the lack of analytical guidance for incorporating buoyant effects, particularly under unstable stratification. This study derives an analytical solution of the standard $k$--$ω$ model for Rayleigh--Bénard convection in an infinite layer, where turbulent kinetic energy is generated solely by buoyancy. The solution provides explicit scaling relations among the Rayleigh ($\mathit{Ra}$), Prandtl ($\mathit{Pr}$), and Nusselt ($\mathit{Nu}$) numbers that capture the simulation trends: $\mathit{Nu} \sim \mathit{Ra}^{1/3}\mathit{Pr}^{1/3}$ for $\mathit{Pr} \ll 1$ and $\mathit{Nu} \sim \mathit{Ra}^{1/3}\mathit{Pr}^{-0.415}$ for $\mathit{Pr} \gg 1$. This framework quantifies the discrepancies in the conventional buoyancy treatment and clarifies their origin. Informed by this analysis, the buoyancy-related modelling terms are reformulated to recover the measured $\mathit{Nu}$--$\mathit{Ra}$--$\mathit{Pr}$ trends. Only two dimensionless algebraic functions are introduced, which vanish in the absence of buoyancy, ensuring full compatibility with the standard closure. The corrected model is validated across a range of buoyancy-driven flows, including two-dimensional Rayleigh--Bénard convection, internally heated convection in two configurations, unstably stratified Couette flow, and vertically heated natural convection with varying aspect ratios. Across all cases, it provides highly accurate predictions.

Figures (17)

• Figure 1: (a) Rayleigh--Bénard convection setup, (b) mean temperature distribution, and (c) turbulent heat flux.
• Figure 2: One-dimensional Rayleigh--Bénard convection predictions of the $k$--$\omega$ model with $C_{\omega b}=1$. (a) Nusselt number as a function of the Prandtl number. Circles, squares, and triangles with dashed lines represent simulation results for $\mathit{Ra} = 10^{6}$, $10^{9}$, and $10^{12}$, respectively, while solid lines denote the analytical prediction. (b) Rayleigh-number dependence at $\mathit{Pr} = 0.7$. The simulation results (circles) and analytical prediction (solid line) are compared with the experimental correlation of niemela2000turbulent (dash-dotted line).
• Figure 3: One-dimensional Rayleigh--Bénard convection predictions of the $k$--$\omega$ model with $C_{\omega b}=1$. (a) Kinetic and (b) thermal boundary-layer thicknesses, (c) the temperature difference across the kinetic boundary layer, and (d) the maximum turbulent viscosity at the domain centre as functions of the Prandtl number. Circles, squares, and triangles with dashed lines denote simulation results for $\mathit{Ra} = 10^{6}$, $10^{9}$, and $10^{12}$, respectively. Solid lines represent the analytical relation.
• Figure 4: One-dimensional Rayleigh--Bénard convection predictions of the $k$--$\omega$ model with $C_{\omega b}=1$. The distribution of $\omega$ at $\mathit{Ra} = 10^{12}$ is shown for (a) $\mathit{Pr} = 10^{-2}$ and (b) $\mathit{Pr} = 10^{2}$. Dashed lines denote simulation results, and solid lines denote the analytical relation. Vertical dotted lines mark the analytical kinetic boundary-layer thicknesses.
• Figure 5: (a) Dependence of $\mathit{Nu}$ on $\mathit{Pr}$ using the applied corrections (Eqs. \ref{['eq: Cwb correction']} and \ref{['eq: Pr Correction integral target']}) in one-dimensional Rayleigh--Bénard convection. Markers denote reference data, and the overlapping solid lines indicate the corrected analytical relations. (b) Reference datasets used for the correction, shown in the $\mathit{Ra}$--$\mathit{Pr}$ space: circles, crosses, diamonds, and pluses denote the datasets of pandey2022convective, xia2002heat, stevens2011prandtl, and li2021effects, respectively.