k−ω
1. Introduction
With hypersonic vehicles emerging as a new frontier in aerodynamic research and engineering, accurate prediction of aerodynamic heating is critical. Recent developments in Reynolds-averaged models allow them to reasonably predict the mean velocity of a cooled hypersonic turbulent boundary layer (HTBL), but they can still overpredict its off-wall temperature maximum (Danis & Durbin Reference Danis and Durbin2022; Durbin & Yin Reference Durbin and Yin2026; Yin Reference Yin2026). The discrepancy is associated with strong wall-normal variations of density and viscosity and with the balance among molecular diffusion, turbulent heat flux and turbulent kinetic-energy transport in the mean energy equation. The present study addresses the closure of the latter two turbulent contributions.
Regarding compressible wall-bounded turbulence, the hypothesis of Morkovin (Reference Morkovin1962) states that turbulence retains an essentially incompressible structure when mean-flow and thermodynamic effects are properly taken into account. Despite its breakdown in the near-wall region, the hypothesis underpins semi-local transformations and the development of compressible turbulence models.
Trettel & Larsson (Reference Trettel and Larsson2016) and Patel et al. (Reference Patel, Boersma and Pecnik2016) established the modern semi-local-transformation framework by showing that wall-normal density and viscosity variations modify the effective near-wall scales and, in cooled-wall cases, can invalidate simpler density-weighted transformations. More recently, accurate compressibility transformations for wall turbulence have been proposed based on eddy-viscosity equivalence or inner-/outer-layer scaling (Griffin et al. Reference Griffin, Fu and Moin2021; Danis & Durbin Reference Danis and Durbin2024), which facilitate model development (Griffin et al. Reference Griffin, Fu and Moin2023; Hu et al. Reference Hu, Huang, Durbin and Yang2025).
Research on turbulence modelling of compressible wall-bounded turbulence has progressed much further than research on turbulent heat transfer. Recently, the viscous-layer and near-wall corrections developed by Danis & Durbin (Reference Danis and Durbin2022) and Hu et al. (Reference Hu, Huang, Durbin and Yang2025) have yielded improved mean-velocity predictions, but the accompanying improvements in temperature prediction remain limited. Yin (Reference Yin2026) proposed an alternative, simpler formulation to Danis & Durbin’s (Reference Danis and Durbin2022) work, although the improvement in temperature prediction is smaller. Durbin & Yin (Reference Durbin and Yin2026) later reformulated the
$\omega$
equation to correct the near-wall behaviour of
$\omega$
, where
$\omega$
is eddy frequency, which improves the prediction of the
$k-\omega$
model in HTBLs. All the above models tend to overpredict the temperature peaks in HTBLs under wall cooling, consistent with Aiken et al.’s (Reference Aiken, Boyd, Duan and Huang2022) report on the Spalart–Allmaras and shear stress transport (SST) models. This persistent discrepancy suggests that, although existing compressibility corrections are effective for momentum fields, significant deficiencies may remain in the modelling of the Reynolds-averaged Navier-Stokes (RANS) energy equation.
Despite plausible progress in compressible turbulence modelling, its counterpart in heat-transfer modelling has made limited progress. Existing research has largely focused on turbulent-Prandtl-number (
$\textit{Pr}_t$
) modelling for incompressible heat-transfer problems, including aerodynamic flows (Cebeci Reference Cebeci1973; Wassel & Catton Reference Wassel and Catton1973; Kays & Crawford Reference Kays and Crawford1993) and low-Prandtl-number fluids (Cheng & il Tak Reference Cheng and il Tak2006; Taler Reference Taler2018).
Although a constant
$\textit{Pr}_t \approx 0.9$
is the most common choice in low-speed aerodynamic flows, it is no longer sufficient at high Mach numbers, especially in the near-wall area with wall cooling (Chen et al. Reference Chen, Gan and Fu2024, Reference Chen, Gan and Fu2025). More recent efforts have extended that line from traditional experimental fitted algebraic closures to scaling- or Reynolds-analogy-based algebraic models for high-speed wall turbulence, guided by semi-local scaling of the heat-transfer process (Cheng & Fu Reference Cheng and Fu2024). These include semi-local-scaling-based eddy conductivity and temperature-scaling constructions (Chen et al. Reference Chen, Lv, Xu, Shi and Yang2022; Huang et al. Reference Huang, Coleman, Spalart and Yang2023), as well as data-driven models for high-speed Reynolds-stress discrepancy learning (Huang et al. Reference Huang, Duan, Wang, Sun and Xiao2017; Wang et al. Reference Wang, Huang, Duan and Xiao2019) and heat-transfer prediction (Parish et al. Reference Parish, Ching, Jordan, Nicholson, Miller, Beresh, Barone, Gupta and Duraisamy2025; Zafar et al. Reference Zafar, Zhou, Roy, Stelter and Xiao2026). Gupta & Duraisamy (Reference Gupta and Duraisamy2026) recently examined machine-learning-augmented turbulence modelling for wall-heat-transfer prediction in hypersonic shock–boundary-layer interactions, highlighting the importance of augmentation form, feature design and generalisability. However, many of these recent closures are either semi-local-scaling-based or ad hoc in nature, which can make implementation and extension to broader flow conditions challenging.
At the same time, recent literature has also made it clear that the heat-transfer problem cannot be sufficiently handled by tuning
$\textit{Pr}_t$
alone. Hasan et al. (Reference Hasan, Elias, Menter and Pecnik2025) and Mo et al. (Reference Mo, Gao, Li and Jiang2026) showed that temperature prediction critically depends on the treatment of the modelled energy equation and that correcting
$\textit{Pr}_t$
alone has limited leverage, while turbulent kinetic energy-diffusion terms become essential in the heat-transfer process of HTBLs with wall cooling. In a related wall-model large-eddy simulation setting, Kumar & Larsson (Reference Kumar and Larsson2025) likewise improved heat-flux prediction by introducing a modelled turbulent kinetic-energy (TKE) diffusion contribution together with a modified thermal-eddy-diffusivity treatment. These studies suggest that the primary closure deficiency should be sought not only in the turbulent heat-flux model, but also in the modelled TKE transport contribution appearing in the energy equation.
Accordingly, the objective of this study is to develop a local closure for the two turbulent contributions identified above within the RANS
$k-\omega _0$
framework (Durbin & Yin Reference Durbin and Yin2026), with turbulent kinetic-energy transport as the central component. Analysis of the HTBL direct numerical simulation (DNS) database of Fan et al. (Reference Fan, Li and Pirozzoli2022) shows that, near the off-wall temperature maximum, the turbulent heat flux vanishes and the normally secondary kinetic-energy transport becomes locally important to the wall-normal energy balance. This contribution, together with several other terms in the Reynolds-averaged energy equation, is often omitted in low-Mach-number computational fluid dynamics (CFD) formulations or evaluated directly from the transported variable
$k$
(Wilcox Reference Wilcox1998). However, the transported
$k$
is calibrated as part of the turbulence closure and need not reproduce the physical turbulent kinetic-energy profile required by the momentum and total-energy equations. It is known that some
$k-\epsilon$
formulations give a closer representation of the turbulent kinetic-energy distribution in wall-bounded flows (Lang & Shih Reference Lang and Shih1991; Biswas et al. Reference Biswas, Durbin and Medic2019). We therefore construct an algebraic surrogate
$\mathcal{K}$
from local RANS variables, with its profile calibrated toward a wall-resolved
$k-\epsilon$
reference. The original
$k$
and
$\omega _0$
transport equations remain unchanged, and the transported
$k$
continues to determine the eddy viscosity. The reconstructed
$\mathcal{K}$
is used only where a physical kinetic energy enters the modelled momentum and total-energy equations, including the kinetic-energy transport contribution that motivates the present closure.
The turbulent heat flux provides the complementary component of the closure and is represented within the scalar gradient-diffusion approximation through a local turbulent-Prandtl-number relation. Before specifying an explicit algebraic form, we infer a spatial distribution of
$\textit{Pr}_t$
using the regularised ensemble-Kalman method. The method is based on the ensemble-Kalman filter (EnKF) (Evensen Reference Evensen1994; Iglesias et al. Reference Iglesias, Law and Stuart2013) and follows its previous application to turbulence model field inversion (Zhang et al. Reference Zhang, Michelén-Ströfer and Xiao2020, Reference Zhang, Xiao and He2022). The inversion is performed in the compressible-channel flows of Modesti & Pirozzoli (Reference Modesti and Pirozzoli2016), for which complete DNS velocity and temperature profiles provide a one-dimensional construction problem constrained by both mean fields. An explicit relation based on local RANS variables is then fitted to the inferred
$\textit{Pr}_t$
profiles. The channel calculations isolate this turbulent heat-flux component, whereas the cooled zero-pressure-gradient boundary layers assess it together with the reconstructed
$\mathcal{K}$
as the complete energy-equation closure. The same algebraic closure is subsequently examined in the specifically defined
$\textit{SST}-\omega _0$
and
$\textit{EARSM}-\omega _0$
implementations based on Menter (Reference Menter1994) and Wallin & Jonasson (Reference Wallin and Jonasson2000), respectively, where
$\textit{EARSM}$
denotes explicit algebraic Reynolds stress model. All new relations use local variables and introduce no additional transport equation.
The paper is organised as follows. The modelling problem is outlined in § 2, including DNS-based physical analysis and a discussion of the RANS modelling framework. Section 3 introduces the numerical methods and details of the flow cases considered in this study. The TKE recovery model is presented in § 4. Section 5 focuses on the data-informed turbulent-Prandtl-number modelling, introducing the inversion algorithm and the validation results of the proposed model. The impact of the new energy modelling closure is also discussed here. Section 6 examines the extension of the closure to other
$k-\omega$
type models. Finally, the main findings are summarised in § 7.
2. Closure of the compressible energy equation
2.1. The energy equation to be closed
Wilcox (Reference Wilcox1998, pp. 243–254) introduced the compressible RANS energy equation, written as
\begin{equation} \begin{aligned} &\frac {\partial } {\partial t}\left [ \bar {\rho }\left (\tilde {e}+\frac {1}{2}\tilde {u}_i\tilde {u}_i\right )+\frac {\bar {\rho }\widetilde {u_i'' u_i''}}{2}\right ]+\frac {\partial }{\partial x_{\!j}}\left [\bar {\rho }\tilde {u}_{\!j}\left (\tilde {h}+\frac {1}{2}\tilde {u}_i \tilde {u}_i\right )+ \tilde {u}_{\!j} \frac {\bar {\rho }\widetilde {u_i'' u_i''}}{2}\right ]\\ &= \frac {\partial }{\partial x_{\!j}}\left [ -\tilde {q} - \widetilde {\rho h'' u_{\!j}''}+\widetilde {\tau _{ij}u_i''}-\frac {\widetilde {\rho u_i'' u_i''u_{\!j}''}}{2}+\tilde {u}_i(\tilde {\tau }_{ij} -\widetilde {\rho u_i'' u_{\!j}''})\right ]\!, \end{aligned} \end{equation}
where
$\bar {(\boldsymbol{\cdot })}$
and
$\tilde {(\boldsymbol{\cdot })}$
denote Reynolds- and Favre-averaging operators,
${\rho }$
is the density,
$e$
is the specific internal energy,
$h$
is the specific enthalpy,
$u''$
is the velocity fluctuation and single and double primes denote the corresponding fluctuations, respectively. The latter operator is defined as
$\tilde {\phi }=\overline {\rho \phi }/\bar {\rho }$
for an arbitrary variable
$\phi$
. For a Newtonian fluid, the stress tensor
$\tilde {\tau }_{ij}$
is given by the constitutive relation
\begin{eqnarray} &\tilde {\tau }_{ij} = 2\mu \tilde {S}_{ij} - \frac {2}{3}\mu \frac {\partial \tilde {u}_k}{\partial x_k}\delta _{ij}, \end{eqnarray}
\begin{eqnarray} &\tilde {S}_{ij}=\frac {1}{2} \left ( \frac {\partial \tilde {u}_i}{\partial x_{\!j}} + \frac {\partial \tilde {u}_{\!j}}{ \partial x_i} \right )\!, \end{eqnarray}
where
$\mu$
is the molecular dynamic viscosity,
$\tilde S_{ij}$
is the strain-rate tensor and
$\delta _{ij}$
is the Kronecker delta. In the current work, Sutherland’s law is used to prescribe viscosity, and the perfect-gas law is used. In (2.1),
$\tilde {q}$
denotes the molecular heat flux and is evaluated from Fourier’s law as
\begin{equation} \tilde {q} = -\lambda \frac {\partial \tilde {T}}{\partial x_{\!j}} = -\alpha \frac {\partial \tilde {h}}{\partial x_{\!j}}, \end{equation}
where
$\lambda$
is thermal conductivity and
$\alpha$
is thermal diffusivity,
$\alpha =\lambda /C_{\!p}$
, where
$C_p$
is the heat capacity. The laminar Prandtl number is then defined as
\begin{equation} \textit{Pr} = \frac {C_{\!p} \mu }{\lambda }. \end{equation}
Here,
$\textit{Pr}$
is the Prandtl number and is set to 0.71 in the current research to match reference cases in the literature and
$0.5\widetilde {u_i'' u_i''}$
represents the TKE.
Among the terms on the right-hand side of (2.1) that require modelling,
$-\widetilde {\rho h'' {u_{\!j}}''}$
, often referred to as the turbulent heat flux or turbulent enthalpy flux, is generally the dominant turbulent contribution to wall-normal thermal transport. By contrast, the turbulent transport of Reynolds stress, i.e.
$({\partial }/{\partial x_{\!j}}) (-0.5\widetilde {\rho u_i'' u_i''u_{\!j}''} )$
, is typically an order of magnitude smaller and is therefore often neglected. Note that, when solved together with eddy-viscosity turbulence models, the turbulent transport of the Reynolds stress is often simplified to the diffusion of TKE. Since this paper solves the Reynolds-averaged equations using eddy-viscosity turbulence models, the modelling problem will be referred to as TKE diffusion modelling in the following to avoid confusion. However, in hypersonic TBL with wall cooling, the thermal field exhibits unique characteristics distinct from those in weakly compressible cases. As a result of the competition between wall heat extraction and thermal production, the mean temperature profiles do not vary monotonically in the wall-normal direction, developing a peak value away from the wall. This temperature overshoot implies the existence of a wall-normal position at which the modelled turbulent heat flux vanishes, as shown in figure 1. The M6Tw025 case is given by Fan et al. (Reference Fan, Li and Pirozzoli2022), with free-stream Mach number
$M_{\infty }=5.84$
and recovery temperature ratio
$T_w/T_r=0.25$
, where
$T_w$
is the wall temperature and
$T_r$
is the recovery temperature. More detailed information on the case can be found in table 3.
The DNS results of the M6Tw025 case. (a) Temperature profile, (b) turbulent heat flux and turbulent transport of Reynolds-stress components in semi-local units. The semi-local scaling method employs the approach proposed by Trettel & Larsson (Reference Trettel and Larsson2016). The black dotted line marks the wall-normal location at which the turbulent heat flux vanishes, and the grey dash-dotted line denotes zero. Data are from Fan et al. (Reference Fan, Li and Pirozzoli2022).
Figure 1. Long description
Two line graphs depict temperature and turbulent heat flux in a hypersonic boundary layer. Panel A: A line graph shows the temperature profile. The x-axis is labeled y* on a logarithmic scale ranging from 10^0 to 10^3. The y-axis is labeled (T - T_w)/(T_w - T_∞) and ranges from -0.4 to 2.4. The temperature profile peaks around y* = 10 and then decreases. Panel B: A line graph shows turbulent heat flux and kinetic-energy transport. The x-axis is labeled y* on a logarithmic scale ranging from 10^0 to 10^3. The y-axis is labeled Turbulent heat flux and kinetic-energy transport (DNS) and ranges from -8 to 2. Two lines are plotted: a blue line representing -ρh’’u2’’ and a red dashed line representing -0.5ρu’’u’’2. The black dotted line marks the wall-normal location at which the turbulent heat flux vanishes, and the grey dash-dotted line denotes zero.
At the wall-normal location marked in figure 1, the turbulent heat flux passes through zero and therefore contributes little to the local wall-normal energy balance. The turbulent transport of kinetic energy, although smaller over most of the boundary layer, consequently becomes non-negligible in this region. Figure 2 compares these two turbulent transport contributions from DNS with their baseline RANS approximations. The
$k-\omega _0$
calculation shifts the zero crossing of the turbulent enthalpy flux towards the wall and predicts a kinetic-energy-diffusion term that differs from the DNS transport in both magnitude and wall-normal variation. These discrepancies occur in the same region in which the baseline calculation overpredicts the off-wall temperature maximum. Figure 2 alone, however, does not establish that either transport error causes the temperature error; rather, it identifies the two closure terms that require separate examination. The baseline turbulence and energy equations underlying these closures are introduced next, before the kinetic-energy reconstruction and turbulent-Prandtl-number relation are developed in §§ 4 and 5.
Predictions of the terms representing turbulent heat flux and kinetic-energy transport from the RANS calculation of M6Tw025 (Durbin & Yin Reference Durbin and Yin2026), compared with the DNS data: (a) full profiles and (b) an enlargement near the wall-normal location at which the turbulent heat flux vanishes. The black dotted line marks this location, and the grey dash-dotted line denotes zero. The DNS data are from Fan et al. (Reference Fan, Li and Pirozzoli2022).
2.2. Transport equations in the RANS framework
2.2.1. The
$k-\omega$
turbulence equations for compressible flow
The
$k-\omega$
two-equation models are widely used in engineering simulations. They solve two additional transport equations for a modelled kinetic-energy variable
$k$
and the specific-dissipation rate
$\omega$
, from which the eddy viscosity is constructed as
$\mu _t=\bar {\rho }k/\omega$
. The transported variable
$k$
is used as part of the turbulence model closure and is distinguished in this study from the physical TKE contribution required in the RANS energy equation; this distinction is discussed in § 4.
The proposed energy-diffusion model closure is evaluated alongside a compressible
$k-\omega$
model, which provides the momentum and eddy-viscosity fields. For the main calculations in this study, the specific-dissipation equation is written using the
$\omega _0$
reformulation proposed by Durbin & Yin (Reference Durbin and Yin2026). Built on the Wilcox
$k-\omega$
model (Wilcox Reference Wilcox1988), the reformulation replaces
$\omega$
with
$\omega _0$
and yields a more regular near-wall variable satisfying
$\omega _0=0$
at the wall; the mapping ensures a consistent high-Mach-number scaling of
$\omega$
. The transport equation of
$\omega _0$
is
\begin{equation} \begin{aligned} \frac {\partial }{\partial t}(\bar {\rho } \omega _0) + \frac {\partial }{\partial x_i}\left (\bar {\rho }\tilde {u}_{\!j} \omega _0 \right ) &= \frac {\partial }{\partial x_{\!j}}\left [(\mu +\sigma _{\omega }\mu _t)\frac {\partial \omega _0}{\partial x_{\!j}}\right ] + \frac {C_{\omega 1} \omega _0}{k}\tau ^t_{ij}\frac {\partial \tilde {u}_i}{\partial x_{\!j}} - C_{\omega 2}\bar {\rho }\omega _0^2,\\ \mathrm{with} \qquad X &= \frac {\omega _0 d^2 \bar {\rho }}{\mu }, \qquad \omega = \frac {C_{\omega 2}X}{6} + e^{-AX}, \end{aligned} \end{equation}
where
$A=0.012$
and
$d$
denotes the wall distance. Durbin & Yin (Reference Durbin and Yin2026) reported that this
$k-\omega _0$
model improved mean-velocity predictions with this formulation for the cooled Mach-6 and Mach-14 boundary layers considered there. We therefore use
$k-\omega _0$
as the momentum-model baseline for the channel and boundary-layer calculations in §§ 4 and 5.
2.2.2. Reynolds-averaged energy equation
Mo et al. (Reference Mo, Gao, Li and Jiang2026) suggested that compressible RANS models, especially for hypersonic simulations, should include TKE diffusion terms. The RANS energy equation with the modelling closure can then be written as
\begin{equation} \begin{aligned} &\frac {\partial } {\partial t}\left [ \bar {\rho }\left (\tilde {e}+\frac {1}{2}\tilde {u}_i\tilde {u}_i\right )+\bar {\rho }k\right ]+\frac {\partial }{\partial x_{\!j}}\left [\bar {\rho }\tilde {u}_{\!j}\left (\tilde {h}+\frac {1}{2}\tilde {u}_i \tilde {u}_i\right )+ \bar {\rho }\tilde {u}_{\!j} k\right ]\\ &= \frac {\partial }{\partial x_{\!j}}\left [\tilde {u}_{\!j}(\tau _{ij}+\tau ^t_{ij})+ \left(\frac {\mu }{\textit{Pr}}+\frac {\mu _t}{\textit{Pr}_t}\right)\frac {\partial \tilde {h}}{\partial x_{\!j}} +(\mu +\sigma _{k,E}\mu _t)\frac {\partial k}{\partial x_{\!j}}\right ]\!. \end{aligned} \end{equation}
When an eddy-viscosity turbulence model is solved, the Reynolds stress
$\tau ^t_{ij}$
is closed by the Boussinesq approximation. The turbulent heat flux is closed here using the scalar gradient-diffusion approximation. The turbulent Prandtl number
$\textit{Pr}_t$
relates the eddy diffusivity of enthalpy to the eddy viscosity, giving
\begin{equation} -\widetilde {\rho h'' u_{\!j}''} = \frac {\mu _t}{\textit{Pr}_t} \frac {\partial \tilde {h}}{\partial x_{\!j}}. \end{equation}
Similarly, for the diffusion of TKE, the molecular diffusion flux and the turbulent flux terms are approximated by the gradient of TKE. When
$k$
is available in a RANS model to represent TKE, the diffusion can be modelled as (2.9)
\begin{equation} \widetilde {\tau _{ij}u_i''} -\frac {\widetilde {\rho u_i'' u_i''u_{\!j}''}}{2} = (\mu + \sigma _{k,E}\mu _t)\frac {\partial k}{\partial x_{\!j}}, \end{equation}
where
$\sigma _{k,E}$
is an empirical coefficient. For the classic
$k-\omega$
model, Wilcox (Reference Wilcox1998) suggested that
$\sigma _{k,E}=0.5$
. However, in this study, it is set to 1.0, similar to the diffusion term in the
$k-\epsilon$
equations of Launder & Spalding (Reference Launder and Spalding1974). This will be explained in the modelling section of § 4.
3. Numerical setting
3.1. One-dimensional compressible-channel flows
A non-dimensional, Python-based one-dimensonal solver is developed to solve the Reynolds-averaged equations for fully developed channel flow, which has been validated in Yin (Reference Yin2026). The
$k-\omega _0$
model is solved in validation cases, setting
$\omega _{0,w}=0$
on the wall boundary. The temperature boundary condition is isothermal, and the wall value of
$ k$
is set to zero. A constant body force drives the flow. The first grid point is placed at
$y^+ \lt 0.1$
. The diffusion terms are discretised by a second-order central scheme. An under-relaxation strategy is adopted to ensure steady-state convergence. The wall temperature and body force are controlled to match the bulk Mach number
$Ma_b$
and the friction Reynolds number
$\textit{Re}_{\tau }$
in DNS data.
Two
$M_b=1.5$
compressible turbulent channel-flow cases provide a well-conditioned inverse problem at
$\textit{Re}_\tau =500\text{ and }1015$
due to their statistically one-dimensional geometry and complete DNS velocity and temperature profiles. Table 1 shows the key non-dimensional parameters. The corresponding DNS simulation results are given by Modesti & Pirozzoli (Reference Modesti and Pirozzoli2016). The compressible-channel-flow cases reported by Trettel & Larsson (Reference Trettel and Larsson2016), with bulk Mach number
$Ma_b$
ranging from 0.7 to 3.0 at different
$\textit{Re}_{\tau }$
, are chosen for model validation within the one-dimensional framework, listed in table 2.
Key non-dimensional parameters of compressible-channel-flow cases for EnKF inversion, listing
$\textit{Re}_{\tau }$
, semi-local friction Reynolds number
$\textit{Re}_{\tau }^*$
, bulk Reynolds number
$\textit{Re}_b$
,
$Ma_b$
, friction Mach number
$Ma_{\tau }$
and centre non-dimensional temperature
$\tilde {T}_c/T_w$
.
Table 1. Long description
A table comparing non-dimensional parameters of compressible-channel-flow cases. The table has 8 rows and 5 columns. Column headers are Case, Re_tau, Re*, Ma_b, and T_c/T_w. Row 1: Case, M0.7R400; Re_tau, 437; Re*, 396; Ma_b, 0.7; T_c/T_w, 1.082. Row 2: Case, M0.7R600; Re_tau, 1015; Re*, 677; Ma_b, 0.7; T_c/T_w, 1.082. Row 3: Case, M1.7R200; Re_tau, 322; Re*, 197; Ma_b, 1.7; T_c/T_w, 1.483. Row 4: Case, M1.7R400; Re_tau, 663; Re*, 406; Ma_b, 1.7; T_c/T_w, 1.481. Row 5: Case, M1.7R600; Re_tau, 972; Re*, 956; Ma_b, 1.7; T_c/T_w, 1.480. Row 6: Case, M3.0R200; Re_tau, 650; Re*, 208; Ma_b, 3.0; T_c/T_w, 2.487. Row 7: Case, M3.0R400; Re_tau, 1232; Re*, 396; Ma_b, 3.0; T_c/T_w, 2.486. Row 8: Case, M3.0R600; Re_tau, 1876; Re*, 601; Ma_b, 3.0; T_c/T_w, 2.491.
Key non-dimensional parameters of compressible-channel-flow cases for model validation, listing
$\textit{Re}_{\tau }$
,
$\textit{Re}_{\tau }^*$
,
$Ma_b$
and
$\tilde {T}_c/T_w$
.
Table 2. Long description
A table comparing non-dimensional parameters of compressible-channel-flow cases. The table has 7 columns and 6 rows. Column headers are Case, Ma∞, U∞ (m/s), ρ∞ (kg/m³), T∞ (K), Tw (K), Tw/Tr, -Bq, δj (mm), (xa - xi)/δj Lx ∗ Ly (m), Nx ∗ Ny. Row 1: Case, M2P5; Ma∞, 2.5; U∞ (m/s), 823.6; ρ∞ (kg/m³), 0.100; T∞ (K), 270.0; Tw (K), 568.0; Tw/Tr, 1.00; -Bq, 0.000; δj (mm), 4.0; (xa - xi)/δj Lx ∗ Ly (m), 0.70 ∗ 0.164; Nx ∗ Ny, 250 ∗ 150. Row 2: Case, M6TW025; Ma∞, 5.84; U∞ (m/s), 869.1; ρ∞ (kg/m³), 0.044; T∞ (K), 55.2; Tw (K), 97.5; Tw/Tr, 0.25; -Bq, 0.140; δj (mm), 1.3; (xa - xi)/δj Lx ∗ Ly (m), 0.31 ∗ 0.075; Nx ∗ Ny, 310 ∗ 200. Row 3: Case, M8TW076; Ma∞, 5.86; U∞ (m/s), 870.4; ρ∞ (kg/m³), 0.043; T∞ (K), 55.0; Tw (K), 300.0; Tw/Tr, 0.76; -Bq, 0.019; δj (mm), 13.8; (xa - xi)/δj Lx ∗ Ly (m), 2.32 ∗ 0.550; Nx ∗ Ny, 345 ∗ 200. Row 4: Case, M8TW048; Ma∞, 7.87; U∞ (m/s), 1155.1; ρ∞ (kg/m³), 0.026; T∞ (K), 51.8; Tw (K), 298.0; Tw/Tr, 0.48; -Bq, 0.060; δj (mm), 20.0; (xa - xi)/δj Lx ∗ Ly (m), 3.10 ∗ 0.824; Nx ∗ Ny, 310 ∗ 200. Row 5: Case, M11TW20; Ma∞, 10.90; U∞ (m/s), 1778.4; ρ∞ (kg/m³), 0.103; T∞ (K), 66.5; Tw (K), 300.0; Tw/Tr, 0.20; -Bq, 0.164; δj (mm), 3.7; (xa - xi)/δj Lx ∗ Ly (m), 1.50 ∗ 0.200; Nx ∗ Ny, 900 ∗ 250. Row 6: Case, M14TW018; Ma∞, 13.64; U∞ (m/s), 1882.2; ρ∞ (kg/m³), 0.017; T∞ (K), 47.4; Tw (K), 300.0; Tw/Tr, 0.18; -Bq, 0.187; δj (mm), 18.8; (xa - xi)/δj Lx ∗ Ly (m), 6.20 ∗ 1.045; Nx ∗ Ny, 510 ∗ 250.
Free-stream conditions, wall temperature, inflow boundary thickness of DNS simulation, streamwise sampling location, computational domain size, and grid resolution of the ZPG TBL cases.
$N_x$
and
$N_y$
denote the total number of grid points in the streamwise and wall-normal directions, respectively.
Table 3. Long description
A table comparing various parameters for two cases, CH15B and CH15C. The table has 2 rows and 7 columns. The columns are labeled as Case, Reτ, Reτ*, Reb, Ma∞, Maτ, and Tc/Tw. Row 1: Case, CH15B; Reτ, 500; Reτ*, 333; Reb, 15334; Ma∞, 1.5; Maτ, 0.072; Tc/Tw, 1.372. Row 2: Case, CH15C; Reτ, 1015; Reτ*, 677; Reb, 34000; Ma∞, 1.5; Maτ, 0.065; Tc/Tw, 1.368.
3.2. Two-dimensional compressible zero-pressure-gradient TBLs
The two-dimensional cases are simulated using the open-source CFD toolbox OpenFOAM-v2312 (Jasak Reference Jasak2009), with the density-based high-speed aerodynamic solver developed by Heyns et al. (Reference Heyns, Oxtoby and Steenkamp2014), which solves the coupled continuity, momentum and energy equations. The convective fluxes are evaluated using the AUSM+-up scheme (Liou Reference Liou2006). The resulting linear systems are solved using the generalised minimal residual method algorithm with a lower-upper symmetric Gauss-Seidel preconditioner. The discretised
$k-\omega _0$
model transport equations are solved separately.
A series of compressible zero-pressure-gradient (ZPG) TBL cases supported by DNS data from Zhang et al. (Reference Zhang, Duan and Choudhari2018) and Huang et al. (Reference Huang, Nicholson, Duan, Choudhari and Bowersox2020) is selected as validation for the new heat-transfer model. The DNS data for the M11Tw020 case are reported by Huang et al. (Reference Huang, Nicholson, Duan, Choudhari and Bowersox2020), and the other cases are from Zhang et al. (Reference Zhang, Duan and Choudhari2018). The M6Tw025 case reported by Fan et al. (Reference Fan, Li and Pirozzoli2022) corresponds to the case reported by Zhang et al. (Reference Zhang, Duan and Choudhari2018), with a more detailed database of higher-order moments. For the RANS simulation, the inflow boundary condition for the model variable
$k$
is specified using a turbulence intensity of
$Tu = 0.1\,\%$
. The Dirichlet condition for
$\omega$
is specified such that the eddy-viscosity ratio
$\mu _{t,\infty }/\mu _{\infty }$
is maintained at 0.009, following the recommendation of Rumsey (Reference Rumsey2010). The other computational settings, including free-stream conditions, wall temperature, computational domain sizes, leading-edge locations and the method for calculating the streamwise sampling location, are used for all cases, consistent with the study of Danis & Durbin (Reference Danis and Durbin2022), as shown in table 3.
$B_q$
represents non-dimensional wall-heat-transfer rate, given by
$B_q = q_w/(\rho _w C_{\!p} u_\tau T_w)$
and should be negative for cooled wall, where
$q_W$
is the wall heat flux.
The mesh generation uses a constant stretching ratio to ensure that the
$y^+$
of the first degree of freedom (DOF) near the solid wall remains below 0.6. To calculate the streamwise sampling location, the streamwise location
$x_i$
at which the boundary-layer thickness in the RANS simulation matches that in the DNS inflow is first identified. Then, the sampling location
$x_a$
is subsequently calculated as
$(x_a - x_i) / \delta _i$
from the DNS database.
It is worth noting that, although both the channel flow and ZPG TBL cases have isothermal wall boundary conditions, their thermal configurations differ. In the fully developed channel cases, since there is no free stream, wall cooling is not characterised by a recovery temperature ratio, and there is no temperature-profile peak near the wall. In the ZPG TBL cases, the flow develops over a single isothermal wall exposed to a prescribed free stream. The cooling strength is measured by
$T_w/T_r$
in table 3, and wall cooling results in a temperature peak near the wall. Throughout the study, a calorically perfect-gas assumption is adopted for all cases, following the settings of the reference DNS databases and the prior RANS studies used as baselines.
Comparison of DNS TKE, the transported variable
$k$
and the recovered
$\mathcal{K}$
in wall units,
$\textit{Re}_{\theta }\approx 2500.$
The DNS data are from Schlatter & Örlü (Reference Schlatter and Örlü2010) and lag elliptic blending (Lag EB)
$k-\epsilon$
data are from Biswas et al. (Reference Biswas, Durbin and Medic2019).
4. Turbulent kinetic-energy-diffusion model
In this section, TKE diffusion is modelled to reduce overprediction of the temperature peak, and explicit model formulae are provided. Validation is performed through simulations of compressible ZPG TBLs.
4.1. Modelling of kinetic-energy distribution and diffusion
Within the present eddy-viscosity and gradient-diffusion framework, the turbulent heat flux is modelled through the turbulent Prandtl number
$\textit{Pr}_t$
, whereas the turbulent transport of Reynolds stress is represented as TKE diffusion, as written in (2.9). In conventional compressible
$k-\omega$
calculations, the contribution of turbulent transport of Reynolds stress is often neglected or directly evaluated using the transported variable
$k$
(Wilcox Reference Wilcox1998). Although
$k$
is structurally associated with TKE, it is a calibrated transported variable that renders a correct eddy viscosity near the wall and therefore does not necessarily retain a clear physical meaning as a direct representation of TKE, as emphasised by Durbin & Bader (Reference Durbin and Bader2021). Its distribution can differ appreciably from the physical TKE distribution, particularly in the near-wall region. This observation motivates the construction of a mapping from the transported
$k$
to a TKE surrogate, denoted
$\mathcal{K}$
, that is more appropriate for the RANS momentum and energy equations.
As indicated by Lang & Shih (Reference Lang and Shih1991), Biswas et al. (Reference Biswas, Durbin and Medic2019) and illustrated in figure 3, the transported
$k$
of the
$k-\epsilon$
model provides a distribution closer to the DNS TKE profile than the original
$k$
predicted by the
$k-\omega _0$
model. Moreover, the transported variable
$k$
in both the
$k-\omega _0$
and
$k-\epsilon$
models follows inner-wall scaling with
$u_\tau ^2$
, whereas the physical TKE usually does not. A direct mapping from
$k$
to DNS TKE would therefore be case-dependent and would not provide a universal local closure. We instead construct a mapping from the
$k-\omega _0$
transported
$k$
to a
$k-\epsilon$
-like RANS surrogate
$\mathcal{K}$
.
It remains to select a local variable for the mapping. The non-dimensional eddy viscosity
$\nu _t^+$
is available from the RANS solution, follows the required inner-wall scaling and is a local quantity suitable for an explicit turbulence closure. It is therefore selected as the primary mapping variable. A standard high-Reynolds-number
$k-\epsilon$
model cannot be integrated reliably to the wall without additional near-wall treatment; accordingly, the wall-resolved Lag EB
$k-\epsilon$
model proposed by Billard & Laurence (Reference Billard and Laurence2012) is used as the reference for the
$k-\epsilon$
-like target profile. The resulting semi-empirical recovery relation is
\begin{equation} \begin{cases} \mathcal{K} = \left [ (1 - \phi _{\mathcal{K}}) \boldsymbol{\cdot }f_{\mathcal{K}} + \phi _{\mathcal{K}} \right ] \boldsymbol{\cdot }k ,\\[9pt] f_{\mathcal{K}} = \left [ \max \left ( 10^{-4}, \left ( 1 - \exp \left ( -\frac {\nu _t^+}{C_1} \right ) \right )^{0.25} \right ) \right ]^{-1}\!, \\[13pt] \phi _{\mathcal{K}} = \tanh \left ( \frac {\nu _t^+}{C_2} \right )\!, \end{cases} \end{equation}
where
$C_1=56$
,
$C_2=70$
are empirical constants. Model parameters were calibrated using an incompressible TBL case, shown in figure 3. This case fixes a RANS-compatible inner-scaled mapping rather than a universal TKE representation that matches DNS at all Reynolds numbers.
The reconstructed
$\mathcal{K}$
is not transported and does not enter either the turbulence transport equation or the eddy-viscosity relation. It replaces
$k$
only where a physical kinetic energy appears in the closed momentum and total-energy equations: the isotropic Reynolds-stress contribution, the kinetic-energy storage and convection terms and the modelled kinetic-energy diffusion. The transported
$k$
remains in the
$k-\omega _0$
equations and in
$\mu _t=\bar {\rho }k/\omega$
. With
$\sigma _{k,E}=1.0$
, the resulting closed equation is
\begin{equation} \begin{aligned} &\frac {\partial } {\partial t}\left [ \bar {\rho }\left (\tilde {e}+\frac {1}{2}\tilde {u}_i\tilde {u}_i\right )+\bar {\rho }\mathcal{K}\right ]+\frac {\partial }{\partial x_{\!j}}\left [\bar {\rho }\tilde {u}_{\!j}\left (\tilde {h}+\frac {1}{2}\tilde {u}_i \tilde {u}_i\right )+ \bar {\rho }\tilde {u}_{\!j} \mathcal{K}\right ]\\[3pt] &= \frac {\partial }{\partial x_{\!j}}\left [\tilde {u}_{\!j}(\tau _{ij}+\tau ^t_{ij})+ \left(\frac {\mu }{\textit{Pr}}+\frac {\mu _t}{\textit{Pr}_t}\right)\frac {\partial \tilde {h}}{\partial x_{\!j}} +(\mu +\sigma _{k,E}\mu _t)\frac {\partial \mathcal{K}}{\partial x_{\!j}}\right ]\!,\\[3pt] {\tau }^t_{ij} &= 2\mu _t \tilde {S}_{ij} - \frac {2}{3}\bar {\rho }\mathcal{K}\delta _{ij}. \end{aligned} \end{equation}
Note that the Reynolds-stress term in the
$k$
equation retains its original form to preserve eddy-viscosity calibration.
4.2. Validation of the turbulent kinetic-energy-diffusion model
Preliminary tests indicate that the TKE diffusion that uses the
$\mathcal{K}$
from (4.2) has no significant impact on the one-dimensional channel-flow cases mentioned. This is expected because, unlike the cooled boundary layers, these channel cases have no off-wall location at which the turbulent enthalpy flux vanishes in association with a temperature maximum. Consequently, validation of the recovery formula is primarily conducted within two-dimensional wall-bounded turbulent flow problems. The commonly used constant value
$\textit{Pr}_t=0.9$
is adopted.
The algebraic recovery relation in (4.1) and its implementation in (4.2) are first examined in M6Tw025, as shown in figure 4. Figure 4(a) shows that the eddy-viscosity field remains close to that of the baseline calculation, consistent with the fact that the transported
$k$
, rather than
$\mathcal{K}$
, is retained in the eddy-viscosity relation. Figure 4(b) shows that the transported
$k$
of the baseline
$k$
–
$\omega _0$
model does not reproduce the near-wall peak of the DNS TKE, whereas the reconstructed
$\mathcal{K}$
provides a closer representation of its peak location and inner-layer variation. Figure 4(c) further shows that the
$\mathcal{K}$
-based gradient-diffusion approximation is closer to the DNS kinetic-energy transport near the off-wall location where the turbulent enthalpy flux changes sign. This comparison provides a diagnostic assessment of the reconstructed transport term but does not, by itself, establish its effect on the mean thermal field. The corresponding a posteriori effects of using
$\mathcal{K}$
in the physical kinetic-energy terms of the momentum and total-energy equations are assessed in figures 5 and 6 and quantified in table 5.
Effect of the recovery model in the M6Tw025 case. (a) Comparison of
$\nu _t^+$
, (b) comparison of DNS TKE, the transported variable
$k$
and the recovered
$\mathcal{K}$
in wall units, (c) comparison of turbulent kinetic-energy transport contribution in semi-local wall units. For the DNS data, TKE profile is given by Zhang et al. (Reference Zhang, Duan and Choudhari2018) and turbulent kinetic-energy transport contribution is given by Fan et al. (Reference Fan, Li and Pirozzoli2022).
Figure 4. Long description
Panel A: A line graph compares the variable vt+ against y* on a logarithmic scale. The graph includes three data series: DNS (represented by circles), k (solid black line), and K (dashed blue line). The lines closely follow the trend of the DNS data points, showing a peak around y* = 100. Panel B: Another line graph compares TKE+, k+, and K+ against y* on a logarithmic scale. The graph includes three data series: DNS TKE+ (circles), k+ (solid black line), and K+ (dashed blue line). The lines and data points show a peak around y* = 10, with K+ having a higher peak than k+. Panel C: A line graph compares the turbulent kinetic-energy transport contribution against y* on a logarithmic scale. The graph includes three data series: DNS (circles), σkE(μt∂k+/μw∂y+) (solid black line), and -0.5ρu′′iu′′j (dashed blue line). The lines and data points show a peak around y* = 10 and a trough around y* = 100.
Figures 5 and 6 show the a posteriori effect of using
$\mathcal K$
in the physical kinetic-energy terms of the momentum and total-energy equations for M6Tw025 and M14Tw018. In both cases, the reconstruction improves the temperature profile between the wall and the temperature maximum while leaving the predicted velocity profile close to that of the baseline momentum model. Because
$\mathcal K$
enters the isotropic Reynolds-stress contribution as well as the kinetic-energy storage, convection and diffusion terms, these comparisons quantify the combined effect of the reconstruction rather than the diffusion term in isolation.
Beyond the temperature maximum and into the wake region, the improvement is limited in both cases. This limitation is also consistent with the model construction, as in the outer layer, the wall-normal energy transport is dominated by turbulent heat flux, so a correction to the
$\mathcal{K}$
diffusion term alone should not be expected to remove the remaining thermal error.
Effect of the TKE diffusion closure in M6Tw025 case. (a) Velocity profile, (b) temperature profile in wall units. The DNS data are from Zhang et al. (Reference Zhang, Duan and Choudhari2018).
Effect of TKE diffusion closure in M14Tw018 case. (a) Velocity profile, (b) non-dimensional temperature profile in wall units. The DNS data are from Zhang et al. (Reference Zhang, Duan and Choudhari2018).
The two cooled-wall cases exhibit different sensitivities to the reconstruction. Because the Reynolds number, Mach number, cooling ratio, property variation and baseline-model error change simultaneously between the cases, the present comparisons do not isolate the parameter responsible for this difference. From the turbulence modelling perspective, the authors emphasise that, because the
$\mathcal{K}$
recovered from (4.1) scales with
$u_\tau ^2$
, it should be self-consistent across different Reynolds numbers. However, such scaling does not hold for the TKE in DNS simulations. Thus, the effectiveness of the TKE recovery formula may vary with the Reynolds numbers. Nevertheless, the difference indicates that the TKE diffusion term exerted a slightly greater influence on heat transfer in M6Tw025 than in M14Tw018, a finding corroborated by the higher peak temperatures observed in the DNS simulation. A broader parametric assessment is therefore required to establish the range over which the recovery relation remains effective.
5. Compressible turbulent-Prandtl-number model
In this section, we describe the methodology of regularised ensemble Kalman filter (REnKF)-inspired
$\textit{Pr}_t$
modelling. The resulting explicit fitting formulae approximate the inversion results. The new
$\textit{Pr}_t$
model is first applied to one-dimensional compressible-channel flows, alongside other existing models for comparison. Subsequently, it is used together with the TKE diffusion model proposed in § 4 to complete the energy-equation closure. The closure is tested in two-dimensional ZPG TBL flows.
5.1. Inferring the Prandtl-number distribution with the REnKF method
5.1.1. The REnKF inversion algorithm
The purpose of the REnKF step is to infer the turbulent-Prandtl-number field from the thermal response, rather than to allow
$\textit{Pr}_t$
to compensate for errors in the momentum closure. Because the standard Wilcox
$k-\omega$
model does not reproduce the DNS mean-velocity profile exactly, residual momentum error could otherwise contaminate the inferred
$\textit{Pr}_t$
. An auxiliary spatial field
$f_c(y)$
is therefore introduced only during the inversion. It multiplies the production and destruction terms in the
$\omega$
equation
\begin{equation} \begin{aligned} \frac {\partial }{\partial t}\left (\bar {\rho }\omega \right ) +\frac {\partial }{\partial x_{\!j}} \left (\bar {\rho }\widetilde {u}_{\!j}\omega \right ) ={}& \frac {\partial }{\partial x_{\!j}} \left [ \left (\mu +\sigma _{\omega }\mu _t\right ) \frac {\partial \omega }{\partial x_{\!j}} \right ] \\ &+f_c\frac {C_{\omega 1}\omega }{k} \tau ^t_{ij} \frac {\partial \widetilde {u}_i}{\partial x_{\!j}} -f_c C_{\omega 2}\bar {\rho }\omega ^2 . \end{aligned} \end{equation}
During the inversion,
$f_c(y)$
and
$\textit{Pr}_t(y)$
are inferred jointly from the DNS velocity and temperature profiles. The field
$f_c$
controls the momentum-side discrepancy, allowing the thermal response to constrain
$\textit{Pr}_t$
with less contamination from velocity profile error. A related use of a spatial correction to the
$\omega$
equation for compressible-TBL velocity profiles was reported by Danis & Durbin (Reference Danis and Durbin2022). The near-wall value is prescribed as
$\omega =60\nu /(C_{\omega 2}y_1^2)$
, where
$y_1$
is the distance from the wall to the centre of the first computational cell. Importantly,
$f_c$
is an auxiliary inversion variable rather than a predictive compressibility correction or a component of the proposed energy-equation closure. It is discarded after the inversion and is used neither in the channel-flow validation nor in any boundary-layer calculation. All forward validation calculations instead use the
$\omega _0$
(2.6) without an
$f_c$
field.
The inversion algorithm is based on the REnKF framework (Ströfer et al. Reference Ströfer, Zhang and Xiao2020; Zhang et al. Reference Zhang, Michelén-Ströfer and Xiao2020). EnKF performs sequential, gradient-free data assimilation via implicit cost-function minimisation, with REnKF extending the algorithm by incorporating regularisation to constrain the inversion. The cost function
$J$
is given as
\begin{equation} J(\boldsymbol {x}_{\!j}) = ||\boldsymbol {x}_{\!j}^a - \boldsymbol {x}_{\!j}^f||^2_{P^{-1}} + ||H\boldsymbol {x}_{\!j}^a-\boldsymbol {y}_{\!j}||^2_{R^{-1}}+\lambda _r||\mathcal{G}[\boldsymbol {x}_{\!j}^a]||^2_{Q^{-1}}, \end{equation}
where
$\boldsymbol{x}$
denotes the state vector to be inferred and
$\boldsymbol{y}$
denotes the observation vector; subscripts
$f$
and
$a$
indicate the forecast and analysis states, respectively;
$H$
is the observation operator mapping
$\boldsymbol{x}$
to the observation space, while
$P$
and
$R$
are the covariance matrices of the state and observation;
$\mathcal{G}$
represents a differentiable regularisation function, with
$Q$
its associated covariance;
$\lambda _r$
is the regularisation weighting coefficient.
The DNS velocity and temperature profiles are both included in the observation vector. The inversion state comprises
$f_c(y)$
and
$\textit{Pr}_t(y)$
, which are inferred jointly. The role of
$f_c$
is limited to controlling the momentum-side discrepancy during inversion so that the thermal observations primarily constrain
$\textit{Pr}_t$
. This inversion uses the transported
$k$
and does not include the algebraic reconstruction
$\mathcal{K}$
. No algebraic model is fitted to the inferred
$f_c$
field, and the field is discarded before all forward validation calculations.
Although channel flow and TBL cases have different global thermal balances, the construction of the present
$\textit{Pr}_t$
model is mainly aimed at the logarithmic region of the boundary layer, where the modelled turbulent heat flux is governed by local eddy diffusivity rather than by the outer thermal boundary condition alone. The channel geometry provides complete one-dimensional velocity and temperature profiles and hence a better-conditioned inversion than a developing boundary layer sampled at a few stations. This advantage does not guarantee transfer between the two thermal configurations. We therefore treat the channel cases as construction data and the cooled ZPG boundary layers as out-of-configuration transfer tests of the fitted local relation.
For the REnKF process, the initial
$f_c$
field is supplied by the model of Yin (Reference Yin2026), whereas the initial
$\textit{Pr}_t$
field is selected independently from the tested turbulent-Prandtl-number relations. Preliminary inversions starting from different
$\textit{Pr}_t$
profiles produce a common peak near
$y^*\approx 20$
. The Ceci et al. (Reference Ceci, Pirozzoli and Spalart2026) relation provides the closest initial representation of this feature and is therefore used as the regularisation baseline. The initial
$f_c$
field does not generate or otherwise determine the initial
$\textit{Pr}_t$
field.
5.1.2. The
$\textit{Pr}_t$
model construction based on inferred distribution
The resulting velocity and temperature profiles from the inversion process are shown in figure 7. The inferred fields bring both channel-flow velocity and temperature profiles closer to their DNS targets than the initial fields, with the remaining differences concentrated toward the channel centreline.
Inferred results of compressible-channel-flow cases. (a) Velocity profile, (b) temperature profile in wall units. Here,
$u^+ = u/u_{\tau }$
and
$T^+=T/T_w$
. Blue denotes the CH15B case, and red denotes the CH15C case. The DNS simulation data are reported by Modesti & Pirozzoli (Reference Modesti and Pirozzoli2016).
The inversion provides a spatially varying
$\textit{Pr}_t$
field, which must be converted into an algebraic relation based on quantities available in a RANS calculation. We therefore parameterise
$\textit{Pr}_t$
in terms of the locally scaled eddy viscosity,
$\nu _t^+=\nu _t/\nu$
. This choice is motivated by the semi-local behaviour of
$\textit{Pr}_t$
and
$\nu _t^+$
reported by Huang et al. (Reference Huang, Coleman, Spalart and Yang2023) and by the related use of
$\nu _t^+$
in the turbulent-diffusivity formulation of Gupta & Duraisamy (Reference Gupta and Duraisamy2026). The two channel-flow inversions are compared in semi-local coordinates to assess whether they can be represented approximately by a common relation.
Based on that, the regularisation parameter
$\lambda _r$
is tuned to drive the CH15B and CH15C cases toward a normalised distribution that is both consistent between cases and physically consistent with free-stream recovery, as shown in figure 8.
Construction of the new
$\textit{Pr}_t$
model. Blue denotes the CH15B case, and red denotes the CH15C case.
Then an ad hoc empirical
$\textit{Pr}_{t}$
formulation for compressible turbulent channel flows may be constructed as
\begin{eqnarray} &\textit{Pr}_{t}^{\textit{new}} = \textit{Pr}_{t,\infty } + (B_1\nu _t^+ - B_2)\boldsymbol{\cdot }\mathrm{exp}\big(-B_3{\nu _t^+}^{B_4}\big), \end{eqnarray}
where
$\textit{Pr}_{t,\infty }= 0.843$
, as suggested by Ceci et al. (Reference Ceci, Pirozzoli and Spalart2026) as an air-specific free-stream value,
$B_1 =0.395$
,
$B_2=0.35$
,
$B_3=0.49$
,
$B_4 = 0.75$
. Although the mathematical form looks cumbersome, its rationale is straightforward. The
$(B_1\nu _t^+ - B_2)$
term creates a linear line that intersects the
$\nu _t^+ = 0$
axis with
$\textit{Pr}_{t}^{\textit{new}}=0.5$
, and quickly increases the
$\textit{Pr}_t$
value in the buffer layer. The factor
$\mathrm{exp}(-B_3{\nu _t^+}^{B_4})$
then smoothly damps the expression towards
$\textit{Pr}_{t,\infty }$
in the logarithmic layer.
However, for a realistic turbulent boundary layer (TBL), (5.3) cannot distinguish the inner boundary layer from the free stream. Thus, an auxiliary parameter
$\phi _{\textit{Pr}_t}$
is introduced to detect the outer boundary of a TBL, calibrated through the M6Tw025 case as shown in figure 9. The result is the complete form
\begin{eqnarray} &\textit{Pr}_{t}^{\textit{new}} = \textit{Pr}_{t,\infty } + \phi _{\textit{Pr}_t}\boldsymbol{\cdot }(B_1\nu _t^+ - B_2)\boldsymbol{\cdot }\mathrm{exp}\big(-B_3{\nu _t^+}^{B_4}\big), \end{eqnarray}
\begin{eqnarray} &\phi _{\textit{Pr}_t} = \mathrm{min}\left(1,\left(\frac {B_5\nu }{\omega d^2}\right)\right)^{B_6}, \end{eqnarray}
where
$B_5 = 2000$
,
$B_6 = 4$
and
$d$
is the wall distance. The blending factor
$\phi _{\textit{Pr}_t}$
originates from the additional function
$F_1$
of the
$k-\omega$
SST model (Menter Reference Menter1994). In a plane turbulent channel flow, (5.4) reduces to (5.3) as
$\phi = 1$
. The present model, which uses (5.4) and (5.5), is referred to thereafter as the new
$\textit{Pr}_t$
model.
Distribution of
$\phi _{\textit{Pr}_t}$
and
$\textit{Pr}_t$
across the boundary layer of M6Tw025 case. Black dotted line indicates the edge of HTBL.
For air, the most widely used option is to use a constant
$\textit{Pr}_t$
value, which is usually
$0.9$
for air (referred to as constant). Several
$\textit{Pr}_t$
models (Kays & Crawford Reference Kays and Crawford1993; Huang et al. Reference Huang, Coleman, Spalart and Yang2023; Ceci et al. Reference Ceci, Pirozzoli and Spalart2026) based on local variables are also employed along with the present model. The detailed formulations of those models are shown in Appendix A.
5.2. Validation of the inferred Prandtl-number model
5.2.1. One-dimensional compressible-channel-flow cases
A series of compressible-channel-flow cases with available DNS reference data is used to assess the turbulent-Prandtl-number relations within the common
$k-\omega _0$
framework. For the cases reported by Trettel & Larsson (Reference Trettel and Larsson2016), the available thermal reference quantity is the centreline temperature. The relative error in the predicted centreline temperature is therefore adopted as the evaluation metric in figure 10.
Relative error in the centreline temperature predicted using the tested
$\textit{Pr}_t$
relations for the compressible-channel-flow cases at (a)
$M_b=0.7$
, (b)
$M_b=1.7$
and (c)
$M_b=3.0$
. The DNS data are from Trettel & Larsson (Reference Trettel and Larsson2016).
The ability of these cases to discriminate among the
$\textit{Pr}_t$
relations depends on both the Reynolds number and the Mach number. At
$\textit{Re}_\tau ^*\approx 200$
, represented by M1.7R200 and M3.0R200, the viscous, buffer and logarithmic regions are insufficiently separated for a clear assessment of the intended wall-layer behaviour. These cases are therefore retained in figure 10 for completeness but are not used as the primary basis for comparing the models. As the Mach number decreases, the thermal variation and the spread among the predicted centreline temperatures also decrease. The
$M_b=3.0$
cases at
$\textit{Re}_\tau ^*\approx 400$
and
$600$
consequently provide the most discriminating comparison within the available channel-flow dataset.
The tested
$\textit{Pr}_t$
relations were originally developed using different turbulence closures, calibration flows and dependent variables. Their implementation here with the same
$k-\omega _0$
momentum closure therefore assesses their transferred performance within the present calculation rather than their intrinsic accuracy in their original formulations. For the two most discriminating cases,
$M_b=3.0$
at
$\textit{Re}_\tau ^*\approx 400$
and
$600$
, the present relation gives among the smaller centreline-temperature errors shown in figure 10. The channel-flow comparison therefore supports that the present implementation helps, but does not show general superiority over the existing
$\textit{Pr}_t$
relations.
5.2.2. Two-dimensional compressible ZPG TBL cases
The ZPG boundary layers test the fitted
$\textit{Pr}_t$
relation in a developing external flow with wall cooling, outside the channel-flow construction configuration. The DNS-derived turbulent Prandtl number can be computed by (5.6)
\begin{equation} \textit{Pr}_t = \frac {\widetilde {\rho {u_1}''{u_2}''} \frac {\partial \tilde {h}}{\partial y}}{\widetilde {\rho h''{u_2}''}\frac {\partial \tilde {u_1}}{\partial y}}. \end{equation}
In the strongly cooled boundary layers, the DNS-implied
$\textit{Pr}_t$
becomes ill conditioned near the wall-normal location at which the turbulent heat flux or the mean temperature gradient passes through zero. Besides the new
$\textit{Pr}_t$
model, the aforementioned models (Kays & Crawford Reference Kays and Crawford1993; Huang et al. Reference Huang, Coleman, Spalart and Yang2023; Ceci et al. Reference Ceci, Pirozzoli and Spalart2026) poorly capture the gradient at this location, exerting only a limited influence on the prediction of the temperature profile compared with the constant model. Therefore, the study of two-dimensional ZPG TBL cases focuses mainly on the effect of the new
$\textit{Pr}_t$
model compared with
$\textit{Pr}_t=0.9$
.
Two cases with a non-monotonic temperature distribution, i.e. M6Tw025 and M14Tw018 (Zhang et al. Reference Zhang, Duan and Choudhari2018), are chosen for validation. The results are shown in figures 11 and 12, respectively. In both cooled boundary layers, the
$\mathcal{K}$
reconstruction provides the dominant reduction in the overpredicted temperature maximum. Adding the new
$\textit{Pr}_t$
relation further reduces the temperature peak. In M6Tw025, it changes the peak value from 2.236 to 2.204, compared with the DNS value 2.206, and changes
$-B_q$
from 0.139 to 0.140, compared with 0.140. In M14Tw018, it changes the peak from 1.885 to 1.849, compared with 1.849, but changes
$-B_q$
from 0.200 to 0.210, whereas the DNS value is 0.187. Thus, the
$\textit{Pr}_t$
component improves some temperature-profile measures, although its effect is secondary and case-dependent.
Effect of present
$\textit{Pr}_t$
model on M6Tw025 case. (a) Velocity profile, (b) non-dimensional temperature profile, (c) turbulent heat flux, (d)
$\textit{Pr}_t$
profile in wall units. The black dotted line marks the zero turbulent heat-flux point and the grey dash-dotted line denotes zero. The DNS data are from Zhang et al. (Reference Zhang, Duan and Choudhari2018) and Fan et al. (Reference Fan, Li and Pirozzoli2022).
Figure 11. Long description
Panel A: A line graph shows the velocity profile with the y-axis labeled u+ and the x-axis labeled y+ on a logarithmic scale. The graph includes data from DNS, k model with Prandtl number 0.9, K model with Prandtl number 0.9, and K model with a new Prandtl number. Panel B: A line graph displays the non-dimensional temperature profile with the y-axis labeled (T-T∞)/(Tw-T∞) and the x-axis labeled y+ on a logarithmic scale. The graph includes data from DNS, k model with Prandtl number 0.9, K model with Prandtl number 0.9, and K model with a new Prandtl number. Panel C: A line graph illustrates the turbulent heat flux with the y-axis labeled Turbulent heat flux term and the x-axis labeled y+. The graph includes data from DNS, (μt/0.9)(∂h+/∂y+)(K), (μt/PrTnew)(∂h+/∂y+)(K), and (μt/0.9)(∂h+/∂y+)(k). The black dotted line marks the zero turbulent heat-flux point and the grey dash-dotted line denotes zero. Panel D: A line graph shows the profile in wall units with the y-axis labeled PrT and the x-axis labeled y+ on a logarithmic scale. The graph includes data from DNS, K model with Prandtl number 0.9, and K model with a new Prandtl number.
The component-wise comparisons in figures 11 and 12 show that the two modifications produce different local responses. The
$\mathcal K$
reconstruction provides the larger change in the temperature profile between the wall and the temperature maximum. In M6Tw025, figure 11(c) shows that replacing
$\textit{Pr}_t=0.9$
by the local relation increases the magnitude of the modelled turbulent heat-flux term on the near-wall side of its zero crossing. This pointwise change is accompanied by a smaller additional reduction in the predicted temperature maximum.
To examine the sensitivity of the proposed closure to wall-normal resolution, M6Tw025 and M14Tw018 are recomputed using the refined meshes listed in table 4. The corresponding velocity and temperature profiles are compared in figure 13.
Grid settings and wall parameters of M6Tw025 and M14Tw018 cases with proposed modelling closure, showing the number of DOF in the wall-normal direction,
$y^+$
of the first DOF nearest to the wall, friction Reynolds number
$\textit{Re}_{\tau }$
, semi-local friction Reynolds number
$\textit{Re}_{\tau }^*$
, friction velocity
$u_{\tau }$
and wall-heat-transfer rate
$B_q$
.
Table 4. Long description
The table presents grid settings and wall parameters for various cases, including M6Tw025, M14Tw018, M2p5, M6Tw076, M8Tw048, and M11Tw020. It has 6 rows and 10 columns. The columns are labeled as Case, Model, Re_theta, Re_tau, Re_tau*, theta, delta (mm), u_tau (m/s), (T_w-T_e)/(T_e-T_inf)_max, and -B_q. The table includes different models such as DNS, k, Pr_t = 0.9, K, Pr_t = 0.9, K, new Pr_t, and their corresponding values for each parameter. The cases are grouped under different headings, and each row provides specific values for the parameters listed in the columns.
Effect of present
$\textit{Pr}_t$
model on M14Tw018 case. (a) Velocity profile, (b) non-dimensional temperature profile, (c)
$\textit{Pr}_t$
profile in wall units. The black dotted line marks the zero turbulent heat-flux point. The DNS data are from Zhang et al. (Reference Zhang, Duan and Choudhari2018).
Figure 12. Long description
Three line graphs compare velocity and temperature profiles in a HTBL. Panel A: A line graph shows the velocity profile with the y-axis labeled u+ and the x-axis labeled y+. The graph includes data points from DNS and three different models: k, Pr_t = 0.9; χ, Pr_t = 0.9; and χ, new Pr_t. Panel B: A line graph displays the non-dimensional temperature profile with the y-axis labeled (T - T∞)/(Tw - T∞) and the x-axis labeled y+. This graph also includes data points from DNS and the same three models as in Panel A. Panel C: A line graph presents the profile in wall units with the y-axis labeled Pr_t and the x-axis labeled y+. The black dotted line marks the zero turbulent heat-flux point. The graph includes data points from DNS and the same three models as in the previous panels.
For M6Tw025, the two meshes give nearly indistinguishable velocity and temperature profiles, and the reported value of
$-B_q$
remains 0.140. For M14Tw018, the temperature profile changes only slightly and
$-B_q$
decreases from 0.210 to 0.207, whereas the velocity profile and the corresponding values of
$\textit{Re}_{\tau }$
and
$u_{\tau }$
show a more visible response to refinement. The thermal predictions are therefore only weakly sensitive to the particular wall-normal refinement examined. Because only two resolutions are considered and a momentum-side change remains visible in M14Tw018, these results should not be interpreted as a general demonstration of grid independence.
The complete energy-equation closure is further evaluated in the additional compressible ZPG TBLs listed in table 5. The corresponding velocity and temperature profiles are shown in figures 14 and 15. Relative to the
$k,\,\textit{Pr}_t=0.9$
baseline, the complete closure reduces the error in the temperature maximum for M6Tw025, M6Tw076, M8Tw048, M11Tw020 and M14Tw018, while producing little change in M2p5. The velocity profiles remain close to those of the
$k-\omega _0$
momentum-model baseline. Changes in
$\textit{Re}_{\tau }^*$
are less uniform: the prediction moves towards the DNS value in M6Tw076, M8Tw048 and M14Tw018 but not in every case. The proposed closure should therefore be interpreted primarily as a modification of the thermal prediction rather than as a general correction to the momentum field.
The TBL parameters given by RANS simulations, listing the momentum thickness Reynolds number
$\textit{Re}_{\theta }$
, friction Reynolds number
$\textit{Re}_{\tau }$
, semi-local friction Reynolds number
$\textit{Re}_{\tau }^*$
, momentum thickness
$\theta$
, boundary-layer thickness
$\delta$
, friction velocity
$u_{\tau }$
, the temperature maximum and wall-heat-transfer rate
$B_q$
. The non-dimensionalisation of
$B_q$
is performed using DNS reference quantities.
Table 5. Long description
A table comparing parameters from RANS simulations for different cases and grid levels. The table has 2 rows and 8 columns. Column headers are Case, Grid level, N subscript y, y superscript plus, Re subscript tau, Re subscript tau superscript x, u subscript tau (m/s), and negative B subscript q. Row labels are M6Tw025 and M14Tw018, each with sub-rows for original and fine grid levels. Row 1: M6Tw025, original, 200, 0.15, 421, 887, 33.2, 0.140. Row 2: M6Tw025, fine, 240, 0.09, 422, 906, 33.2, 0.140. Row 3: M14Tw018, original, 250, 0.15, 750, 5644, 66.5, 0.210. Row 4: M14Tw018, fine, 300, 0.09, 726, 5694, 65.9, 0.207.
Effect of mesh refinement on M6Tw025 and M14Tw018 cases with proposed modelling closure. (a) Velocity profile, (b) non-dimensional temperature profile in wall units. The DNS data are from Zhang et al. (Reference Zhang, Duan and Choudhari2018).
Figure 13. Long description
Two line graphs compare velocity and temperature profiles for different mesh refinements. Panel A: The line graph shows the velocity profile with the vertical axis labeled u+ and the horizontal axis labeled y+. The graph includes data series for M6Tw025 and M14Tw018 cases, both with original and fine grid refinements, and DNS data. The velocity profiles show an increasing trend with y+ and converge at higher y+ values. Panel B: The line graph shows the non-dimensional temperature profile in wall units with the vertical axis labeled (T - T∞)/(Tw - T∞) and the horizontal axis labeled y+. The graph includes data series for M6Tw025 and M14Tw018 cases, both with original and fine grid refinements, and DNS data. The temperature profiles show a peak and then a decreasing trend with y+.
The response of the reported wall-heat-transfer measure also varies by case. The complete closure moves
$-B_q$
closer to the DNS value in M6Tw025, M8Tw048 and M14Tw018, gives a comparable absolute error in M6Tw076 and moves it slightly away from the DNS value in M11Tw020; the values in M2p5 are close to zero for both DNS and the complete closure. Moreover, in M14Tw018, the
$\mathcal{K},\,\textit{Pr}_t=0.9$
result remains closer to the DNS wall-heat-transfer value than the complete closure. Taken together, the additional boundary-layer cases provide evidence of improved temperature maximum and temperature-profile shape, and also some improvements in wall heat flux.
Velocity profiles of other ZPG TBL cases: (a) M2p5, (b) M6Tw076, (c) M8Tw048, (d) M11Tw020. The DNS data are from Zhang et al. (Reference Zhang, Duan and Choudhari2018) and Huang et al. (Reference Huang, Nicholson, Duan, Choudhari and Bowersox2020).
Non-dimensional temperature profiles of other ZPG TBL cases: (a) M2p5, (b) M6Tw076, (c) M8Tw048, (d) M11Tw020. The DNS data are from Zhang et al. (Reference Zhang, Duan and Choudhari2018) and Huang et al. (Reference Huang, Nicholson, Duan, Choudhari and Bowersox2020).
Figure 15. Long description
Four line graphs depict non-dimensional temperature profiles for various zero pressure gradient TBL cases. Each graph represents a different case: (a) M2p5, (b) M6Tw076, (c) M8Tw048, and (d) M11Tw020. The x-axis is labeled y+ on a logarithmic scale ranging from 10^0 to 10^3, and the y-axis is labeled (T - T∞)/(Tw - T∞) ranging from -0.2 to 1.4. Each graph includes three data series: DNS data represented by black circles, k, Prt = 0.9 represented by a solid black line, and a new Prt represented by a dashed red line. The graphs show how the temperature profiles vary across different cases, with the DNS data points providing a reference for the other two models. The trends indicate that the new Prt model generally follows the DNS data more closely than the k, Prt = 0.9 model, especially in the higher y+ ranges.
6. Applicability to the
$\boldsymbol{k-\omega}$
model family
To assess the transferability of the proposed energy-equation closure, it is tested within two
$k-\omega$
-based RANS family models: the SST model of Menter (Reference Menter1994), and the EARSM of Wallin & Jonasson (Reference Wallin and Jonasson2000) built on the
$k-\omega$
formulation recalibrated by Hellsten (Reference Hellsten2005). To keep the numerical wall treatment and the momentum-field baseline as consistent as possible across all transfer tests, the same
$k-\omega _0$
reformulation (2.6) is also introduced into these standard models. Apart from this change, their blending functions, stress closures and model constants are otherwise left unchanged; all the occurrences of
$\omega$
in these models are evaluated from the recovered
$\omega _0$
. For simplicity, these cases are therefore referred to as
$\textit{SST}-\omega _0$
and
$\textit{EARSM}-\omega _0$
. This terminology is used to avoid confusing these baselines with unmodified standard SST or EARSM models. The cases M6Tw025 and M14Tw018 are chosen as test cases.
Figures 16 and 17 examine the closure in the specifically defined
$\textit{SST}-\omega _0$
and
$\textit{EARSM}-\omega _0$
implementations. In both M6Tw025 and M14Tw018, adding
$\mathcal{K}$
and the new
$\textit{Pr}_t$
lowers the predicted temperature maximum while leaving the corresponding velocity profile close to its modified momentum-model baseline. The outer-layer velocity discrepancy inherited from the
$\omega _0$
reformulation remains. These calculations show that the algebraic closure can be inserted into the two tested implementations to improve temperature prediction without retuning their remaining model constants.
Effect of energy modelling closure within
$k-\omega$
family baselines in the M6Tw025 case. (a) Velocity profile, (b) non-dimensional temperature profile in wall units. The DNS data are from Zhang et al. (Reference Zhang, Duan and Choudhari2018).
Effect of energy modelling closure within
$k-\omega$
family baselines in the M14Tw018 case. (a) Velocity profile, (b) non-dimensional temperature profile in wall units. The DNS data are from Zhang et al. (Reference Zhang, Duan and Choudhari2018).
Figure 17. Long description
Two line graphs compare velocity and temperature profiles in wall units. Panel A: The line graph shows the velocity profile. The x-axis is labeled y+ on a logarithmic scale, and the y-axis is labeled u+ on a linear scale. The graph includes multiple data series: DNS (black circles), SST - ω0, k, Prt = 0.9 (blue dashed line), SST - ω0, k, new Prt (blue solid line), EARSM - ω0, k, Prt = 0.9 (red dashed line), and EARSM - ω0, k, new Prt (red solid line). The data series show a similar trend, increasing rapidly at lower y+ values and leveling off at higher y+ values. Panel B: The line graph shows the non-dimensional temperature profile in wall units. The x-axis is labeled y+ on a logarithmic scale, and the y-axis is labeled (T - T∞)/(Tw - T∞) on a linear scale. The graph includes the same data series as Panel A. The data series show a peak around y+ = 10 and then decrease, with some variations among the different models.
7. Summary and conclusions
This study develops a local algebraic closure for two modelled turbulent contributions to the Reynolds-averaged energy equation for attached compressible wall turbulence. The first component is an algebraic RANS-compatible surrogate
$\mathcal K(k,\nu _t^+)$
for the physical kinetic-energy terms in the momentum and total-energy equations. The transported variable
$k$
remains in the
$k-\omega _0$
transport equations and in the eddy-viscosity relation, whereas
$\mathcal{K}$
is used in the physical kinetic-energy terms of the momentum and total-energy equations, including the isotropic Reynolds-stress contribution and the kinetic-energy storage, convection and diffusion terms. The mapping is calibrated against a single incompressible Lag–EB boundary-layer calculation and is subsequently tested in cooled compressible boundary layers. It is neither a transported
$k-\epsilon$
equation nor a direct reconstruction of DNS TKE.
The DNS analysis of M6Tw025 identifies a discrepancy in the baseline kinetic-energy transport near the off-wall location where the turbulent heat flux passes through zero. Consistent with this diagnosis, replacing
$k$
by
$\mathcal{K}$
in the physical kinetic-energy terms accounts for most of the reduction in the overpredicted temperature maximum. In M6Tw025, the reported non-dimensional temperature maximum decreases from 2.336 with the baseline model to 2.236 with
$\mathcal{K}$
, compared with the DNS value of 2.206. In M14Tw018, the corresponding value decreases from 1.976 to 1.885, compared with the DNS value of 1.849. The effect is concentrated between the wall and the temperature maximum, whereas discrepancies remain in the wake region, where the turbulent heat flux provides the dominant turbulent transport contribution.
The second component retains the scalar gradient-diffusion approximation for the turbulent heat flux and introduces a local
$\textit{Pr}_t(\nu _t^+)$
relation. Regularised REnKF inversions of two
$M_b=1.5$
channel-flow cases are used to construct this relation; these cases provide the inversion data and are not treated as evidence of hypersonic performance. In the channel-flow tests, the proposed relation gives among the smaller centreline-temperature errors for the tested
$M_b=3.0$
cases at moderate
$\textit{Re}_\tau ^*$
. In the cooled boundary layers, however, its incremental effect is smaller than that of the
$\mathcal{K}$
reconstruction and depends on the observable. For M6Tw025, adding the new
$\textit{Pr}_t$
relation changes the temperature maximum from 2.236 to 2.204, compared with the DNS value of 2.206, and changes
$-B_q$
from 0.139 to 0.140, compared with the DNS value of 0.140. For M14Tw018, it changes the temperature maximum from 1.885 to 1.849, equal to the reported DNS value, but changes
$-B_q$
from 0.200 to 0.210, whereas the DNS value is 0.187. The
$\textit{Pr}_t$
relation therefore provides a secondary and metric-dependent correction and does not uniformly improve the wall heat flux.
Across the additional ZPG boundary layers, the complete closure reduces the temperature-maximum error relative to the
$k,\,\textit{Pr}_t=0.9$
baseline in M6Tw076, M8Tw048, M11Tw020 and M14Tw018, while producing little change in M2p5. The corresponding velocity profiles remain close to those of the momentum-model baseline. The two wall-normal resolutions tested for M6Tw025 and M14Tw018 give similar temperature profiles and wall quantities, indicating limited sensitivity to the particular refinement examined rather than general mesh independence. Calculations with the specifically defined
$\textit{SST}-\omega _0$
and
$\textit{EARSM}-\omega _0$
implementations show similar temperature trends for M6Tw025 and M14Tw018 without retuning the tested model constants. These results demonstrate compatibility with the two implementations examined, but do not establish transferability to all
$k-\omega$
-family models.
The demonstrated scope is limited to the attached channel flows and ZPG boundary layers examined here under the calorically perfect-gas and Sutherland-viscosity assumptions of the reference DNS calculations. In particular, the Mach-14 comparison is a validation against a controlled perfect-gas benchmark and does not assess high-enthalpy thermochemical effects. The present evidence does not cover pressure gradients, curvature, separation, shock–boundary-layer interaction, misalignment between the turbulent enthalpy flux and the mean enthalpy gradient or vibrational and thermochemical non-equilibrium. The closure introduces no additional transport equation, but its accuracy and numerical behaviour outside the tested flow configurations remain to be established.
Acknowledgements
The numerical computations were performed using Siyuan at the Centre for High-Performance Computing, Shanghai Jiao Tong University.
Funding
The authors gratefully acknowledge financial support from the National Natural Science Foundation of China (Grant No. 12472226).
Declaration of interests
The authors report no conflict of interest.
Appendix A. The
$\textit{Pr}_t$
models based on local variables
Kays & Crawford (Reference Kays and Crawford1993) proposed a universal dynamic model applicable to different incompressible fluids based on
$\nu _t^+$
and
$\textit{Pr}$
of the fluid. For air, the model equation is written as
\begin{equation} \textit{Pr}_t = \left [0.5882 + 0.2280\nu _t^+-0.04414{\nu _t^+}^2\left\{1-\mathrm{exp}\left(\frac {-5.165}{\nu _t^+}\right)\right\}\right ]^{-1} \!. \end{equation}
This model is referred to as Kays.
Ceci et al. (Reference Ceci, Pirozzoli and Spalart2026) constructed a
$\tilde {\nu }$
-based
$\textit{Pr}_t$
model for the Spalart–Allmaras (SA) model and applied it in RANS simulation of incompressible wall-bounded flows. This model follows the framework of eddy-viscosity damping function in the SA model and is calibrated by DNS simulation results of pipe flow reported by Pirozzoli & Modesti (Reference Pirozzoli and Modesti2024), written as
\begin{equation} \textit{Pr}_t = {\textit{Pr}_{t,\infty }}^3 \chi \frac {\chi _{\theta }^3+c_{v,\theta }^2}{\chi ^3 +c_{v1}^3}, \end{equation}
where
$\chi =\tilde {\nu }/\nu$
,
$\chi _{\theta }=\chi /\textit{Pr}_{t,\infty }$
,
$c_{v,\theta }=10$
and
$c_{v1}=7.1$
. For air,
$\textit{Pr}_{t,\infty }$
is
$0.843$
. To remove the dependence of (A2) on the SA model, an interpolation table that reflects
$\textit{Pr}_t$
vs
$\nu _t^+$
is constructed from an incompressible-channel-flow case at
$\textit{Re}_{\tau }=1000$
. This interpolation table is referred to as Ceci.
Huang et al. (Reference Huang, Coleman, Spalart and Yang2023) claimed that for compressible flows,
$\textit{Pr}_t$
can be written as a function of semi-local non-dimensional wall distance
$y^*$
, written as
\begin{equation} \textit{Pr}_t = 1.05 - 0.2 \mathrm{tanh}^3\left (\frac {y^*}{A^*_{Pr}}\right )\!, \end{equation}
where
$A^*_{\textit{Pr}}=70$
.
$\nu _t^+$
also follows the semi-local wall scaling, written as
\begin{equation} \nu _t^+ = \kappa y^* \left [1-\mathrm{exp}\left(-\frac {y^*}{A^*}\right)\right ]^2 \!, \end{equation}
where
$\kappa \approx 0.41$
and
$A^*=17$
. Therefore, by combining the inverse function of (A4) with (A3), a new local
$\textit{Pr}_t$
model can be constructed. In the following, this model is referred to as inv-Huang.
In compressible TBL, the analytical solution of
$\textit{Pr}_t$
is given as
\begin{equation} \textit{Pr}_t = \frac {\widetilde {\rho {u_1}''{u_2}''} \frac {\partial \tilde {h}}{\partial y}}{\widetilde {\rho h''{u_2}''}\frac {\partial \tilde {u_1}}{\partial y}}. \end{equation}
Following the construction methodology of the Kays model, a similar ad hoc fitting can be built based on existing DNS simulation results of compressible-channel flows reported by Lusher & Coleman (Reference Lusher and Coleman2022), Alcántara-Ávila & Hoyas (Reference Alcántara-Ávila and Hoyas2021) and Gerolymos & Vallet (Reference Gerolymos and Vallet2023) as a function of
$\nu _t^+$
, shown in figure 18. The fitted model is denoted as channel-DNS-fitted profile. The explicit fitting formulae are written as (A6)
\begin{equation} \begin{cases} \textit{Pr}_{t}^{vis} = 0.98 - 0.04 \ln (\nu _t^+), \\[9pt] \textit{Pr}_t^{log} = 1.05 - 0.2 \tanh ^3 \left ( \frac {\nu _t^+}{30} \right ) \!,\\[9pt] \phi = \tanh \left ( \frac {\nu _t^+}{12} \right ) \!,\\[9pt] \textit{Pr}_t = \textit{Pr}_t^{vis}(1 - \phi ) + \textit{Pr}_t^{log} \phi .\end{cases} \end{equation}
Visualisation of channel-DNS-fitted profile. The black line denotes the fitting function, and the grey lines denote
$\textit{Pr}_t$
profiles of the DNS simulation. The reference cases are aiD, aiD2, aiE and aiF2 reported by Lusher & Coleman (Reference Lusher and Coleman2022), Retau500, Retau1000, Retau2000 and Retau5000 reported by Alcántara-Ávila & Hoyas (Reference Alcántara-Ávila and Hoyas2021) and Retau965 and Retau985 reported by Gerolymos & Vallet (Reference Gerolymos and Vallet2023).
A comparison of the proposed new formulation with existing incompressible and compressible turbulent-Prandtl-number models is shown in figure 19.
Comparison of the new formulation with the Kays (Kays & Crawford Reference Kays and Crawford1993), Ceci (Ceci et al. Reference Ceci, Pirozzoli and Spalart2026), Huang (Huang et al. Reference Huang, Coleman, Spalart and Yang2023) and DNS-fitted correlations as functions of
$\nu _t^+$
. Profiles are given by an incompressible-channel-flow case at
$\textit{Re}_{\tau }=5000$
. Solid lines indicate that the models are calibrated using incompressible flow cases, while dashed lines indicate that the models are calibrated using compressible ones.
Reτ
Reτ∗
Reb
Mab
Maτ
T~c/Tw
Reτ
Reτ∗
Mab
T~c/Tw
Nx
Ny
k
K
Reθ≈2500.
k−ϵ
νt+
k
K
u+=u/uτ
T+=T/Tw
Prt
ϕPrt
Prt
Prt
Mb=0.7
Mb=1.7
Mb=3.0
Prt
Prt
y+
Reτ
Reτ∗
uτ
Bq
Prt
Prt
Reθ
Reτ
Reτ∗
θ
δ
uτ
Bq
Bq
k−ω
k−ω
Prt
νt+
Reτ=5000
